Changeset 1179

Timestamp:

06/09/07 13:18:05 (18 months ago)

Author:

gary

Message:

most files change because hard lower and upper limits to range of allowed temperatures changed to macro set in cddefines.h - this defines the following:
const double TEMP_STOP_DEFAULT = 4000.;
/** highest temperature to ever allow */
const double TEMP_LIMIT_HIGH = 1e10;
/** lowest temperature to ever allow */
const double TEMP_LIMIT_LOW = 2.8;

optimize_func.cpp - do not gather optimization quality metric in case where grid is being run. some grid models may abort due to T < 2.8K in which case it is not possible to gather this info anyway.

Location:

trunk/source

Files:

• abundances.cpp (modified) (1 diff)
• cddefines.h (modified) (1 diff)
• conv_init_solution.cpp (modified) (4 diffs)
• conv_temp_eden_ioniz.cpp (modified) (2 diffs)
• date.h (modified) (1 diff)
• init_defaults_preparse.cpp (modified) (2 diffs)
• map.cpp (modified) (3 diffs)
• mole_h2.cpp (modified) (1 diff)
• optimize_func.cpp (modified) (1 diff)
• parse_commands.cpp (modified) (1 diff)
• parse_constant.cpp (modified) (2 diffs)
• parse_coronal.cpp (modified) (2 diffs)
• parse_stop.cpp (modified) (1 diff)
• phycon.h (modified) (1 diff)
• prt_comment.cpp (modified) (1 diff)
• stopcalc.h (modified) (5 diffs)

trunk/source/abundances.cpp

@@ -223 +223 @@
         /* if stop temp set below default then we are going into cold and possibly molecular
          * gas - check some parameters in this case */
-        if( called.lgTalk && (StopCalc.tend < TEND || 
+        if( called.lgTalk && (StopCalc.tend < TEMP_STOP_DEFAULT || 
                 /* thermal.ConstTemp def is zero, set pos when used */
-                (thermal.ConstTemp > 0. && thermal.ConstTemp < TEND ) ) )
+                (thermal.ConstTemp > 0. && thermal.ConstTemp < TEMP_STOP_DEFAULT ) ) )
         {
 

trunk/source/cddefines.h

@@ -175 +175 @@
 const int NHYDRO_MAX_LEVEL = 401;
 
-/** dimension of vector to do with number of iterations stored
- * so, this is the the most that can possibly be performed */
-/*const int ITR DIM = 200;*/
-
 /** number of hydrogen species in hmole */
 const int N_H_MOLEC = 8;
 
-/** stopping temperature */
-const double TEND = 4000.;
+/** The default value of the stopping temperature */
+const double TEMP_STOP_DEFAULT = 4000.;
+/** highest temperature to ever allow */
+const double TEMP_LIMIT_HIGH = 1e10;
+/** lowest temperature to ever allow */
+const double TEMP_LIMIT_LOW = 2.8;
 
 /** this is used to add to depth to prevent div or log of zero */

trunk/source/conv_init_solution.cpp

@@ -4 +4 @@
 #include "cddefines.h"
 #include "trace.h"
-#include "stopcalc.h"
 #include "struc.h"
 #include "rfield.h"
@@ -328 +327 @@
                                   phycon.te, thermal.htot, thermal.ctot, Cool_min_Heat );
                         }
-                        /* check range 10-4000K, TeLowest usually 2.8K, set with SET LOWEST command
+                        /* check range 10-4000K, TEMP_LIMIT_LOW is lowest possible temperature
                          * following are series of tests that determine the factor by which
                          * which can change the temperature.  must be VERY small when 
                          * ice formation on grains is important due to exponential sublimation
                          * rate.  Also small when molecules are important, to keep chemistry
-                         * within limits of linearization solver */
+                         * within limits of linearized solver */
                         fac2 = 0.8;
-                        while( phycon.te*fac2 > StopCalc.TeLowest && Cool_min_Heat >= 0. && !lgAbort )
+                        while( phycon.te*fac2 > TEMP_LIMIT_LOW && Cool_min_Heat >= 0. && !lgAbort )
                         {
                                 double oxy_in_grains;
@@ -474 +473 @@
                                 ShowMe();
                                 
+                                /* this is an error return, calculation will immediately stop */
                                 DEBUG_EXIT( "ConvInitSolution()" );
+                                return(1);
+                        }
+                        else if( thermal.ctot > thermal.htot )
+                        {
+                                /* even at lowest temp, too much cooling; no solution possible */
+                                fprintf( ioQQQ, 
+                                        "PROBLEM The initial kinetic temperature was below %.1e K,"
+                                        " the lowest possible temperature I can consider."
+                                        "\nSorry, but I cannot go on.\n\n", 
+                                  TEMP_LIMIT_LOW );
+
                                 /* this is an error return, calculation will immediately stop */
-                                return(1);
-                        }
-                        else if( thermal.ctot > thermal.htot )
-                        {
-                                /* even at lowest temp, too much cooling; no solution possible */
-                                fprintf( ioQQQ, " Equilibrium temperature below T=%8.2e.  If this is OK, then reset with LOWEST command.  This is ConvInitSolution.  TeLowest was%9.2e\n", 
-                                  phycon.te, StopCalc.TeLowest );
-                                
                                 DEBUG_EXIT( "ConvInitSolution()" );
-                                /* this is an error return, calculation will immediately stop */
                                 return(1);
                         }
@@ -505 +507 @@
 
                 /* coming up to here te is either 4,000 (usually) or 10^6 */
-                t1 = StopCalc.TeLowest/1.001;
+                t1 = TEMP_LIMIT_LOW/1.001;
 
                 /* >>chng 96 dec 28, added div by factor to prevent checking on temp > 1e10 */
-                t2 = StopCalc.TeHighest*1.001/factor;
+                t2 = TEMP_LIMIT_HIGH*1.001/factor;
                 if( trace.nTrConvg )
                 {

trunk/source/conv_temp_eden_ioniz.cpp

@@ -14 +14 @@
 #include "colden.h"
 #include "h2.h"
-#include "stopcalc.h"
 #include "pressure.h"
 #include "dense.h"
@@ -943 +942 @@
 {
         bool lgRet;
-        /* >>chng 03 sep 19, add check on StopCalc.TeLowest */
+        /* >>chng 03 sep 19, add check on TEMP_LIMIT_LOW */
         if( fabs(thermal.htot - thermal.ctot)/thermal.htot <= conv.HeatCoolRelErrorAllowed ||
-                thermal.lgTSetOn || phycon.te <= StopCalc.TeLowest )
+                thermal.lgTSetOn || phycon.te <= TEMP_LIMIT_LOW )
         {
                 /* announce that temp is converged if relative heating - cooling mismatch

trunk/source/date.h

@@ -9 +9 @@
 #define YEAR    107
 #define MONTH   5
-#define DAY     7
+#define DAY     9

trunk/source/init_defaults_preparse.cpp

@@ -143 +143 @@
         StopCalc.AV_point = 1e30f;
         /* highest allowed temperature */
-        StopCalc.T2High = 1e12f;
-
-        /* highest and lowest temperatures to ever allow */
-        StopCalc.TeHighest = 1e10f;
-        StopCalc.TeLowest = 2.8f;
+        StopCalc.T2High = (float)TEMP_LIMIT_HIGH;
 
         /* the floor sets a limit to the temperature in the calculation -
@@ -155 +151 @@
 
         /* this is in cddefines.h and is 4000 */
-        StopCalc.tend = (float)TEND;
+        StopCalc.tend = (float)TEMP_STOP_DEFAULT;
         /* ending column densities */
         StopCalc.HColStop = COLUMN_INIT;

trunk/source/map.cpp

@@ -6 +6 @@
 #include "thermal.h"
 #include "cooling.h"
-#include "stopcalc.h"
 #include "called.h"
 #include "dense.h"
@@ -125 +124 @@
         {
                 /* create some sort of map of heating-cooling */
-                tlowst = MAX2(map.RangeMap[0],StopCalc.TeLowest);
+                tlowst = MAX2(map.RangeMap[0],TEMP_LIMIT_LOW);
                 tmin = tlowst*0.998;
-                tmax = MIN2(map.RangeMap[1],StopCalc.TeHighest)*1.002;
+                tmax = MIN2(map.RangeMap[1],TEMP_LIMIT_HIGH)*1.002;
 
                 /* we want about nMapStep (=20) steps between high and low temperature */
@@ -137 +136 @@
                         factor = 1./factor;
                         /* TeHighest is highest possible temperature, 1E10 */
-                        phycon.te = (MIN2(map.RangeMap[1],StopCalc.TeHighest)/factor);
+                        phycon.te = (MIN2(map.RangeMap[1],TEMP_LIMIT_HIGH)/factor);
                 }
 

trunk/source/mole_h2.cpp

@@ -2831 +2831 @@
                         conv.HeatCoolRelErrorAllowed/2.
                         /* >>chng 04 may 09, do not check on error in heating if constant temperature */
-                        /*&& !(thermal.lgTSetOn || phycon.te <= StopCalc.TeLowest  )*/ )
+                        /*&& !(thermal.lgTSetOn || phycon.te <= TEMP_LIMIT_LOW  )*/ )
                 {
                         /* default on HeatCoolRelErrorAllowed is 0.02 */

trunk/source/optimize_func.cpp

@@ -130 +130 @@
                 GridInitialize();
         }
-
-        /* Line fluxes now in commn blocks, so extract and
-         * compare with observations */
-        chisq = 0.0;
-        for( i=0; i < MAXCAT; i++ )
-        {
-                nobs_cat[i] = 0;
-                chi2_cat[i] = 0.0;
-        }
-
-        if( LineSave.ipNormWavL < 0 )
-        {
-                fprintf( ioQQQ, 
-                        " Normalization line array index is bad.  What has gone wrong?\n" );
-                puts( "[Stop in optimize_func]" );
-                cdEXIT(EXIT_FAILURE);
-        }
-        snorm = LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent];
-
-        if( snorm == 0. )
-        {
-                fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity.  What has gone wrong?\n" );
-                fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
-                puts( "[Stop in optimize_func]" );
-                cdEXIT(EXIT_FAILURE);
-        }
-
-        /* first print all warnings */
-        cdWarnings(ioQQQ);
-
-        /* cycle through the observational values */
-        nfound = 0;
-        
-        /* first is to optimize relative emission line spectrum */
-        if( optimize.lgOptLin )
-        {
-                /* set pointers to all optimized lines if first call */
-                if( !lgLinSet )
-                {
-                        lgHIT = true;
-                        lgLinSet = true;
+        else
+        {
+                /* this branch optimizing, not grid 
+                / * extract line fluxes and compare with observations */
+                chisq = 0.0;
+                for( i=0; i < MAXCAT; i++ )
+                {
+                        nobs_cat[i] = 0;
+                        chi2_cat[i] = 0.0;
+                }
+
+                if( LineSave.ipNormWavL < 0 )
+                {
+                        fprintf( ioQQQ, 
+                                " Normalization line array index is bad.  What has gone wrong?\n" );
+                        puts( "[Stop in optimize_func]" );
+                        cdEXIT(EXIT_FAILURE);
+                }
+
+                if( (snorm = LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) == 0. )
+                {
+                        fprintf( ioQQQ, "\n\n PROBLEM Normalization line has zero intensity.  What has gone wrong?\n" );
+                        fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
+                        puts( "[Stop in optimize_func]" );
+                        cdEXIT(EXIT_FAILURE);
+                }
+
+                /* first print all warnings */
+                cdWarnings(ioQQQ);
+
+                /* cycle through the observational values */
+                nfound = 0;
+                
+                /* first is to optimize relative emission line spectrum */
+                if( optimize.lgOptLin )
+                {
+                        /* set pointers to all optimized lines if first call */
+                        if( !lgLinSet )
+                        {
+                                lgHIT = true;
+                                lgLinSet = true;
+                                for( i=0; i < optimize.nlobs; i++ )
+                                {
+                                        double temp1, temp2;
+                                        cap4(chFind , (char*)optimize.chLineLabel[i]);
+                                        /* j = 0; */
+
+                                        /* >> chng 06 may 04, use cdLine instead of ad hoc treatment.
+                                         * no need to complain, cdLine will do it automatically.  */
+                                        /* this is an intensity, get the line, returns false if could not find it */
+                                        j=cdLine( chFind , optimize.wavelength[i] , &temp1, &temp2 );
+                                        if( j<=0 )
+                                        {
+                                                fprintf( ioQQQ, "\n" );
+                                                lgHIT = false;
+                                        }
+                                        else
+                                        {
+                                                ipobs[i] = j;
+                                        }
+                                }
+
+                                /* we did not find the line */
+                                if( !lgHIT )
+                                {
+                                        fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" );
+                                        fprintf( ioQQQ, " Sorry.\n");
+                                        puts( "[Stop in optimize_func]" );
+                                        cdEXIT(EXIT_FAILURE);
+                                }
+                        }
+
+                        for( i=0; i < 10; i++ )
+                        {
+                                optimize.SavGenericData[i] = 0.;
+                        }
+
                         for( i=0; i < optimize.nlobs; i++ )
                         {
-                                double temp1, temp2;
-                                cap4(chFind , (char*)optimize.chLineLabel[i]);
-                                /* j = 0; */
-
-                                /* >> chng 06 may 04, use cdLine instead of ad hoc treatment.
-                                 * no need to complain, cdLine will do it automatically.  */
-                                /* this is an intensity, get the line, returns false if could not find it */
-                                j=cdLine( chFind , optimize.wavelength[i] , &temp1, &temp2 );
-                                if( j<=0 )
+                                /* and find corresponding model value by straight search */
+                                nfound += 1;
+                                scld = (double)LineSv[ipobs[i]].sumlin[LineSave.lgLineEmergent]/(double)snorm*LineSave.ScaleNormLine;
+                                chi1 = chi2_func(scld,(double)optimize.xLineInt_Obs[i],(double)optimize.xLineInt_error[i]);
+                                cat = 0;
+                                nobs_cat[cat]++;
+                                chi2_cat[cat] += chi1;
+
+                                fprintf( ioQQQ, " %4.4s ", 
+                                  LineSv[ipobs[i]].chALab);
+
+                                prt_wl( ioQQQ,LineSv[ipobs[i]].wavelength);
+
+                                fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity", 
+                                  scld, 
+                                  optimize.xLineInt_Obs[i], 
+                                  optimize.xLineInt_error[i], 
+                                  chi1 );
+                                
+                                fprintf( ioQQQ, "\n" );
+
+                                if( i<10 )
                                 {
-                                        fprintf( ioQQQ, "\n" );
-                                        lgHIT = false;
+                                        optimize.SavGenericData[i] = chi1;
                                 }
-                                else
+                        }
+                }
+
+                /* >>chng 06 may 04, moved this from above so that it only
+                 * gets printed if all lines are found */
+                /*if( optimize.lgOptLum || optimize.lgOptCol || optimize.lgOptTemp || optimize.lgOptLin )*/
+                if( optimize.lgOptimize )
+                {
+                        fprintf( ioQQQ, "  ID               Model    Observed       error      chi**2     Type\n" );
+                }
+                else
+                {
+                        ASSERT( grid.lgGrid );
+                }
+
+                /* this is to optimize a mean temperature */
+                if( optimize.lgOptTemp )
+                {
+                        for( i=0; i < optimize.nTempObs; i++ )
+                        {
+                                if( cdTemp(/*(char*)*/optimize.chTempLab[i],optimize.ionTemp[i], &temp_theory, optimize.chTempWeight[i]) )
                                 {
-                                        ipobs[i] = j;
+                                        /* did not find column density */
+                                        fprintf(ioQQQ," optimizer did not find column density %s %li \n",
+                                                optimize.chTempLab[i],optimize.ionTemp[i] );
+                                        puts( "[Stop in optimize_func]" );
+                                        cdEXIT(EXIT_FAILURE);
                                 }
-                        }
-
-                        /* we did not find the line */
-                        if( !lgHIT )
-                        {
-                                fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" );
-                                fprintf( ioQQQ, " Sorry.\n");
-                                puts( "[Stop in optimize_func]" );
-                                cdEXIT(EXIT_FAILURE);
-                        }
-                }
-
-                for( i=0; i < 10; i++ )
-                {
-                        optimize.SavGenericData[i] = 0.;
-                }
-
-                for( i=0; i < optimize.nlobs; i++ )
-                {
-                        /* and find corresponding model value by straight search */
+                                nfound += 1;
+                                chi1 = chi2_func(temp_theory,(double)optimize.temp_obs[i],(double)optimize.temp_error[i]);
+                                cat = 3;
+                                nobs_cat[cat]++;
+                                chi2_cat[cat] += chi1;
+
+                                fprintf( ioQQQ, " %4.4s %2ld        ",
+                                  optimize.chTempLab[i], 
+                                  optimize.ionTemp[i] );
+                                PrintE82( ioQQQ, temp_theory );
+                                fprintf( ioQQQ, "    ");
+                                PrintE82( ioQQQ, optimize.temp_obs[i] );
+                                fprintf( ioQQQ, "    %.5f   %.2e", 
+                                        optimize.temp_error[i],  chi1 );
+                                fprintf( ioQQQ, " Temperature\n");
+                        }
+                }
+
+                /* option to optimize column densities */
+                if( optimize.lgOptCol )
+                {
+                        for( i=0; i < optimize.ncobs; i++ )
+                        {
+                                if( cdColm((char*)optimize.chColDen_label[i],optimize.ion_ColDen[i], &theocl) )
+                                {
+                                        /* did not find column density */
+                                        fprintf(ioQQQ," optimizer did not find column density %s %li \n",
+                                                optimize.chColDen_label[i],optimize.ion_ColDen[i] );
+                                        puts( "[Stop in optimize_func]" );
+                                        cdEXIT(EXIT_FAILURE);
+                                }
+                                nfound += 1;
+                                chi1 = chi2_func(theocl,(double)optimize.ColDen_Obs[i],(double)optimize.chColDen_error[i]);
+                                cat = 1;
+                                nobs_cat[cat]++;
+                                chi2_cat[cat] += chi1;
+
+                                fprintf( ioQQQ, " %4.4s%6ld%12.4e%12.4e%12.5f%12.2e Column density\n", 
+                                  optimize.chColDen_label[i], optimize.ion_ColDen[i], theocl, 
+                                  optimize.ColDen_Obs[i], optimize.chColDen_error[i], chi1 );
+                        }
+                }
+
+                /* option to optimize line flux */
+                if( optimize.lgOptLum )
+                {
                         nfound += 1;
-                        scld = (double)LineSv[ipobs[i]].sumlin[LineSave.lgLineEmergent]/(double)snorm*LineSave.ScaleNormLine;
-                        chi1 = chi2_func(scld,(double)optimize.xLineInt_Obs[i],(double)optimize.xLineInt_error[i]);
-                        cat = 0;
+                        if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0.f )
+                        {
+                                predin = log10(LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten;
+                                help = pow(10.,predin-(double)optimize.optint);
+                                chi1 = chi2_func(help,1.,(double)optimize.optier);
+                        }
+                        else
+                        {
+                                predin = -999.99999;
+                                chi1 = (double)FLT_MAX;
+                        }
+                        cat = 2;
                         nobs_cat[cat]++;
                         chi2_cat[cat] += chi1;
 
-                        fprintf( ioQQQ, " %4.4s ", 
-                          LineSv[ipobs[i]].chALab);
-
-                        prt_wl( ioQQQ,LineSv[ipobs[i]].wavelength);
-
-                        fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity", 
-                          scld, 
-                          optimize.xLineInt_Obs[i], 
-                          optimize.xLineInt_error[i], 
+                        fprintf( ioQQQ, " %4.4s%6f%12.5f%12.5f%12.5f%12.2e Line intensity\n", 
+                          LineSv[LineSave.ipNormWavL].chALab, 
+                          LineSv[LineSave.ipNormWavL].wavelength, 
+                          predin, 
+                          optimize.optint, 
+                          optimize.optier, 
                           chi1 );
-                        
-                        fprintf( ioQQQ, "\n" );
-
-                        if( i<10 )
-                        {
-                                optimize.SavGenericData[i] = chi1;
-                        }
-                }
-        }
-
-        /* >>chng 06 may 04, moved this from above so that it only
-         * gets printed if all lines are found */
-        if( optimize.lgOptLum || optimize.lgOptCol || optimize.lgOptTemp || optimize.lgOptLin )
-        {
-                fprintf( ioQQQ, "  ID               Model    Observed       error      chi**2     Type\n" );
-        }
-        else
-        {
-                ASSERT( grid.lgGrid );
-        }
-
-        /* this is to optimize a mean temperature */
-        if( optimize.lgOptTemp )
-        {
-                for( i=0; i < optimize.nTempObs; i++ )
-                {
-                        if( cdTemp(/*(char*)*/optimize.chTempLab[i],optimize.ionTemp[i], &temp_theory, optimize.chTempWeight[i]) )
-                        {
-                                /* did not find column density */
-                                fprintf(ioQQQ," optimizer did not find column density %s %li \n",
-                                        optimize.chTempLab[i],optimize.ionTemp[i] );
-                                puts( "[Stop in optimize_func]" );
-                                cdEXIT(EXIT_FAILURE);
-                        }
-                        nfound += 1;
-                        chi1 = chi2_func(temp_theory,(double)optimize.temp_obs[i],(double)optimize.temp_error[i]);
-                        cat = 3;
-                        nobs_cat[cat]++;
-                        chi2_cat[cat] += chi1;
-
-                        fprintf( ioQQQ, " %4.4s %2ld        ",
-                          optimize.chTempLab[i], 
-