Changeset 2440 for branches/newmole

Timestamp:

11/04/08 14:25:16 (2 months ago)

Author:

rjrw

Message:

Merged from trunk r2405:2439.

Location:

branches/newmole

Files:

• docs/latex/Hazy1 (copied) (copied from trunk/docs/latex/Hazy1)
• source/atmdat_readin.cpp (modified) (3 diffs)
• source/atoms.h (modified) (4 diffs)
• source/atom_feii.cpp (modified) (1 diff)
• source/atom_oi.cpp (modified) (8 diffs)
• source/cddefines.h (modified) (1 diff)
• source/conv_base.cpp (modified) (1 diff)
• source/date.h (modified) (1 diff)
• source/dynamics.cpp (modified) (1 diff)
• source/dynamics.h (modified) (1 diff)
• source/ion_solver.cpp (modified) (1 diff)
• source/iso_level.cpp (modified) (1 diff)
• source/lines_service.cpp (modified) (2 diffs)
• source/lines_service.h (modified) (1 diff)
• source/mole_h2.cpp (modified) (1 diff)
• source/prt_final.cpp (modified) (1 diff)
• source/rt.h (modified) (1 diff)
• source/rt_line_all.cpp (modified) (13 diffs)
• source/rt_line_one.cpp (modified) (15 diffs)
• tsuite/auto/blr_f92.in (modified) (2 diffs)
• tsuite/auto/blr_n12_p19.in (modified) (1 diff)
• tsuite/auto/blr_n12_p19_Z20.in (modified) (1 diff)
• tsuite/auto/blr_n13_p18.in (modified) (1 diff)
• tsuite/auto/dynamics_orion_flow.in (modified) (1 diff)
• tsuite/auto/nova_photos.in (modified) (6 diffs)
• tsuite/slow/feii_blr_n11_p20_Z20.in (modified) (1 diff)
• tsuite/slow/feii_blr_n12_p19.in (modified) (1 diff)
• tsuite/slow/feii_blr_n12_p19_Z20.in (modified) (1 diff)
• tsuite/slow/feii_blr_n13_p18_Z20.in (modified) (1 diff)
• tsuite/slow/feii_blr_n13_p22.in (modified) (2 diffs)

branches/newmole/source/atmdat_readin.cpp

@@ -1050 +1050 @@
         }
 
-        nUTA = UTALines.size();
+        nUTA = (long)UTALines.size();
 
         /* option to dump UTA lines */
@@ -1640 +1640 @@
                         UTA.back().Hi->IonStg = ion+1;
 
-                        UTA.back().Hi->g = level[ind_hi].g;
-                        UTA.back().Lo->g = level[ind_lo].g;
+                        UTA.back().Hi->g = (realnum)level[ind_hi].g;
+                        UTA.back().Lo->g = (realnum)level[ind_lo].g;
 
                         double WavNum = edif*RYD_INF;
@@ -1659 +1659 @@
                         /* Badnell gives UPWARD transition rate Alu, an unusual notation,
                          * convert it here to the normal downward transition rate Aul */
-                        UTA.back().Emis->Aul = Bij*UTA.back().Lo->g/UTA.back().Hi->g;
+                        UTA.back().Emis->Aul = (realnum)(Bij*UTA.back().Lo->g/UTA.back().Hi->g);
                         UTA.back().Emis->gf =
                                 (realnum)GetGF( UTA.back().Emis->Aul, UTA.back().EnergyWN, UTA.back().Hi->g );

branches/newmole/source/atoms.h

@@ -4 +4 @@
 #ifndef _ATOMS_H_
 #define _ATOMS_H_
-
 
  /**
@@ -219 +218 @@
 
 /** number of levels in OI atom */
-#define N_OI_LEVELS     6
+const int N_OI_LEVELS = 6;
+const long LIMLEVELN = 20L;
 
 EXTERN struct t_atoms {
@@ -240 +240 @@
         /** populations from OI fluorescense problem */
         realnum popoi[N_OI_LEVELS];
-
-        double pmph31,
-                esch31;
 
         realnum pmpo51, 
@@ -252 +249 @@
                 popmg2;
 
-#define LIMLEVELN       20L
         /**
          * this stores most recently evaluated level populations

branches/newmole/source/atom_feii.cpp

@@ -1584 +1584 @@
                                 /*RT_line_one do rt for emission line structure - 
                                  * calls RT_line_static or RT_line_wind */
-                                RT_line_one( &Fe2LevN[ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,true);
+                                RT_line_one( &Fe2LevN[ipHi][ipLo], lgDoEsc, lgUpdateFineOpac, true, 0.f, 1.f );
                         }
                 }

branches/newmole/source/atom_oi.cpp

@@ -18 +18 @@
 /*oi_level_pops get OI level population with Ly-beta pumping */
 STATIC void oi_level_pops(double abundoi, 
-  double *coloi);
+                          double *coloi);
 
 /*atom_oi drive the solution of OI level populations, Ly-beta pumping */
 void atom_oi_calc(double *coloi)
 {
-        long int i;
-        double esin, 
+        double esab,
           eslb, 
           esoi, 
@@ -31 +30 @@
           opaclb, 
           opacoi, 
-          tin, 
-          tout, 
           xlb, 
           xoi;
-        static double esab;
-        static double aoi = TauLines[ipTO1025].Emis->Aul;
-        static double alb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul;
+        double aoi = TauLines[ipTO1025].Emis->Aul;
+        double alb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul;
 
         DEBUG_ENTRY( "atom_oi_calc()" );
 
+        fixit();
+        // The code below should be calculating the O I 1025 pumping by H Ly beta, as well as
+        // the inverse process (this can become important in hydrogen-deficient environments).
+        // It now uses the Elitzur & Netzer (1985, ApJ, 291, 464) theory, which is no longer
+        // valid since the line overlap code prevents us from getting at the escape probability
+        // of individual lines.
+
         /* A's from Pradhan; OI pump line; Ly beta, 8446 */
-
-        /* Notation;
-         * PMPH31 net rate hydrogen level 3 depopulated to 1
-         * PMPO15 and PMPO51 are similar, but for oxygen
-         * ESCH31 is emission in 1025 transition */
 
         /* called by LINES before calc really start, protect here
@@ -52 +50 @@
         if( dense.xIonDense[ipOXYGEN][0] <= 0. )
         {
-                for( i=0; i < 6; i++ )
+                for( int i=0; i < 6; i++ )
                 {
                         atoms.popoi[i] = 0.;
                 }
-                /* return zero */
                 *coloi = 0.;
                 atoms.pmpo15 = 0.;
                 atoms.pmpo51 = 0.;
-                atoms.pmph31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul*
-                        (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc +
-                  Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc);
-
-                atoms.esch31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul*
-                        (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc +
-                        Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc);
-
-                /* all trace output turned on with "trace ly beta command' */
-                if( trace.lgTr8446 && trace.lgTrace )
-                {
-                        fprintf( ioQQQ, 
-                                "       P8446 called for first time, finds A*escape prob 3-1 =%10.2e\n", 
-                          atoms.esch31 );
-                }
                 return;
         }
 
-        if( iteration > 1 )
-        {
-                /* two sided escape prob */
-                tin = TauLines[ipTO1025].Emis->TauIn;
-                esin = esc_CRDwing_1side(tin,1e-4);
-                tout = (TauLines[ipTO1025].Emis->TauTot)*0.9 - 
-                  TauLines[ipTO1025].Emis->TauIn;
-
-                if( trace.lgTr8446 && trace.lgTrace )
-                {
-                        fprintf( ioQQQ, "       P8446 tin, tout=%10.2e%10.2e\n", tin, tout );
-                }
-
-                if( tout > 0. )
-                {
-                        tout = TauLines[ipTO1025].Emis->TauTot - TauLines[ipTO1025].Emis->TauIn;
-
-                        /* do not update esab if we overran optical depth scale */
-                        esab = 0.5*(esin + esc_CRDwing_1side(tout,1e-4));
-                }
-        }
-        else
-        {
-                /* one-sided escape probability */
-                esab = esc_CRDwing_1side(TauLines[ipTO1025].Emis->TauIn,1e-4);
-        }
-
-        esoi =TauLines[ipTO1025].Emis->Pelec_esc + TauLines[ipTO1025].Emis->Pesc;
-        eslb =Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + 
+        // line overlap code makes this the escape probability of the combined lines
+        esab = TauLines[ipTO1025].Emis->Pelec_esc + TauLines[ipTO1025].Emis->Pesc;
+
+        // these two are no longer correct, the line overlap code makes it impossible
+        // to get at the escape probabilities of the individual lines...
+        esoi = TauLines[ipTO1025].Emis->Pelec_esc + TauLines[ipTO1025].Emis->Pesc;
+        eslb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + 
                 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc;
 
@@ -118 +78 @@
 
         /* find relative opacities for two lines */
-        opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/DoppVel.doppler[7];
-        opaclb = 1.22e-8*dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/DoppVel.doppler[0];
+        opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/DoppVel.doppler[ipOXYGEN];
+        opaclb = 1.22e-8*dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/DoppVel.doppler[ipHYDROGEN];
 
         /* these are x sub a (OI) and x sub b (ly beta) defined in Elitz+Netz */
@@ -128 +88 @@
         foi = MIN2(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[ipOXYGEN];
         flb = MIN2(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[ipHYDROGEN]*
-          MAX2(0.,1.- TauLines[ipTO1025].Emis->Pelec_esc - TauLines[ipTO1025].Emis->Pesc);
+                MAX2(0.,1.- TauLines[ipTO1025].Emis->Pelec_esc - TauLines[ipTO1025].Emis->Pesc);
 
         if( trace.lgTr8446 && trace.lgTrace )
@@ -139 +99 @@
         /* pumping of OI by L-beta - this goes into OI matrix as 1-5 rate
          * lgInducProcess set false with no induced command, usually true */
+        /* Notation: pmpo15 net rate oxygen level 5 populated from 1 */
         if( rfield.lgInducProcess )
         {
-                atoms.pmpo15 = (realnum)(flb*alb*(StatesElem[ipH_LIKE][ipHYDROGEN][ipH3p].Pop*dense.xIonDense[ipHYDROGEN][1])*xoi*(1. - esab)/
-                  dense.xIonDense[ipOXYGEN][0]);
+                atoms.pmpo15 = (realnum)(flb*alb*dense.xIonDense[ipHYDROGEN][1]*
+                                         StatesElem[ipH_LIKE][ipHYDROGEN][ipH3p].Pop*
+                                         xoi*(1. - esab)/dense.xIonDense[ipOXYGEN][0]);
                 /* net decay rate from upper level */
-                atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - 
-                  esab)*foi));
+                atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - esab)*foi));
         }
         else
@@ -155 +116 @@
         /* find level populations for OI */
         oi_level_pops(dense.xIonDense[ipOXYGEN][0],coloi);
-
-        /* escape term from n=3 of H; followed by pump term to n=3 */
-        /** \todo 2 this is not used, should it be? */
-        //atoms.pmph31 = alb*(1. - (1. - flb)*(1. - eslb) - xlb*(1. - esab)*flb);
-
-        /*pmph13 = xlb*(1. - esab)*foi*aoi*atoms.popoi[4];*/
-        atoms.pmph31 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul*
-                (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc);
-
-        /* following is actual emission rate, used to predict Ly-beta in lines.for */
-        atoms.esch31 = alb*(eslb*(1. - flb) + esab*flb);
-
-        /* all trace output turned on with "trace ly beta command' */
-        if( trace.lgTr8446 && trace.lgTrace )
-        {
-                fprintf( ioQQQ, "       P8446 PMPH31=%10.2e EscP3-1%10.2e ESCH31=%10.2e\n", 
-                  atoms.pmph31, atoms.pmph31/alb, atoms.esch31/alb );
-        }
 
         /* continuum pumping due to J=1, 0 sub states.
@@ -220 +163 @@
 /*oi_level_pops get OI level population with Ly-beta pumping */
 STATIC void oi_level_pops(double abundoi, 
-  double *coloi)
+                          double *coloi)
 {
         bool lgNegPop;

branches/newmole/source/cddefines.h

@@ -1559 +1559 @@
 }
 
+/**get_ptr attribute shim to get raw pointer to contained data with
+ * correct type */
+template<class T> inline T* get_ptr(valarray<T> &v)
+{
+        return &v[0];
+}
+template<class T> inline T* get_ptr(vector<T> &v)
+{
+  return &v[0];
+}
+template<class T> inline const T* get_ptr(const valarray<T> &v)
+{
+        return const_cast<const T*>(&const_cast<valarray<T>&>(v)[0]);
+}
+template<class T> inline const T* get_ptr(const vector<T> &v)
+{       
+  return const_cast<const T*>(&const_cast<vector<T>&>(v)[0]);
+}
+
+
 /**************************************************************************
  *

branches/newmole/source/conv_base.cpp

@@ -629 +629 @@
                                                          * activate.  when it passes rm this block and 
                                                          * change err_tol in following block to a small number */
-                                        if(0 &  conv.lgConvIoniz && 
+                                        if( 0 && conv.lgConvIoniz && 
                                                 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 &&
                                                 dense.xIonDense[nelem][nelem+1-ipISO]>SMALLFLOAT &&

branches/newmole/source/date.h

@@ -16 +16 @@
 #define MONTH   9
 /* day is correct */
-#define DAY     9
+#define DAY     26
 
 #endif /* _DATE_H_ */

branches/newmole/source/dynamics.cpp

@@ -1655 +1655 @@
         dynamics.lgRecom = false;
 
+        /* don't force populations to equilibrium levels */
+        dynamics.lgEquilibrium = false;
+
         /* set true if time dependent calculation is finished */
         dynamics.lgStatic_completed = false;

branches/newmole/source/dynamics.h

@@ -165 +165 @@
         realnum dDensityDT;
 
+        /** Enforce equilibrium populations */
+        bool lgEquilibrium;
 
         /* set true with trace option on time command */

branches/newmole/source/ion_solver.cpp

@@ -376 +376 @@
 
         /* chng 03 jan 13 rjrw, add in dynamics if required here,
-         * last test - only do advection if we have not overrun the radius scale */
-        if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax+1 
+         * - only do advection if we have not overrun the radius scale */
+        if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection && !dynamics.lgEquilibrium
                 /*&& radius.depth < dynamics.oldFullDepth*/ )
         {

branches/newmole/source/iso_level.cpp

@@ -328 +328 @@
                 /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
                 /*>>>chng 02 oct 01, only include if lgAdvection is set */
-                if( dynamics.lgAdvection && dynamics.lgISO[ipISO])
+                if( dynamics.lgAdvection && !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
                 {
                         /* add in advection - these terms normally zero */

branches/newmole/source/lines_service.cpp

@@ -23 +23 @@
 /*EmLineZero set all elements of transition struc to zero */
 /*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
-/*lgTauGood returns true is we have not overrun optical depth scale */
 /*OccupationNumberLine - derive the photon occupation number at line center for any line */
 /*outline - adds line photons to reflin and outlin */
@@ -1154 +1153 @@
 }
 
-/* returns true is we have not overrun optical depth scale */
-bool lgTauGood( transition * t)
-{
-        bool lgGoodTau;
-
-        DEBUG_ENTRY( "lgTauGood()" );
-
-        if( (t->Emis->TauTot*0.9 - t->Emis->TauIn < 0. && t->Emis->TauIn > 0.) && 
-            ! fp_equal( t->Emis->TauTot, opac.taumin ) )
-        {
-                /* do not do anything if we have overrun the scale, */
-                lgGoodTau = false;
-        }
-        else
-        {
-                lgGoodTau = true;
-        }
-        return lgGoodTau;
-}
-
 /*gbar0 compute g-bar gaunt factor for neutrals */
 STATIC void gbar0(double ex, 

branches/newmole/source/lines_service.h

@@ -220 +220 @@
 \param *t
 */
-bool lgTauGood( transition * t);
+inline bool lgTauGood( transition* t )
+{
+        return ( iteration == 1 || t->Emis->TauIn <= 0. || t->Emis->TauTot*0.9 - t->Emis->TauIn >= 0. );
+}
 
 /**MakeCS compute collision strength by g-bar approximations 

branches/newmole/source/mole_h2.cpp

@@ -663 +663 @@
                                                                  * variation in the line pumping rate, which made H2 abundance fluctuate due to
                                                                  * Solomon process having slight dr-caused mole. */
-                                                                RT_line_one( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo], lgDoEsc, lgUpdateFineOpac, false );
+                                                                RT_line_one( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo], lgDoEsc, lgUpdateFineOpac, false, 0.f, 1.f );
                                                         }
                                                 }

branches/newmole/source/prt_final.cpp

@@ -802 +802 @@
                 fprintf( ioQQQ, "%s" , chIntensityType[ipEmType] );
                 fprintf( ioQQQ, " line intensities\n" );
+                // caution about emergent intensities when outward optical
+                // depths are not yet known
+                if( ipEmType==1 && iteration==1 )
+                        fprintf(ioQQQ," Caution: emergent intensities are not reliable on the "
+                        "first iteration.\n");
 
                 /* option to only print brighter lines */

branches/newmole/source/rt.h

@@ -32 +32 @@
                 this can cause pump to depend on zone thickness, and leads to unstable
                 feedback in some models with the large H2 molecule, due to Solomon
-                process depending on zone thickness and level populations. 
-*/
-void RT_line_one(transition * t , 
-        bool lgDoEsc , 
-        bool lgUpdateFineOpac,
-        bool lgShield_this_zone );
+                process depending on zone thickness and level populations.
+\param pestrk Stark escape probability
+\param efac Exponential correction for single-flight escape probability? usually 1
+*/
+void RT_line_one(transition* t, 
+                 bool lgDoEsc, 
+                 bool lgUpdateFineOpac,
+                 bool lgShield_this_zone,
+                 realnum pestrk,
+                 realnum efac);
 
 /**rt_continuum_shield_fcn computing continuum shielding due to single line 

branches/newmole/source/rt_line_all.cpp

@@ -33 +33 @@
                 nelem;
         long ipHi , ipLo;
-        double factor,
-                coloi;
+        double factor;
         double SaveLyaPesc[NISO][LIMELM] , 
                 SaveLyaPdest[NISO][LIMELM];
@@ -95 +94 @@
         /*if( lgUpdateFineOpac )
                 fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac in rt_line_all\n");*/
-        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
+        for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
         {
                 /* loop over all iso-electronic sequences */
-                for( nelem=ipISO; nelem<LIMELM; ++nelem )
+                for( nelem=ipISO; nelem < LIMELM; ++nelem )
                 {
                         /* parent ion stage, for H is 1, for He is 1 for He-like and 
@@ -108 +107 @@
                                 continue;
                         /* need we consider this ion? */
-                        if( ion <=dense.IonHigh[nelem] )
+                        if( ion <= dense.IonHigh[nelem] )
                         {
                                 /* save escape and destruction probs for each Lya since
@@ -115 +114 @@
                                  * evaluating iso sequence that does not exist for an element,
                                  * as in He-like hydrogen */
-                                if( ipISO<=nelem )
-                                {
-                                        SaveLyaPesc[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc;
-                                        SaveLyaPdest[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest;
-                                }
+                                SaveLyaPesc[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc;
+                                SaveLyaPdest[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest;
 
                                 /* convert pops to per unit vol rather than per ion */
@@ -153 +149 @@
                                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot);*/
 
+                                                fixit(); // is this factor needed? should it be for he-like lya only?
+                                                realnum efac;
+                                                if( ipISO > 0 && ipHi == iso.nLyaLevel[ipISO] && ipLo == 0 )
+                                                        /* do not let too much Lya escape outward since so important 
+                                                         * H Lya is special case done elsewhere */
+                                                        efac = opac.ExpmTau[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
+                                                else
+                                                        efac = 1.f;
+
                                                 /* generate escape prob, pumping rate, destruction prob, 
                                                  * inward outward fracs */
-                                                RT_line_one(&Transitions[ipISO][nelem][ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,true);
-
-
+                                                RT_line_one( &Transitions[ipISO][nelem][ipHi][ipLo],
+                                                             lgDoEsc,
+                                                             lgUpdateFineOpac,
+                                                             true,
+                                                             (realnum)iso.pestrk[ipISO][nelem][ipLo][ipHi],
+                                                             efac );
+
+                                                /* set true to print pump rates*/
+                                                enum {DEBUG_LOC=false};
+              &nbs