Changeset 2034 for branches/c08_branch/source/iso_cool.cpp

Timestamp:

05/10/08 09:03:02 (6 months ago)

Author:

peter

Message:

Merging all changes from mainline upto r2033, except r1902.

Files:

• branches/c08_branch/source/iso_cool.cpp (modified) (13 diffs)

branches/c08_branch/source/iso_cool.cpp

@@ -12 +12 @@
 #include "cooling.h"
 #include "iso.h"
-/*lint -e661 Possible access of out-of-bounds pointer*/
-/*lint -e662  (Warning -- Possible creation of out-of-bounds pointer  */
 
 /* HP cc cannot compile this routine with any optimization */ 
@@ -61 +59 @@
         char chLabel[NCOLNT_LAB_LEN+1];
 
-        /* these are labels for iso line cooling, sum of ALL hydrogen lines,
-         * along the iso sequence
-        static char chHlin[LIMELM][5]=
-        {"Hl 0","Hl 1","Hl 2","Hl 3","Hl 4","Hl 5","Hl 6","Hl 7","Hl 8","Hl 9",
-         "Hl10","Hl11","Hl12","Hl13","Hl14","Hl15","Hl16","Hl17","Hl18","Hl19",
-         "Hl20","Hl21","Hl22","Hl23","Hl24","Hl25","Hl26","Hl27","Hl28","Hl29"}; */
-
-        /* these are labels for iso collisional ionization cooling,
-         * along the iso sequence
-        static char chHion[LIMELM][5]=
-        {"Hi 0","Hi 1","Hi 2","Hi 3","Hi 4","Hi 5","Hi 6","Hi 7","Hi 8","Hi 9",
-         "Hi10","Hi11","Hi12","Hi13","Hi14","Hi15","Hi16","Hi17","Hi18","Hi19",
-         "Hi20","Hi21","Hi22","Hi23","Hi24","Hi25","Hi26","Hi27","Hi28","Hi29"}; */
-
-        /* these are labels for iso photoionization heating,
-         * along the iso sequence 
-        static char chHp[LIMELM][5]=
-        {"Hp 0","Hp 1","Hp 2","Hp 3","Hp 4","Hp 5","Hp 6","Hp 7","Hp 8","Hp 9",
-         "Hp10","Hp11","Hp12","Hp13","Hp14","Hp15","Hp16","Hp17","Hp18","Hp19",
-         "Hp20","Hp21","Hp22","Hp23","Hp24","Hp25","Hp26","Hp27","Hp28","Hp29"};*/
-
-        /* these are labels for iso induced recombination heating,
-         * along the iso sequence
-        static char chHn[LIMELM][5]=
-        {"Hn 0","Hn 1","Hn 2","Hn 3","Hn 4","Hn 5","Hn 6","Hn 7","Hn 8","Hn 9",
-         "Hn10","Hn11","Hn12","Hn13","Hn14","Hn15","Hn16","Hn17","Hn18","Hn19",
-         "Hn20","Hn21","Hn22","Hn23","Hn24","Hn25","Hn26","Hn27","Hn28","Hn29"}; */
-
         DEBUG_ENTRY( "iso_cool()" );
 
@@ -160 +130 @@
          * and its coupling to continuum can be large, at various times code 
          * had to ignore effects of very highest level, but starting in mid
-         * 20006 all levels have been included */
-        for( n=0; n < iso.numLevels_local[ipISO][nelem]; ++n )
+         * 20006 all levels have been included 
+         * 2008 apr 18, do not include highest - when only 2 collapsed levels
+          * are used several density 13 BLR sims have serious convergence
+          * problems */
+        for( n=0; n < iso.numLevels_local[ipISO][nelem]-1; ++n )
         {
                 /* total collisional ionization cooling less three body heating */
@@ -168 +141 @@
                   (StatesElem[ipISO][nelem][n].Pop -iso.PopLTE[ipISO][nelem][n]*dense.eden);
                 CollisIonizatCoolingTotal += CollisIonizatCooling;
-                /*fprintf(ioQQQ,"DEBUG CollisIonizatCoolingTotal %.2e %2li %2li %2li %.2e %.4e %.4e %.3f\n",
-                        CollisIonizatCoolingTotal,
-                        ipISO , nelem , n ,
-                        CollisIonizatCooling , 
-                        StatesElem[ipISO][nelem][n].Pop ,
-                        iso.PopLTE[ipISO][nelem][n]*dense.eden ,
-                        fabs(StatesElem[ipISO][nelem][n].Pop - iso.PopLTE[ipISO][nelem][n]*dense.eden ) / 
-                            (iso.PopLTE[ipISO][nelem][n]*dense.eden) );*/
 
                 /* the derivative of the cooling 
                  * need extra factor of temp of 1 Ryd since div by square of temp in Ryd */
                 CollisIonizatCoolingDT = CollisIonizatCooling*
-                        (iso.xIsoLevNIonRyd[ipISO][nelem][n]/TE1RYD/POW2(phycon.te_ryd)/* - 
-                  thermal.halfte*/);
+                        (iso.xIsoLevNIonRyd[ipISO][nelem][n]/TE1RYD/POW2(phycon.te_ryd));
+
                 dCollisIonizatCoolingTotalDT += CollisIonizatCoolingDT;
                 // save values for debug printout
@@ -280 +245 @@
                  thin = iso.RadRecomb[ipISO][nelem][0][ipRecRad]*
                         /* arg is the scaled temperature, T * n^2 / Z^2, 
-                         * n is principle quantum number, Z is charge, 1 for H */
+                         * n is principal quantum number, Z is charge, 1 for H */
                         HCoolRatio( 
                         phycon.te * POW2( (double)StatesElem[ipISO][nelem][0].n / (double)(nelem+1-ipISO) ))*
@@ -297 +262 @@
                  thin = iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
                         /* arg is the scaled temperature, T * n^2 / Z^2, 
-                         * n is principle quantum number, Z is charge, 1 for H */
+                         * n is principal quantum number, Z is charge, 1 for H */
                         HCoolRatio( 
                         phycon.te * POW2( (double)StatesElem[ipISO][nelem][n].n / (double)(nelem+1-ipISO) ))*
@@ -316 +281 @@
                 if( DEBUG_LOC  )
                 {
-                        if( nelem==1 && ipISO==ipHE_LIKE)
+                        if( nelem==ipISO )
                         {
                                 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
@@ -526 +491 @@
 
                         iso.xLineTotCool[ipISO][nelem] += hlone;
-                        {
-                                /* debug block for collisional cooling between two levels */
-                                enum {DEBUG_LOC=false};
-                                if( DEBUG_LOC && ipISO==ipHE_LIKE && nelem==ipHELIUM /**/ )
-                                {
-                                        fprintf(ioQQQ,"iso_coolll %li %li %.2e %.2e %.2e %.2e\n",
-                                                ipHi , ipLo, 
-                                                StatesElem[ipISO][nelem][ipHi].Pop ,
-                                                Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL , 
-                                                Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs , 
-                                                hlone/SDIV( thermal.htot ) );
-                                }
-                        }
 
                         /* next get derivative */
@@ -569 +521 @@
                 if( DEBUG_LOC )
                 {
-                        if( ipISO==ipHE_LIKE && nelem==ipHELIUM )
+                        if( nelem==ipISO )
                                 fprintf(ioQQQ,"%li %li %.2e\n", nlo_heat_max, nhi_heat_max, heat_max );
                 }
@@ -655 +607 @@
         if( iso.xLineTotCool[ipISO][nelem] > 0. )
         {
-                /* hydrogen is a net coolant label gives iso sequence, "wavelength" gives element */
+                /* species is a net coolant label gives iso sequence, "wavelength" gives element */
                 CoolAdd(chLabel,(realnum)nelem,iso.xLineTotCool[ipISO][nelem]);
                 thermal.dCooldT += dCdT_all;
@@ -662 +614 @@
         else
         {
-                /* hydrogen is a net heat source, thermal.heating[0][23]was set to 0 in CoolEvaluate*/
+                /* species is a net heat source, thermal.heating[0][23]was set to 0 in CoolEvaluate*/
                 thermal.heating[0][23] -= iso.xLineTotCool[ipISO][nelem];
                 CoolAdd(chLabel,(realnum)nelem,0.);
@@ -681 +633 @@
                                         iso.cBal_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ,
                                         iso.cRest_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] );
-
-#                               if 0
-                                if( iso.cLya_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] > 0.9 )
-                                fprintf(ioQQQ," 1s 2s 2p %.2e %.2e %.2e\n",
-                                        StatesElem[ipISO][nelem][ipH1s].Pop , 
-                                        StatesElem[ipISO][nelem][ipH2s].Pop ,
-                                        StatesElem[ipISO][nelem][ipH2p].Pop);
-
-#                               endif
                         }
                 }
@@ -756 +699 @@
 #pragma OPTIMIZE ON
 #endif
-/*lint +e661 Possible access of out-of-bounds pointer*/
-/*lint +e662 (Warning -- Possible creation of out-of-bounds pointer  */