Ticket #135: rt.diff

atom_level2.cpp

@@ -12 +12 @@
 #include "rfield.h"
 #include "thermal.h"
 #include "cooling.h"
+#include "conv.h"
+#include "doppvel.h"
+#include "rt.h"
 #include "atoms.h"
 
 void atom_level2(transition * t)
@@ -93 +96 @@
                 return;
         }
 
+        // update RT quantities to keep them synched with level populations
+        if( conv.nTotalIoniz && !conv.lgLastSweepThisZone )
+                RT_line_one( t, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+
         /* net rate down A21*(escape + destruction) */
         a21 = t->Emis->Aul*(t->Emis->Pesc+ t->Emis->Pdest + t->Emis->Pelec_esc);
 

atom_level3.cpp

@@ -4 +4 @@
 #include "cddefines.h"
 #include "phycon.h"
 #include "physconst.h"
+#include "conv.h"
+#include "doppvel.h"
+#include "rt.h"
 #include "dense.h"
 #include "transition.h"
 #include "thermal.h"
@@ -168 +171 @@
                 TotalInsanity();
         }
 
+        long nelem;
         /* abundances from the elements grid
          * one of these must be a true line */
         if( g010 > 0. )
         {
                 /* put null terminated line label into chLab using line array*/
                 chIonLbl(chLab, t10);
+                nelem = t10->Hi->nelem;
                 AbunxIon = dense.xIonDense[ t10->Hi->nelem -1][t10->Hi->IonStg-1];
         }
 
@@ -181 +186 @@
         {
                 /* put null terminated line label into chLab using line array*/
                 chIonLbl(chLab, t21);
+                nelem = t21->Hi->nelem;
                 AbunxIon = dense.xIonDense[t21->Hi->nelem -1][t21->Hi->IonStg-1];
         }
 
@@ -198 +204 @@
         c = t20->EnergyK*phycon.teinv;
 
         if( c == 0. )
-        {
                 c = a + b;
-        }
 
         /* if still neg at end, then success!, so possible to
          * to check why zero returned */
@@ -209 +213 @@
         /* if two of the lines are below the plasma frequency, hand the third in atom_level2 */
         if( ( t10->EnergyErg/EN1RYD < rfield.plsfrq && t20->EnergyErg/EN1RYD < rfield.plsfrq ) )
         {
+                // update RT quantities to keep them synched with level populations
+                if( conv.nTotalIoniz && !conv.lgLastSweepThisZone )
+                        RT_line_one( t21, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
                 atom_level2( t21 );
                 return;
         }
         else if( ( t10->EnergyErg/EN1RYD < rfield.plsfrq && t21->EnergyErg/EN1RYD < rfield.plsfrq ) )
         {
+                // update RT quantities to keep them synched with level populations
+                if( conv.nTotalIoniz && !conv.lgLastSweepThisZone )
+                        RT_line_one( t20, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
                 atom_level2( t20 );
                 return;
         }
         else if( ( t20->EnergyErg/EN1RYD < rfield.plsfrq && t21->EnergyErg/EN1RYD < rfield.plsfrq ) )
         {
+                // update RT quantities to keep them synched with level populations
+                if( conv.nTotalIoniz && !conv.lgLastSweepThisZone )
+                        RT_line_one( t10, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
                 atom_level2( t10 );
                 return;
         }
@@ -314 +327 @@
         ASSERT( (t20->Hi->nelem * t21->Hi->nelem == 0) || (t20->Hi->nelem == t21->Hi->nelem) );
 
         ASSERT( o10 > 0. && o21 > 0. && o20 > 0. );
-
         /*end sanity checks */
 
         /* net loss of line escape and destruction */
@@ -335 +347 @@
                 temp12 = t21->EnergyK;
                 temp02 = t20->EnergyK;
                 temp01 = temp02 - temp12;
+
+                // update RT quantities to be in sync with level populations
+                RT_line_one( t21, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                RT_line_one( t20, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
         }
 
         else if( a21 == 0. )
@@ -348 +364 @@
                 temp02 = t20->EnergyK;
                 temp01 = t10->EnergyK;
                 temp12 = temp02 - temp01;
+                // update RT quantities to be in sync with level populations
+                if( conv.nTotalIoniz && !conv.lgLastSweepThisZone )
+                {
+                        RT_line_one( t10, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                        RT_line_one( t20, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                }
         }
 
         else if( a20 == 0. )
@@ -361 +383 @@
                 temp01 = t10->EnergyK;
                 temp12 = t21->EnergyK;
                 temp02 = temp01 + temp12;
+                // update RT quantities to be in sync with level populations
+                if( conv.nTotalIoniz && !conv.lgLastSweepThisZone )
+                {
+                        RT_line_one( t21, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                        RT_line_one( t10, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                }
         }
 
         else
@@ -375 +403 @@
                 temp02 = t20->EnergyK;
                 temp01 = t10->EnergyK;
                 temp12 = t21->EnergyK;
+                // update RT quantities to be in sync with level populations
+                if( conv.nTotalIoniz && !conv.lgLastSweepThisZone )
+                {
+                        RT_line_one( t21, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                        RT_line_one( t10, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                        RT_line_one( t20, true,0.f, GetDopplerWidth(dense.AtomicWeight[nelem-1]) );
+                }
         }
 
         /* check all energies positive */

cool_dima.cpp

@@ -53 +53 @@
                         }
                         /* now put the cs into the line array */
                         PutCS(cs,&TauLine2[i] );
-                        RT_line_one( &TauLine2[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[TauLine2[i].Hi->nelem-1]) );
                         atom_level2(&TauLine2[i] );
                 }
         }