Changeset 2276

Timestamp:

07/28/08 01:45:32 (6 months ago)

Author:

gary

Message:

Lya maser drives negative solution to 21 cm levels - add ticket on this - current solution is to consider only non-negative Lya photon occupation numbers

Location:

trunk/source

Files:

• atom_hyperfine.cpp (modified) (4 diffs)
• lines_service.cpp (modified) (3 diffs)

trunk/source/atom_hyperfine.cpp

@@ -10 +10 @@
 /*H21cm_proton - evaluate proton spin changing H atom collision rate, */
 #include "cddefines.h"
+#include "conv.h"
 #include "lines_service.h"
 #include "phycon.h"
@@ -60 +61 @@
         ASSERT( a21>0. );
 
-        /* >>chng 04 dec 01, add hyperfine.lgLya_pump_21cm, option to turn off Lya pump
-         * of 21 cm, with no 21cm lya pump */
+        /* hyperfine.lgLya_pump_21cm is option to turn off Lya pump
+         * of 21 cm, with no 21cm lya pump command - note that this
+         * can be negative if Lya mases - can occur during search phase */
         occnu_lya = OccupationNumberLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ) *
                 hyperfine.lgLya_pump_21cm;
-
-        /* >>chng 05 apr 20, GS, the lya occupation number for the hyperfine levels 0S1/2 and 1S1/2 are different*/
+        if( occnu_lya<0 )
+        {
+                static bool lgCommentDone = false;
+                /* Lya is masing - could be due to very bad solution in search phase */
+                if( !conv.lgSearch && !lgCommentDone )
+                {
+                        fprintf(ioQQQ,
+                        "NOTE Lya masing will invert 21 cm, occupation number set zero\n");
+                        lgCommentDone = true;
+                }
+                occnu_lya = 0.;
+        }
+
+        /* Lya occupation number for the hyperfine levels 0S1/2 and 1S1/2 are different*/
         texc = TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );
-        /* >>chng 05 apr 21, GS, Energy difference between 2p1/2 and 2p3/2 is taken from NSRDS */
+        /* Energy difference between 2p1/2 and 2p3/2 taken from NSRDS */
         if( texc > 0. )
         {
-                /* convert to boltz factor, which will applied to occupation number of hiher energy transition */
+                /* convert to Boltzmann factor, which will applied to occupation 
+                 * number of higher energy transition */
                 texc1 = sexp(0.068/texc);
                 texc2 = sexp(((82259.272-82258.907)*T1CM)/texc);
@@ -134 +149 @@
         HFLines[0].Hi->Pop = (x/(1.+x))* PopTot;
         HFLines[0].Lo->Pop = (1./(1.+x))* PopTot;
+        ASSERT( HFLines[0].Hi->Pop >0. );
+        ASSERT( HFLines[0].Lo->Pop >0. );
 
         /* the population with correction for stimulated emission */
@@ -140 +157 @@
         /* number of escaping line photons, used elsewhere for outward beam */
         HFLines[0].Emis->phots = a21*HFLines[0].Hi->Pop;
-
+        ASSERT( HFLines[0].Emis->phots >= 0. );
         /* intensity of line */
         HFLines[0].Emis->xIntensity = HFLines[0].Emis->phots*HFLines[0].EnergyErg;

trunk/source/lines_service.cpp

@@ -142 +142 @@
         ASSERT( t->ipCont > 0 );
 
-        /* routine to evaluate line photon occupation number */
-        if( t->Lo->Pop > SMALLFLOAT )
+        /* routine to evaluate line photon occupation number - 
+         * return negative number if line is maser */
+        if( fabs(t->Emis->PopOpc) > SMALLFLOAT )
         {
                 /* the lower population with correction for stimulated emission */
-                double PopLo_corr = t->Lo->Pop / t->Lo->g - t->Hi->Pop / t->Hi->g;
-                OccupationNumberLine_v = ( t->Hi->Pop / t->Hi->g )/SDIV(PopLo_corr ) *
+                OccupationNumberLine_v = ( t->Hi->Pop / t->Hi->g ) /
+                        ( t->Emis->PopOpc / t->Lo->g )  *
                         (1. - t->Emis->Pesc);
         }
@@ -154 +155 @@
                 OccupationNumberLine_v = 0.;
         }
+        /* return value is not guaranteed to be positive - negative if
+         * line mases */
         return( OccupationNumberLine_v );
 }
@@ -407 +410 @@
   transition * t)
 {
-
         DEBUG_ENTRY( "PutCS()" );
 
-        /* collision strength must not be negative, had been test for being positive,
-         * but called with zero - did not check why 98 jul 5 */
-        ASSERT( cs >= 0. );
+        /* collision strength must be negative */
+        ASSERT( cs > 0. );
 
         t->Coll.col_str = (realnum)cs;
+
         return;
 }