Changeset 1196

Timestamp:

06/15/07 19:09:54 (1 year ago)

Author:

rporter

Message:

More iso refactoring. Still no botches. Biggest changes are as follows:

helike_level.cpp and hydrolevelpop.cpp have been deleted and replaced with new iso_level.cpp

In quantumState structure, spin variable s has been replaced by S, often used to denote multiplicity. This way the number will always be an integer, which, of course, is not the case for spin. Many of the files only changed as a reflection of this change. In some places, 2*s+1 were changed to just S.

prt_He_like_DeparCoef and prt_He_like_Pops were merged and renamed to iso_prt_pops, which now takes a third flag that causes departure coefficients to be printed if true. H-like analogs prt_H_like_DeparCoef and prt_H_like_Pops are now only sticking around for ease of comparison with old output. Will go away when new treatment is fully functional. Won't be long now!

Files:

• branches/porter/source/assert_results.cpp (modified) (1 diff)
• branches/porter/source/cddefines.h (modified) (1 diff)
• branches/porter/source/conv_base.cpp (modified) (2 diffs)
• branches/porter/source/helike.cpp (modified) (10 diffs)
• branches/porter/source/helike.h (modified) (4 diffs)
• branches/porter/source/helike_create.cpp (modified) (22 diffs)
• branches/porter/source/helike_cs.cpp (modified) (8 diffs)
• branches/porter/source/helike_cs.h (modified) (3 diffs)
• branches/porter/source/helike_einsta.cpp (modified) (16 diffs)
• branches/porter/source/helike_level.cpp (deleted)
• branches/porter/source/helike_recom.cpp (modified) (7 diffs)
• branches/porter/source/hydrogenic.cpp (modified) (4 diffs)
• branches/porter/source/hydrogenic.h (modified) (1 diff)
• branches/porter/source/hydrolevel.cpp (modified) (5 diffs)
• branches/porter/source/hydrolevelpop.cpp (deleted)
• branches/porter/source/init_coreload.cpp (modified) (1 diff)
• branches/porter/source/ion_recomb_Badnell.cpp (modified) (4 diffs)
• branches/porter/source/iso.h (modified) (5 diffs)
• branches/porter/source/iso_create.cpp (modified) (23 diffs)
• branches/porter/source/iso_level.cpp (added)
• branches/porter/source/iso_radiative_recomb.cpp (modified) (1 diff)
• branches/porter/source/prt_lines.cpp (modified) (1 diff)
• branches/porter/source/prt_lines_helium.cpp (modified) (3 diffs)
• branches/porter/source/prt_zone.cpp (modified) (3 diffs)
• branches/porter/source/rt_tau_init.cpp (modified) (1 diff)
• branches/porter/source/rt_tau_reset.cpp (modified) (1 diff)
• branches/porter/source/sanity_check.cpp (modified) (1 diff)

branches/porter/source/assert_results.cpp

@@ -2380 +2380 @@
-                                                                                StatesElem[ipHE_LIKE][ipHELIUM][ipHi].n,
-                                                                                StatesElem[ipHE_LIKE][ipHELIUM][ipHi].l,
-s
+S
-                                                                                StatesElem[ipHE_LIKE][ipHELIUM][ipLo].n,
-                                                                                StatesElem[ipHE_LIKE][ipHELIUM][ipLo].l,
-s
+S
-                                                                        break;
-                                                                }

branches/porter/source/cddefines.h

@@ -1031 +1031 @@
-        double ConBoltz;
-
-
+
+
-
-        double lifetime;

branches/porter/source/conv_base.cpp

@@ -423 +423 @@
-                /* do all hydrogenic species, and fully do hydrogen itself, including molecules */
-                /*fprintf(ioQQQ,"DEBUG h2\t%.2f\t%.3e\t%.3e", fnzone,hmi.H2_total,hmi.Hmolec[ipMH2g]);*/
-
+ ipH_LIKE 
-                /*fprintf(ioQQQ,"\t%.3e\n", hmi.Hmolec[ipMH2g]);*/
-
@@ -437 +437 @@
-                /* >>chng 02 mar 08 move helike to before helium */
-                /* do all he-like species */
-
+ ipHE_LIKE 
-
-                /*>>chng 04 may 09, add option to abort here, inspired by H2 pop failures

branches/porter/source/helike.cpp

@@ -1 +1 @@
-/* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland
- * For conditions of distribution and use see copyright notice in license.txt */
-
-
-
+
+
-/*AGN_He1_CS routine to punch table needed for AGN3 - collision strengths of HeI */
-/*defect - calculate quantum defect. */
-/*helike_energy - calculates energy of a given level.   */
-
-#include "cddefines.h" 
-
-/* various useful print statements occur throughout this file - search on PRINTIT */
-/* 
-  
-  Energy order within 2 3P
-
-  The order of the levels within the 2 3P level of atomic helium is opposite
-  from the order within astrophysically abundant ions.  The indices below
-  consistently point to the correct level, and the energies are correct,
-  so the J levels within 2 3P are not in increasing energy order for helium itself. 
-  This is ok since the atomic data is correct, and the difference in energies is so small.
-
-*/
-
-/* 
-the following block of defines occurs in cddefines.h 
-const int ipHe1s1S = 0;
-const int ipHe2s3S = 1;
-const int ipHe2s1S = 2;
-const int ipHe2p3P0 = 3;
-const int ipHe2p3P1 = 4;
-const int ipHe2p3P2 = 5;
-const int ipHe2p1P = 6;
-
-level 3
-const int ipHe3s3S = 7;
-const int ipHe3s1S = 8;
-const int ipHe3p3P = 9;
-const int ipHe3d3D = 10;
-const int ipHe3d1D = 11;
-const int ipHe3p1P = 12;
-
-level 4
-const int ipHe4s3S = 13;
-const int ipHe4s1S = 14;
-const int ipHe4p3P = 15;
-const int ipHe4d3D = 16;
-const int ipHe4d1D = 17;
-const int ipHe4f3F = 18;
-const int ipHe4f1F = 19;
-const int ipHe4p1P = 20;
-*/
-
-/*const int ipRad = 0;
-const int ipCollis = 1;
-const int ipEnergy = 2;*/
-
-
-#include "dense.h"
-#include "elementnames.h"
@@ -62 +10 @@
-#include "helike_cs.h"
-#include "helike_recom.h"
+#include "hydrogenic.h"
-#include "iso.h"
-#include "phycon.h"
@@ -71 +20 @@
-/*lint -e661 creation of  out of bound pointer */
-
-/*static long ipLev,globalZ;*/
-
-/* Ionization potentials (in eV) for each ion in the iso-seq up to Z=30. These are 
- * exactly what you get if you take Verner's numbers in atmdat_ph1(0,1,nelem,0) 
- * and multiply by 0.9998787, exactly as is done elsewhere in Cloudy.   */
-double EionEV[29] =       
-        {24.59077,75.630824,153.881326,259.368529,392.052444,552.033006,739.210311,
-        953.784316,1195.854925,1464.822296,1761.786269,2085.746968,2437.704271,
-        2816.658298,3223.608929,3658.556163,4120.500123,4610.440686,5128.377852,
-        5674.311623,6248.241996,6851.168852,7481.092433,8140.012497,8827.929042,
-        9543.842191,10288.751823,11058.658422,11868.560169};
-
-/*He-like main routine to call HeLevel and determine model he-like species atom level balance,
- * called by ionize */
-void
+long ipISO
-{
-
-
-
-
+
+
+
+
+
-        static bool lgFinitePop[LIMELM];
-        static bool lgMustInit=true;
@@ -134 +72 @@
-
-                                /* evaluate collisional rates */
-
+
+
+
+
+
+
-
-                                /* evaluate photoionization rates */
@@ -147 +90 @@
-
-                                /* solve for the level populations */
-HeLikeL
+iso_l
-
-                                /* say that we have set the populations */
@@ -154 +97 @@
-                                if( trace.lgTrace )
-                                {
-Helium-like Z=%2ld H2OVH1="
+iso=%2ld Z=%2ld H2OVH1=",ipISO
-                                        fprintf( ioQQQ,PrintEfmt("%10.3e", iso.pop_ion_ov_neut[ipISO][nelem]));
-                                        fprintf( ioQQQ, " simple=");
@@ -430 +373 @@
-}
-
-He_like_DeparCoef print departure coefficients for he-like specie
-prt_He_like_DeparCoef( long ipISO, long nelem
+iso_Pops print out iso sequence populations or departure coefficient
+iso_prt_pops( long ipISO, long nelem, long lgPrtDeparCoef
-{
-        long int in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-                ipFirstCollapsed= iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem];
@@ -512 +410 @@
-        
-                /* give element number and spin */
-population
+%
-                        elementnames.chElementSym[nelem],
-
+
+
-
-                /* header with the l states */
@@ -525 +424 @@
-
-                /* loop over prin quant numbers, one per line, with l across */
-
-
+
+
+
+
+
-                        fprintf(ioQQQ," %2ld           ",in);
-
-                        for( il = 0; il < in; ++il)
-                        {
-(in==2) && (il==1) && (is==1
-
-StatesElem[ipISO][nelem][ipHe2p3P0].Pop
-StatesElem[ipISO][nelem][ipHe2p3P1].Pop
-StatesElem[ipISO][nelem][ipHe2p3P2].Pop
+ipISO==ipHE_LIKE && (in==2) && (il==1) && (is==3
+
+ITEM_TO_PRINT(ipHe2p3P0) 
+ITEM_TO_PRINT(ipHe2p3P1) 
+ITEM_TO_PRINT(ipHe2p3P2) 
-                                }
-                                else
-                                {
-                                        ipLo = iso.QuantumNumbers2Index[ipISO][nelem][in][il][is];
-StatesElem[ipISO][nelem][ipLo].Pop
+ITEM_TO_PRINT(ipLo) 
-                                }
-                        }
@@ -552 +454 @@
-                /* prin quan number of collapsed levels */
-                fprintf(ioQQQ," %2ld           ",in);
-StatesElem[ipISO][nelem][il].Pop
+ITEM_TO_PRINT(il) 
-                fprintf(ioQQQ,"\n");
-        }
-
-prt_He_like_P
+iso_prt_p
-        return;
-}

branches/porter/source/helike.h

@@ -22 +22 @@
-#define CASEABMAGIC             (51214)
-
-
-
-
+
+
+
+
+
-
-LikeLevel level populations  
+Create create he-like series
-\param ipISO
-\param nelem
-*/
-void HeLikeLevel( long int ipISO, long int nelem);
-
-/** create he-like series */
-void HeCreate( long ipISO );
-
-
+HeCollidSetup
-void HeCollidSetup( void );
-
-
+helike_energy
-\param nelem
-\param ipLev 
@@ -44 +42 @@
-double helike_energy(long nelem, long ipLev );
-
- get helike quantum defect 
+helike_quantum_defect get quantum defect for helium-like level
-\param nelem
-\param ipLev
@@ -50 +48 @@
-double helike_quantum_defect( long int nelem, long int ipLev );
-
- get helike transition probability in s-1 
+helike_transprob get transition probability for helium-like transition [s-1]
-\param nelem
-\param ipHi
@@ -57 +55 @@
-float helike_transprob( long nelem, long ipHi, long ipLo );
-
-
-
-
-
-
-
-
-
-
-
-
+
-\param *ioPun
-*/
-void AGN_He1_CS( FILE *ioPun );
-
-
-
+
+
+
-

branches/porter/source/helike_create.cpp

@@ -1 +1 @@
-/* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland
- * For conditions of distribution and use see copyright notice in license.txt */
-
+
+
-#include "physconst.h" 
-#include "lines_service.h"
@@ -17 +18 @@
-#include "helike_einsta.h"
-
+/*
+  Energy order within 2 3P
+
+  The order of the levels within the 2 3P level of atomic helium is opposite
+  from the order within astrophysically abundant ions.  The indices below
+  consistently point to the correct level, and the energies are correct,
+  so the J levels within 2 3P are not in increasing energy order for helium itself. 
+  This is ok since the atomic data is correct, and the difference in energies is so small.
+*/
+
+/* 
+the following block of defines occurs in cddefines.h 
+const int ipHe1s1S = 0;
+const int ipHe2s3S = 1;
+const int ipHe2s1S = 2;
+const int ipHe2p3P0 = 3;
+const int ipHe2p3P1 = 4;
+const int ipHe2p3P2 = 5;
+const int ipHe2p1P = 6;
+
+level 3
+const int ipHe3s3S = 7;
+const int ipHe3s1S = 8;
+const int ipHe3p3P = 9;
+const int ipHe3d3D = 10;
+const int ipHe3d1D = 11;
+const int ipHe3p1P = 12;
+
+level 4
+const int ipHe4s3S = 13;
+const int ipHe4s1S = 14;
+const int ipHe4p3P = 15;
+const int ipHe4d3D = 16;
+const int ipHe4d1D = 17;
+const int ipHe4f3F = 18;
+const int ipHe4f1F = 19;
+const int ipHe4p1P = 20;
+*/
+
-/* the designation of the levels, chLevel[n][string] */
-static char **chLevel;
-
-               */
-/
+
+ 
-static double EionWN[LIMELM] =       
-        /* the first, or H-like, element is not defined for he-like species */
- 
- 
- 
- 
- 
- 
- 
+
+
+
+
+
+
+
-        
-*/
-/
+
+ 
-static double EionRYD[LIMELM] =       
-        /* the first, or H-like, element is not defined for he-like species */
@@ -42 +82 @@
-        549.850208,598.279945,648.840883,701.459535,756.209388,812.796486,872.323172};
-
+/* Ionization potentials (in eV) for each ion in the iso-seq up to Z=30. These are 
+ * exactly what you get if you take Verner's numbers in atmdat_ph1(0,1,nelem,0) 
+ * and multiply by 0.9998787, exactly as is done elsewhere in Cloudy.   */
+double EionEV[LIMELM] =       
+        /* the first, or H-like, element is not defined for he-like species */
+        {-DBL_MAX,
+        24.59077,75.630824,153.881326,259.368529,392.052444,552.033006,739.210311,
+        953.784316,1195.854925,1464.822296,1761.786269,2085.746968,2437.704271,
+        2816.658298,3223.608929,3658.556163,4120.500123,4610.440686,5128.377852,
+        5674.311623,6248.241996,6851.168852,7481.092433,8140.012497,8827.929042,
+        9543.842191,10288.751823,11058.658422,11868.560169};
-
-/* experimental energies, in wavenumbers, for atomic helium */
@@ -297 +348 @@
-        long int n = StatesElem[ipHE_LIKE][nelem][ipLo].n;
-        long int lqn = StatesElem[ipHE_LIKE][nelem][ipLo].l;
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-        DEBUG_ENTRY( "helike_quantum_defect()" );
@@ -410 +474 @@
-        {
-                Eff_n = N_(ipLev) - helike_quantum_defect( nelem, ipLev );
-
+
+
-
-                /* energies (in wavenumbers) that correspond to quantum defect */
@@ -430 +495 @@
-        double **SumAPerN, **RobbinsC;
-
-ipFirstCollapsed, 
+
-
-        static int nCalled = 0;
-
-        char chSpin[2]={'1','3'};
-        char chL[6]={'S','P','D','F','G','H'};
-
@@ -454 +518 @@
-        }
-
- first find the largest number of levels in any element in the iso sequence */
-max_num_levels = 0;
+* find the largest n ever explicitly considered in the he-like iso sequence 
+ * resolved or collapsed is irrelevant.  */
-        max_n = 0;
-        for( nelem=ipISO; nelem < LIMELM; nelem++ )
@@ -462 +526 @@
-                if( nelem == ipISO || dense.lgElmtOn[nelem] )
-                {
-                        max_num_levels = MAX2( max_num_levels, iso.numLevels_max[ipISO][nelem] );
-                        max_n = MAX2( max_n, iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] );
-                }
@@ -509 +572 @@
-                        }
-
-                        /* this is the number of resolved levels, so first collapsed level is [ipFirstCollapsed] */
-                        ipFirstCollapsed = iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; 
-                        
-                        /* Initialize some ground state stuff, easier here than in loops.       */
-                        RobbinsC[ipHe1s1S][ipHe1s1S] = 1.;
-                        StatesElem[ipISO][nelem][ipHe1s1S].lifetime = 0.;
-                        iso.CascadeProb[ipISO][nelem][ipHe1s1S][ipHe1s1S] = 1.;
-                        if( iso.lgRandErrGen[ipISO] )
@@ -522 +581 @@
-                        }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
-                        for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
-                        {
@@ -570 +589 @@
-                                RobbinsC[ipHi][ipHi] = 1.;
-                                iso.CascadeProb[ipISO][nelem][ipHi][ipHi] = 1.;
+                                iso.BranchRatio[ipISO][nelem][ipHi][ipHe1s1S] = 0.;
+
-                                if( iso.lgRandErrGen[ipISO] )
-                                {
@@ -575 +596 @@
-                                        iso.SigmaCascadeProb[ipISO][nelem][ipHi][ipHi] = 0.;
-                                }
-
+
-                                for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
-                                {
-                                        RobbinsC[ipHi][ipLo] = 0.;
-                                        iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.;
+                                        iso.BranchRatio[ipISO][nelem][ipHi][ipLo] = 
+                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis.Aul *
+                                                StatesElem[ipISO][nelem][ipHi].lifetime;
+                                        /* >>chng 06 mar 08, had been <= 1.0000001, increased to
+                                         * when Aul converted to float */
+                                        /* \todo 2 find out what's wrong here. */
+                                        ASSERT( iso.BranchRatio[ipISO][nelem][ipHi][ipLo] <= 1.007 );
+
+                                        SumAPerN[nelem][N_(ipHi)] += Transitions[ipISO][nelem][ipHi][ipLo].Emis.Aul;
-
-                                        /* there are some negative energy transitions, where the order
@@ -587 +617 @@
-                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis.Aul <= iso.SmallA );
-
-
-
-
-
+
-                                        {
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
-                                        }
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-                                        {
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
-                                                
-                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis.gf = 
-                                                        (float)(GetGF(Transitions[ipISO][nelem][ipHi][ipLo].Emis.Aul,
-                                                        Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN,
-                                                        Transitions[ipISO][nelem][ipHi][ipLo].Hi.g));
-                                                ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis.gf > 0.);
-
-                                                /* derive the abs coef, call to function is Emis.gf, wl (A), g_low */
-                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis.opacity = 
-                                                        (float)(abscf(Transitions[ipISO][nelem][ipHi][ipLo].Emis.gf,
-                                                        Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN,
-                                                        Transitions[ipISO][nelem][ipHi][ipLo].Lo.g));
-                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis.opacity = 
-                                                        MAX2(0.f, Transitions[ipISO][nelem][ipHi][ipLo].Emis.opacity );
-                                        }
-
+
+
-                                        {
-
-