Changeset 2125 for branches/newmole

Timestamp:

06/24/08 16:39:58 (7 months ago)

Author:

rjrw

Message:

Merged from trunk r2111:2124

Location:

branches/newmole/source

Files:

• cloudy.cpp (modified) (1 diff)
• cont_createmesh.cpp (modified) (1 diff)
• conv_base.cpp (modified) (2 diffs)
• conv_init_solution.cpp (modified) (1 diff)
• dynamics.cpp (modified) (1 diff)
• iso_level.cpp (modified) (2 diffs)
• iter_end_chk.cpp (modified) (1 diff)
• mole_h2_io.cpp (modified) (2 diffs)
• parse_punch.cpp (modified) (2 diffs)
• punch_do.cpp (modified) (2 diffs)
• rfield.h (modified) (3 diffs)
• rt_diffuse.cpp (modified) (1 diff)

branches/newmole/source/cloudy.cpp

@@ -223 +223 @@
 
                         /* converge the pressure-temperature-ionization solution for this zone 
-                         * NB ignoring return value - should be ok (return 1 for disaster) */
+                         * NB ignoring return value - should be ok (return 1 for abort) 
+                         * we can't break here in case of abort since there is still housekeeping
+                         * that must be done in following routines */
                         ConvPresTempEdenIoniz();
 

branches/newmole/source/cont_createmesh.cpp

@@ -147 +147 @@
         /* nend_max+1 because search phase is zone 0, first zone at illumin face is 1 */
         rfield.ConEmitLocal = (realnum**)MALLOC( (size_t)(geometry.nend_max+1)*sizeof(realnum *)  );
+        rfield.ConSourceFcnLocal = (realnum**)MALLOC( (size_t)(geometry.nend_max+1)*sizeof(realnum *)  );
         for( i=0;i<(geometry.nend_max+1); ++i )
         {
                 rfield.ConEmitLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum)  );
+                rfield.ConSourceFcnLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum)  );
+        }
+        for( i=0;i<(geometry.nend_max+1); ++i )
+        {
+                for( long j=0; j<rfield.nupper; ++j)
+                {
+                        rfield.ConSourceFcnLocal[i][j] = 1.;
+                }
         }
 

branches/newmole/source/conv_base.cpp

@@ -855 +855 @@
         PresTotCurrent();
 
-        if( trace.lgTrace )
-        {
-                fprintf( ioQQQ, 
-                        "   ConvBase return. %.2f Te:%.3e  HI:%.3e  HII:%.3e  H2:%.3e  Ne:%.3e  htot:%.3e  CSUP:%.2e Conv?%c\n", 
-                  fnzone,
-                  phycon.te, 
-                  dense.xIonDense[ipHYDROGEN][0], 
-                  dense.xIonDense[ipHYDROGEN][1], 
-                  hmi.H2_total,
-                  dense.eden, 
-                  thermal.htot, 
-                  secondaries.csupra[ipHYDROGEN][0] ,
-                  TorF(conv.lgConvIoniz) );
-        }
-
         ASSERT(lgElemsConserved());
 
@@ -879 +864 @@
          * reset to zero each time ConvPresTempEdenIoniz is called */
         ++conv.nPres2Ioniz;
+
+        if( trace.lgTrace )
+        {
+                fprintf( ioQQQ, 
+                        "   ConvBase return. fnzone %.2f nPres2Ioniz %li Te:%.3e  HI:%.3e  HII:%.3e  H2:%.3e  Ne:%.3e  htot:%.3e  CSUP:%.2e Conv?%c reason:%s\n", 
+                        fnzone,
+                        conv.nPres2Ioniz,
+                        phycon.te, 
+                        dense.xIonDense[ipHYDROGEN][0], 
+                        dense.xIonDense[ipHYDROGEN][1], 
+                        hmi.H2_total,
+                        dense.eden, 
+                        thermal.htot, 
+                        secondaries.csupra[ipHYDROGEN][0] ,
+                        TorF(conv.lgConvIoniz) ,
+                        conv.chConvIoniz);
+        }
 
         /* this is abort option set with SET PRESIONIZ command,

branches/newmole/source/conv_init_solution.cpp

@@ -319 +319 @@
                 PresTotCurrent();
 
-                /* now get final temperature */
+                /* now get new initial temperature */
                 if( ConvPresTempEdenIoniz() )
                 {

branches/newmole/source/dynamics.cpp

@@ -1850 +1850 @@
         if( lgWind )
         {
-                /* increse precision of solution */
+                /* increase precision of solution */
                 conv.EdenErrorAllowed = 1e-3;
                 /* the actual relative error is relative to the total heating and cooling,

branches/newmole/source/iso_level.cpp

@@ -202 +202 @@
         else
         {
+                /* this branch is main solver, full level populations 
+                 * assert since this code must change if NISO ever increased */
                 ASSERT( NISO == 2 );
                 long SpinUsed[NISO] = {2,1};
@@ -455 +457 @@
                  *
                  * ==================================================================== */
-
                 /* save matrix, this allocates SaveZ */
                 SaveZ = z;

branches/newmole/source/iter_end_chk.cpp

@@ -161 +161 @@
         }
 
+        else if( (StopCalc.nTotalIonizStop>0) && conv.nTotalIoniz> StopCalc.nTotalIonizStop )
+        {
+                /* stop if exceed number of calls to conv base set with 
+                 * stop nTotalIonizStop command */
+                lgDone = true;
+                strncpy( StopCalc.chReasonStop, "nTotalIonizStop reached.", sizeof(StopCalc.chReasonStop) );
+        }
+
         /* this flag says that 21cm line optical depth is the stop quantity */
         else if( StopCalc.lgStop21cm && (HFLines[0].Emis->TauCon >=  StopCalc.tauend) )

branches/newmole/source/mole_h2_io.cpp

@@ -272 +272 @@
                 strcpy( punch.chPunch[punch.npunch], "H2ln" );
                 sprintf( chHeader, 
-                        "#H2 line\tEhi\tVhi\tJhi\tElo\tVlo\tJlo\twl(mic)\twl(lab)\tlog L or I\tI/Inorm\tExcit(hi, K)\tg_u h\\nu * Aul\n" );
+                        "#H2 line\tEhi\tVhi\tJhi\tElo\tVlo\tJlo\twl(mic)\twl(lab)\tlog L or I\tI/Inorm\tExcit(hi, K)\tg_u h nu * Aul\n" );
                 /* first optional number changes the threshold of weakest line to print*/
                 /* fe2thresh is intensity relative to normalization line,
@@ -2192 +2192 @@
                                                                                 /* excitation energy of upper level in K */
                                                                                 fprintf( io, "\t%.3f", energy_wn[iElecHi][iVibHi][iRotHi]*T1CM );
-                                                                                /* the product h nu * Aul */
+                                                                                /* the product g_hi h nu * Aul */
                                                                                 ASSERT( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 );
                                                                                 fprintf( io, "\t%.3e", H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul* 

branches/newmole/source/parse_punch.cpp

@@ -2007 +2007 @@
         else if( nMatch("SOUR",chCard) )
         {
+
+                /* check for keyword UNITS on line, then scan wavelength or energy units if present,
+                * units are copied into punch.chConPunEnr */
+                ChkUnits(chCard);
+
                 if( nMatch("DEPT",chCard) )
                 {
@@ -2019 +2024 @@
                         strcpy( punch.chPunch[punch.npunch], "SOUS" );
                         sprintf( chHeader, 
-                                "#continuum source function\n" );
+                                "#continuum source function nu(Ryd)\tConEmitLocal/widflx"
+                                "\tabs opac\tConSourceFcnLocal\tConSourceFcnLocal[nzone][j]/plankf\n" );
                 }
                 else

branches/newmole/source/punch_do.cpp

@@ -453 +453 @@
 
                                                 /* the outward emitted continuum */
-                                                conem = ((rfield.ConEmitOut[j])/
+                                                conem = (rfield.ConEmitOut[j]/
                                                   rfield.widflx[j]*rfield.anu2[j] + punch.PunchLWidth*
                                                   rfield.outlin[j]*rfield.anu[j])*radius.r1r0sq*
@@ -2901 +2901 @@
                                         {
                                                 fprintf( punch.ipPnunit[ipPun], 
-                                                        "%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 
-                                                  rfield.anu[j], 
+                                                        "%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 
+                                                  AnuUnit(rfield.AnuOrg[j]),
                                                   rfield.ConEmitLocal[nzone][j]/rfield.widflx[j], 
                                                   opac.opacity_abs[j], 
-                                                  rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j]), 
-                                                  rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]/SDIV( opac.opacity_abs[j])/plankf(j) );
+                                                  rfield.ConSourceFcnLocal[nzone][j], 
+                                                  rfield.ConSourceFcnLocal[nzone][j]/plankf(j) );
                                         }
                                 }

branches/newmole/source/rfield.h

@@ -124 +124 @@
         realnum **ConEmitLocal/* [depth][energy]*/;
 
+        /** the local source function - diffuse emission, photons Ryd-1 s-1 */
+        realnum **ConSourceFcnLocal/* [depth][energy]*/;
+
         /** reflected diffuse emission continuum */
         realnum *ConEmitReflec;
@@ -129 +132 @@
         /** outward diffuse emission continuum (not the interactive one),
          * this is incremented in radinc because of interplay between absorption
-         * and emission - get the outward bremsstrahlung right */
+         * and emission - get the outward bremsstrahlung right 
+         * photons cell-1 cm-2 s-1 */
         realnum *ConEmitOut;
 
@@ -168 +172 @@
                 **otssav/*[NC_ELL][2]*/;
 
-        /** outward directed line emission */
+        /** outward directed line emission photons cm-2 s-1 */
         realnum *outlin,
                 *outlin_noplot;

branches/newmole/source/rt_diffuse.cpp

@@ -738 +738 @@
                 /* >>chng 00 oct 03, add diffuse continua */
                 /* local diffuse continua */
-                rfield.OccNumbDiffCont[nu] = rfield.ConEmitLocal[nzone][nu]*rfield.convoc[nu];
+                rfield.OccNumbDiffCont[nu] = 
+                        rfield.ConEmitLocal[nzone][nu]*rfield.convoc[nu];
+
+                rfield.ConSourceFcnLocal[nzone][nu] = 
+                        rfield.ConEmitLocal[nzone][nu]/rfield.widflx[nu]/
+                        (realnum)SDIV( opac.opacity_abs[nu]);
 
                 /* confirm that all are non-negative */