Changeset 2376 for branches/newmole

Timestamp:

09/30/08 16:36:15 (3 months ago)

Author:

rjrw

Message:

Merged from trunk r2345:2375

Location:

branches/newmole/source

Files:

• atmdat_adfa.cpp (modified) (1 diff)
• atmdat_ligbar.cpp (modified) (2 diffs)
• atom_feii.cpp (modified) (5 diffs)
• cddefines.h (modified) (2 diffs)
• container_classes.h (modified) (1 diff)
• cont_createmesh.cpp (modified) (1 diff)
• cool_iron.cpp (modified) (1 diff)
• date.h (modified) (1 diff)
• ionbal.h (modified) (2 diffs)
• ion_recomb.cpp (modified) (4 diffs)
• ion_recomb_Badnell.cpp (modified) (10 diffs)
• iter_startend.cpp (modified) (2 diffs)
• parse_grain.cpp (modified) (1 diff)
• parse_punch.cpp (modified) (1 diff)
• parse_set.cpp (modified) (7 diffs)
• punch_do.cpp (modified) (1 diff)
• radius_increment.cpp (modified) (4 diffs)
• rfield.h (modified) (2 diffs)
• rt_diffuse.cpp (modified) (12 diffs)
• rt_line_one.cpp (modified) (4 diffs)
• stars.cpp (modified) (2 diffs)
• temp_change.cpp (modified) (8 diffs)
• zero.cpp (modified) (3 diffs)

branches/newmole/source/atmdat_adfa.cpp

@@ -863 +863 @@
         ASSERT( ipLo < ipHi );
 
-#ifndef NDEBUG
+#       if !defined NDEBUG || defined ASSERTDEBUG
         long ipISO = ipH_LIKE;
         long nelem = ipHYDROGEN;
-#endif
+#       endif
         ASSERT( N_(ipLo) < N_(ipHi) );
         ASSERT( N_(ipHi) <= 5 );

branches/newmole/source/atmdat_ligbar.cpp

@@ -130 +130 @@
         /* tarray(ipLnGF) = gf; tarray(ipLnBolt) excit temp kelvin */
         /*
-        *cs2s2p = gbar*197.47*EVDEGK*t2s2p->Lo->gf/t2s2p->EnergyK;
-        */
-         *cs2s2p = gbar*197.47*EVDEGK*t2s2p->Emis->gf/t2s2p->EnergyK;
+         *cs2s2p = gbar*197.47*EVDEGK*t2s2p->Lo->gf/t2s2p->EnergyK;
+         */
+        *cs2s2p = gbar*197.47*EVDEGK*t2s2p->Emis->gf/t2s2p->EnergyK;
         /* small correction factors to CMcW83 2s-2p fits:
          * fits 3.57% too small compared to R-matrix calc. for Mg X.
@@ -167 +167 @@
         /* tarray(ipLnGF) = gf */
         /**cs2s3p = gbar*197.47*EVDEGK*t2s3p->Lo->gf/t2s3p->EnergyK;*/
-     *cs2s3p = gbar*197.47*EVDEGK*t2s3p->Emis->gf/t2s3p->EnergyK;
+        *cs2s3p = gbar*197.47*EVDEGK*t2s3p->Emis->gf/t2s3p->EnergyK;
         /* cs2s3p = gbar * 197.47*eVdegK *  GF2/(e2/1.60184e-12)
          * */

branches/newmole/source/atom_feii.cpp

@@ -480 +480 @@
         fclose(ioDATA);
 
-        /* now get radiation data.
-         * These are from the following data sources:
-         >>refer        fe2     as      Nahar, S., 1995, A&A 293, 967
-         >>refer        fe2     as      Quinet, P., LeDourneuf, M., & Zeipppen C.J., 1996, A&AS, 120, 361 
-         >>refer        fe2     as      Furh, J.R., Martin, G.A., & Wiese, W.L., 1988; J Phys Chem Ref Data 17, Suppl 4
-         >>refer        fe2     as      Giridhar, S., & Arellano Ferro, A., 1995; Ref Mexicana Astron Astrofis 31, 23
-         >>refer        fe2     as      Kurucz, R.L., 1995, SAO CD ROM 23
-         */
+        /* transition probabilities */
 
         if( trace.lgTrace )
-        {
                 fprintf( ioQQQ," FeIICreate opening fe2rad.dat:");
-        }
 
         ioDATA = open_data( "fe2rad.dat", "r" );
@@ -504 +495 @@
         /* scan off three integers */
         sscanf( chLine ,"%ld%ld%ld",&lo, &ihi, &m1);
-        if( lo!=3 || ihi!=1 || m1!=28 )
-        {
-                fprintf( ioQQQ, " fe2rad.dat has the wrong magic numbers, expected 3 1 28\n" );
+        const int nYrFe2Rad=8, nMoFe2Rad=8, nDyFe2Rad=24;
+        if( lo!=nYrFe2Rad || ihi!=nMoFe2Rad || m1!=nDyFe2Rad )
+        {
+                fprintf( ioQQQ, "DISASTER fe2rad.dat has the wrong magic numbers, expected "
+                        "%2i %2i %2i and got %2li %2li %2li\n" ,
+                        nYrFe2Rad, nMoFe2Rad, nDyFe2Rad,
+                        lo, ihi, m1);
                 cdEXIT(EXIT_FAILURE);
         }
@@ -534 +529 @@
                                 &lo, &ihi, &m1, &m2 ,
                                 &save[0], &save[1] , &m3);
-                        /* the pointeres ihi and lo within this array were 
+                        /* the indices ihi and lo within this array were 
                          * read in from the line above */
                         Fe2LevN[ihi-1][lo-1].Lo->g = (realnum)m1;
@@ -545 +540 @@
                         else
                                 Fe2LevN[ihi-1][lo-1].EnergyWN = Fe2Energies[ihi-1]-Fe2Energies[lo-1];
-                        if( Fe2LevN[ihi-1][lo-1].Emis->Aul == 1e3f )
-                        {
-                                /*realnum xmicron;
-                                xmicron = 1e4f/ Fe2LevN[ihi-1][lo-1].EnergyWN;*/
-                                /* these are made-up intercombination lines, set gf to 1e-5
-                                 * then get better A */
+
+                        /* Aul == -1 indicates intercombination line with no real Aul */
+                        if( fp_equal( Fe2LevN[ihi-1][lo-1].Emis->Aul , -1.0f ) )
+                        {
+                                /* these are made-up intercombination lines, set gf to 1e-5 */
                                 Fe2LevN[ihi-1][lo-1].Emis->gf = 1e-5f;
-
+                                
+                                /* get corresponding A */
                                 Fe2LevN[ihi-1][lo-1].Emis->Aul = Fe2LevN[ihi-1][lo-1].Emis->gf*(realnum)TRANS_PROB_CONST*
                                         POW2(Fe2LevN[ihi-1][lo-1].EnergyWN)*Fe2LevN[ihi-1][lo-1].Lo->g/Fe2LevN[ihi-1][lo-1].Hi->g;
-
-                                /*fprintf(ioQQQ," %e from 1e3 to %e\n",xmicron , Fe2LevN[ihi-1][lo-1].Emis->Aul );*/
-                        }
-                        /* this is the last column of fe2rad.dat, and is 1, 2, or 3.  
+                        }
+
+                        /* the last column of fe2rad.dat, and is 1, 2, or 3.  
                          * 1 signifies that transition is permitted,
                          * 2 is semi-forbidden
@@ -564 +558 @@
                          * transition is from 64 */
                         ncs1[ihi-1][lo-1] = (int)m3;
-                        /*if( m3==1 )
-                                fprintf(ioQQQ,"DEBUG fe2 new old energ\t%.5e\t%.5e\t%.3e\n",
-                                Fe2Energies[ihi-1]-Fe2Energies[lo-1],
-                                Fe2LevN[ihi-1][lo-1].EnergyWN ,
-                                (Fe2Energies[ihi-1]-Fe2Energies[lo-1]-Fe2LevN[ihi-1][lo-1].EnergyWN)/
-                                Fe2LevN[ihi-1][lo-1].EnergyWN
-                                );*/
                 }
         }

branches/newmole/source/cddefines.h

@@ -39 +39 @@
 #include <memory>
 #include <stdexcept>
+#include <algorithm>
 #ifdef DMALLOC
 #include <dmalloc.h>
@@ -595 +596 @@
 #undef  ASSERT
 #ifndef STD_ASSERT
-#       ifdef NDEBUG
+#       ifdef ASSERTDEBUG
+#               define ASSERT(exp) if (!(exp)) MyAssert(__FILE__, __LINE__)
+#       elif NDEBUG
 #               define ASSERT(exp) ((void)0)
-#       elif defined ASSERTDEBUG
-#               define ASSERT(exp) if (!(exp)) MyAssert(__FILE__, __LINE__)
 #       else
 /* the do { ... } while ( 0 ) prevents bugs in code like this:

branches/newmole/source/container_classes.h

@@ -66 +66 @@
         }
 public:
+        typedef random_access_iterator_tag iterator_category;
+        typedef T            value_type;
+        typedef T&           reference;
+        typedef const T&     const_reference;
+        typedef T*           pointer;
+        typedef const T*     const_pointer;
+        typedef size_t       size_type;
+        typedef ptrdiff_t    difference_type;
+
         // constructors
         basic_pntr( T* p0, T* p1, T* p2 )

branches/newmole/source/cont_createmesh.cpp

@@ -538 +538 @@
         rfield.SummedOcc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
         rfield.ConOTS_local_photons = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
+        rfield.DiffuseEscape = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
         rfield.TotDiff2Pht = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
         rfield.ConOTS_local_OTS_rate = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );

branches/newmole/source/cool_iron.cpp

@@ -327 +327 @@
                 {
                         /* spacing needed to get proper trace convergence output */
-                        fprintf( ioQQQ, "           FeIILevelPops5 returned cool=%.2e heat=%.2e derivative=%.2e\n ",
+                        fprintf( ioQQQ, "           FeIILevelPops5 returned cool=%.2e heat=%.2e derivative=%.2e\n",
                                         FeII.Fe2_large_cool,FeII.Fe2_large_heat ,FeII.ddT_Fe2_large_cool);
                 }

branches/newmole/source/date.h

@@ -14 +14 @@
 #define YEAR    108
 /* month, January is 0, December is 11 */
-#define MONTH   7
+#define MONTH   8
 /* day is correct */
 #define DAY     9

branches/newmole/source/ionbal.h

@@ -220 +220 @@
                 **lgRR_Badnell_rate_coef_exist;
 
-        /** use the mean Badnell rates in place of existing hacks? */
-        bool lg_use_DR_Badnell_rate_coef_mean_ion;
-
         /** do we use new Badnell rates? */
         bool lgDR_recom_Badnell_use,
@@ -267 +264 @@
         bool lgSupDie[2];
 
-        /** this is fudge factor for scaled Nussbaumer and Storey dielectronic 
-         recombination, set with dielectronic kludge command, usually 1 */
-        realnum GuessDiel[4];
-
         /** following all for 3-body recombination */
         /** lgNoCota flag set with no three body recombination */

branches/newmole/source/ion_recomb.cpp

@@ -43 +43 @@
           fac2, 
           factor,
-          DielRecomRateCoef_HiT[LIMELM] ,
-          DielRecomRateCoef_LowT[LIMELM] ,
-          ChargeTransfer[LIMELM] ,
-          t4m1, 
-          tefac;
+          DielRecomRateCoef_HiT[LIMELM],
+          DielRecomRateCoef_LowT[LIMELM],
+          ChargeTransfer[LIMELM];
 
         /* these are used for adding noise to rec coef */
@@ -72 +70 @@
 
         atmdat.nsbig = MAX2(dense.IonHigh[nelem]+1,atmdat.nsbig);
-        t4m1 = 1e4/phycon.te;
         fac2 = 1e-14*phycon.sqrte;
 
@@ -198 +195 @@
                 DielRecomRateCoef_LowT[ion] = 0.;
                 if( ((n_bnd_elec_after_recomb-1) !=  2) && 
-                        ((n_bnd_elec_after_recomb-1) != 10) && 
-                        ((n_bnd_elec_after_recomb-1) != 18) )
-                {
-                        tefac = ff[ion]*t4m1;
+                    ((n_bnd_elec_after_recomb-1) != 10) && 
+                    ((n_bnd_elec_after_recomb-1) != 18) )
+                {
                         /* do not do iron here since all dn = 0 DR either Badnell or a guess */
-                        if( ff[ion] != 0. && nelem!=ipIRON )
-                        {
+                        if( ff[ion] != 0. && nelem != ipIRON )
+                        {
+                                double t4m1 = 1e4/phycon.te;
+                                double tefac = ff[ion]*t4m1;
                                 /* >>chng 06 feb 14, O+3 ff[ion] is very negative, as a result exp goes to +inf
                                  * at very low temperatures.  This is a error in the Nussbaumer & Storey fits to DR.
@@ -221 +219 @@
                         else if( ionbal.lg_guess_coef )
                         {
-                                /* first guesses are those based on Nussbaumer & Storey and are from
-                                 * >>refer      ion     DR      Ali, B., Blum, R. D., Bumgardner, T. E., Cranmer, S. R., 
-                                 * >>refercon   Ferland, G. J., Haefner, R. I., & Tiede, G. P. 1991, PASP, 103, 1182*/
-                                if( (ff[ion] == 0.) && (ion <= 3)  )
-                                {
-                                        /* these are guessed dielectronic rec rates, taken from 
-                                         * >>refer      all     dielectronic    Ali, B., Blum, R. D., Bumgardner, T. E., Cranmer, S. R., Ferland, G. J., 
-                                         * >>refercon Haefner, R. I., & Tiede, G. P. 1991, PASP, 103, 1182 */
-                                        static double cludge[4]={3e-13,3e-12,1.5e-11,2.5e-11};
-
-                                        /* >>chng 01 jun 30, make GuessDiel array */
-                                        DielRecomRateCoef_LowT[ion] = ionbal.DielSupprs[1][ion]*cludge[ion]*
-                                                /* this is optional scale factor on set dielectronic rec command 
-                                                 * >>chng 06 feb 07, add Boltzmann factor to induce behavior
-                                                 * like in Nussbaumer & Storey */
-                                                ionbal.GuessDiel[ion] * sexp( 1e3 / phycon.te );
-                                }
-                                /* this implements the low-T kludge dielectronic command */
-                                else
-                                {
-                                        /* assume the recombination rate is the coefficient scanned off the steve option, times the charge
-                                         * before recombination raised to a power */
-                                        double fitcoef[3][3] = 
-                                        {
-                                                /* L- shell, Z=17-23 */
-                                                {-52.5073,+5.19385,-0.126099} ,
-                                                /*  M-shell (Z=9-15) */
-                                                {-10.9679,1.66397,-0.0605965} ,
-                                                /* N shell */
-                                                {-3.95599,1.61884,-0.194540} 
-                                        };
-
-                                        /* these implement the guesses in 
-                                         * Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ, 604, 556  */
-                                        if( nelem==ipIRON )
-                                        {
-                                                int nshell;
-                                                if( (n_bnd_elec_after_recomb>=4) && (n_bnd_elec_after_recomb<=11) )
-                                                        nshell = 0;
-                                                else if( (n_bnd_elec_after_recomb>=12) && (n_bnd_elec_after_recomb<=19 ))
-                                                        nshell = 1;
-                                                else
-                                                        nshell = 2;
-                                                /* n_bnd_elec_after_recomb is the number of bound electrons */
-                                                DielRecomRateCoef_LowT[ion] = fitcoef[nshell][0] +
-                                                        fitcoef[nshell][1]*(ion+1) + 
-                                                        fitcoef[nshell][2]*POW2(ion+1.);
-                                                DielRecomRateCoef_LowT[ion] = 1e-10*pow(10.,DielRecomRateCoef_LowT[ion]);
-
-                                        }
-                                        else
-                                        {
-                                                /* this is guess for all other elements presented in 
-                                                 * >>refer      all     DR      Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ, 604, 556 */
-                                                DielRecomRateCoef_LowT[ion] = 3e-12*pow(10.,(double)(ion+1)*0.1);
-                                        }
-                                        /* >>chng 06 feb 02, add option to use mean of Badnell DR in place
-                                         * of these hacks */
-                                        if( ionbal.lg_use_DR_Badnell_rate_coef_mean_ion )
-                                                DielRecomRateCoef_LowT[ion] = ionbal.DR_Badnell_rate_coef_mean_ion[ion];
-                                }
+                                /* use mean of Badnell DR rates */
+                                DielRecomRateCoef_LowT[ion] = ionbal.DR_Badnell_rate_coef_mean_ion[ion];
                                 /* include optional noise here 
                                  * >>chng 06 feb 07, move noise down to here so that use for both

branches/newmole/source/ion_recomb_Badnell.cpp

@@ -273 +273 @@
         char *chs;
 
-#define BIGGEST_INDEX_TO_USE    103
+        const int BIGGEST_INDEX_TO_USE = 103;
 
         /* Declaration of data file name array - done by Kausalya */
@@ -974 +974 @@
 
         /* save mean rec coefficients - used to derive rates for those ions with none */
-        for( ion=0; ion<LIMELM; ++ion )
-        {
-                ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 0.;
+        for( ion=0; ion < LIMELM; ++ion )
                 nsumrec[ion] = 0;
-        }
+
         TeUsed = phycon.te;
         /* NB - for following loop important to go over all elements and all
         * ions, not just active ones, since mean DR is used as the guess for
-        * DR rates for unknown species.  */
-        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
-        {
-                for( ion=0; ion<nelem+1; ++ion )
+        * DR rates for unknown species. */
+        multi_arr<double,2> DR_Badnell_rate_stack( LIMELM, LIMELM );
+        for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
+        {
+                for( ion=0; ion < nelem+1; ++ion )
                 {
                         long int n_bnd_elec_before_recom ,
@@ -992 +991 @@
                         n_bnd_elec_before_recom = nelem-ion;
                         n_bnd_elec_after_recom = nelem-ion+1;
-#                       define DR2SMALL 1e-15
                         /* Badnell dielectronic recombination rate coefficients */
                         if( (ionbal.DR_Badnell_rate_coef[nelem][ion] = 
@@ -1003 +1001 @@
                         {
                                 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
-                                if( ionbal.DR_Badnell_rate_coef[nelem][ion] > DR2SMALL)
-                                {
-                                        /* keep track of mean DR rate for this ion */
-                                        ++nsumrec[ion];
-                                        ionbal.DR_Badnell_rate_coef_mean_ion[ion] +=
-                                                log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
-                                }
+                                /* keep track of mean DR rate for this ion */
+                                DR_Badnell_rate_stack[ion][nsumrec[ion]++] =
+                                        ionbal.DR_Badnell_rate_coef[nelem][ion];
                         }
                         else
@@ -1113 +1107 @@
         ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
         ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
-        if( ionbal.DR_Badnell_rate_coef[nelem][ion]>DR2SMALL )
-        {
-                ++nsumrec[ion];
-                ionbal.DR_Badnell_rate_coef_mean_ion[ion] += 
-                        log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
-        }
+        DR_Badnell_rate_stack[ion][nsumrec[ion]++] =
+                ionbal.DR_Badnell_rate_coef[nelem][ion];
 
         /* >>chng 06 jun 21 by Mitchell Martin, added new DR rate coefficient
@@ -1138 +1128 @@
                 ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] = true;
                 ionbal.DR_Badnell_rate_coef[nelem][ion] = fitSum / phycon.te32;
-                if( ionbal.DR_Badnell_rate_coef[nelem][ion] > DR2SMALL )
-                {
-                        ++nsumrec[ion];
-                        ionbal.DR_Badnell_rate_coef_mean_ion[ion] += 
-                                log10(ionbal.DR_Badnell_rate_coef[nelem][ion]);
-                }
+                DR_Badnell_rate_stack[ion][nsumrec[ion]++] =
+                        ionbal.DR_Badnell_rate_coef[nelem][ion];
         }
 
         /* this is the temperature in eV evaluated to the 3/2 power */
-        double te_eV32 = pow( phycon.te_eV, 1.5);
+        double te_eV32 = pow( phycon.te_eV, 1.5 );
 
         /* >>chng 06 jul 07 by Mitchell Martin, added DR rate coefficient 
@@ -1165 +1151 @@
                         ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion+5] = true;
                         ionbal.DR_Badnell_rate_coef[nelem][ion+5] = fitSum / te_eV32;
-                        if( ionbal.DR_Badnell_rate_coef[nelem][ion+5] > DR2SMALL )
-                        {
-                                ++nsumrec[ion+5];
-                                ionbal.DR_Badnell_rate_coef_mean_ion[ion+5] +=
-                                        log10(ionbal.DR_Badnell_rate_coef[nelem][ion+5]);
-                        }
+                        DR_Badnell_rate_stack[ion+5][nsumrec[ion+5]++] =
+                                ionbal.DR_Badnell_rate_coef[nelem][ion+5];
                 }
         }
@@ -1176 +1158 @@
 
         /* now get mean of the DR rates - may be used for ions with no DR rates */
-        for( ion=0; ion<LIMELM; ++ion )
-        {
+        for( ion=0; ion < LIMELM; ++ion )
+        {
+                ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 0.;
                 if( nsumrec[ion] > 0 )
-                        ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 
-                        pow(10., ionbal.DR_Badnell_rate_coef_mean_ion[ion]/nsumrec[ion]);
+                {
+                        // we have collected all the DR rates above, now we calculate the logarithmic average
+                        // of the top 3/4 of the entries. this procedure has been checked not to result in
+                        // discontinuities in the average rate as a function of temperature.
+                        // NB NB if this code is changed, check again that the resulting average rates
+                        // are continuous, otherwise discontiuities in the ionization balance and heating-
+                        // cooling function may result.
+
+                        // first sort the stack so that we can easily determine the largest rates
+                        sort( DR_Badnell_rate_stack.ptr( ion, 0 ), DR_Badnell_rate_stack.ptr( ion, nsumrec[ion] ) );
+
+                        double maxval = DR_Badnell_rate_stack[ion][nsumrec[ion]-1];
+                        if( maxval > 0. )
+                        {
+                                // now calculate logarithmic average of top three quarters of entries
+                                // we discard the lower quarter because many of these are very low or 0
+                                // NB NB we MUST use a fixed number of values at any temperature to
+                                // calculate the average since discontinuities would result otherwise!
+                                double sum = 0.;
+                                for( int j = nsumrec[ion]/4; j < nsumrec[ion]; ++j )
+                                {
+                                        double rate_one = log(max(DR_Badnell_rate_stack[ion][j],1e-5*maxval));
+                                        sum += rate_one;
+                                }
+                                sum /= (double)(nsumrec[ion] - nsumrec[ion]/4);
+                                ionbal.DR_Badnell_rate_coef_mean_ion[ion] = exp(sum);
+                        }
+                }
                 /*fprintf(ioQQQ,"DEBUG %li %.2e\n", ion , ionbal.DR_Badnell_rate_coef_mean_ion[ion] );*/
         }
@@ -1190 +1199 @@
         for( ion=0; ion<5; ion++ )
         {
-
                 /* only do this rate if not already done by a previous approximation
                 * so following used for ion <= 3 */
@@ -1289 +1297 @@
                 cdEXIT( EXIT_SUCCESS );
         }
-
         return;
 }