Changeset 1751

Timestamp:

01/22/08 13:18:24 (10 months ago)

Author:

rporter

Message:

trunk/source/grid_do.cpp - minor changes in grid/optimizer printout to make difference more clear.

trunk/source/iso_continuum_lower.cpp - fix bug when continuum was lowered in previous zone and then is not lowered.

trunk/source/iso_create.cpp - change lifetime calculation for upper level of extra lyman lines. Dramatically reduces startup time in fast models.

trunk/source/optimize_func.cpp - minor printout change

trunk/source/parse.h - declare ParseGrid?

trunk/source/parse_commands.cpp - call ParseGrid? instead of ParseOptimize? on GRID keyword

trunk/source/parse_grid.cpp - new file, extracted from ParseOptimize?

trunk/source/parse_optimize.cpp - remove GRID keyword as above

trunk/source/punch_fits.cpp - add fixit about FITS files problem with all double precision

trunk/tsuite/auto/dynamics_wind.in - increase itrzn - this has apparently been botching for a month

Location:

trunk

Files:

• source/grid_do.cpp (modified) (4 diffs)
• source/iso_continuum_lower.cpp (modified) (3 diffs)
• source/iso_create.cpp (modified) (4 diffs)
• source/optimize_func.cpp (modified) (1 diff)
• source/parse.h (modified) (1 diff)
• source/parse_commands.cpp (modified) (1 diff)
• source/parse_grid.cpp (added)
• source/parse_optimize.cpp (modified) (1 diff)
• source/punch_fits.cpp (modified) (1 diff)
• tsuite/auto/dynamics_wind.in (modified) (1 diff)

trunk/source/grid_do.cpp

@@ -237 +237 @@
         }
 
-        fprintf( ioQQQ, " Up to%5ld iterations will be performed,\n", 
-          optimize.nIterOptim );
-        fprintf( ioQQQ, " and the final version of the input file will be written to the file %s\n", 
-          chOptimFileName );
 
         if( strcmp(optimize.chOptRtn,"PHYM") == 0 )
         {
+                fprintf( ioQQQ, " Up to%5ld iterations will be performed,\n", 
+                  optimize.nIterOptim );
+                fprintf( ioQQQ, " and the final version of the input file will be written to the file %s\n", 
+                  chOptimFileName );
+
                 fprintf( ioQQQ, " The optimize_phymir method will be used" );
                 if( optimize.lgParallel ) {
@@ -255 +256 @@
         else if( strcmp(optimize.chOptRtn,"SUBP") == 0 )
         {
+                fprintf( ioQQQ, " Up to%5ld iterations will be performed,\n", 
+                  optimize.nIterOptim );
+                fprintf( ioQQQ, " and the final version of the input file will be written to the file %s\n", 
+                  chOptimFileName );
+
                 fprintf( ioQQQ, " Subplex method will be used.\n" );
         }
@@ -260 +266 @@
         else if( strcmp(optimize.chOptRtn,"XSPE") == 0 )
         {
-                fprintf( ioQQQ, " Producing XSPEC output.\n" );
+                fprintf( ioQQQ, " Producing grid output.\n" );
         }
 
@@ -267 +273 @@
                 fprintf( ioQQQ, " I do not understand what method to use.\n" );
                 fprintf( ioQQQ, " Sorry.\n" );
-                puts( "[Stop in lgOptimize_do]" );
+                puts( "[Stop in grid_do]" );
                 cdEXIT(EXIT_FAILURE);
         }
 
-        fprintf( ioQQQ, "\n%4ld parameters will be varied.  The first lines, and the increments are:\n", 
+        fprintf( ioQQQ, "\n%4ld parameter(s) will be varied.  The first lines, and the increments are:\n", 
           optimize.nvary );
 

trunk/source/iso_continuum_lower.cpp

@@ -79 +79 @@
         {
                 iso.lgLevelsLowered[ipISO][nelem] = true;
+                iso.lgLevelsEverLowered[ipISO][nelem] = true;
                 hydro.lgReevalRecom = true;
                 iso.n_HighestResolved_local[ipISO][nelem] = nc;
@@ -94 +95 @@
         {
                 iso.lgLevelsLowered[ipISO][nelem] = true;
+                iso.lgLevelsEverLowered[ipISO][nelem] = true;
                 hydro.lgReevalRecom = true;
                 iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem];
@@ -118 +120 @@
                 iso.nCollapsed_local[ipISO][nelem] = iso.nCollapsed_max[ipISO][nelem];
                 iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem];
+                
+                /* if levels were lowered on last past but are not now, must reeval */
+                if( iso.lgLevelsLowered[ipISO][nelem] )
+                {
+                        hydro.lgReevalRecom = true;
+                }
+                else
+                {
+                        hydro.lgReevalRecom = false;
+                }
+
                 iso.lgLevelsLowered[ipISO][nelem] = false;
-                hydro.lgReevalRecom = false;
         }
 

trunk/source/iso_create.cpp

@@ -32 +32 @@
 STATIC void FillExtraLymanLine( transition *t, long ipISO, long nelem, long nHi );
 
+/* calculate radiative lifetime of an individual iso state */
+STATIC double iso_state_lifetime( long ipISO, long nelem, long n, long l );
+
 /* calculate cascade probabilities, branching ratios, and associated errors. */
 STATIC void iso_cascade(void);
@@ -1236 +1239 @@
 {
         double Enerwn, Aul;
-        long nLo;
 
         DEBUG_ENTRY( "FillExtraLymanLine()" );
@@ -1279 +1281 @@
         t->Emis->Aul = (realnum)Aul;
 
-        t->Hi->lifetime = t->Emis->Aul;
-
-        /* Now loop over lower levels to get sum of As out, and invert for lifetime. */
-        for( nLo=2; nLo < nHi; nLo++ )
-        {
-                double ea;
-
-                ea = H_Einstein_A( nHi, 1, nLo, 0, nelem+1-ipISO );
-
-                ASSERT( ea > 0. );
-                t->Hi->lifetime += (realnum)ea;
-        }
-
-        /* Now invert to get actual lifetime. */
-        t->Hi->lifetime = 1./t->Hi->lifetime;
+        t->Hi->lifetime = iso_state_lifetime( ipISO, nelem, nHi, 1 );
+
+        if( ipISO == ipHE_LIKE )
+        {
+                /* iso_state_lifetime is not spin specific, must exclude helike triplet here. */
+                t->Hi->lifetime /= 3.;
+                /* there is also necessary to correct the helike lifetimes */
+                t->Hi->lifetime *= 1.1722 * pow( (double)nelem, 0.1 );
+        }
+
+        /* would probably need a new lifetime algorithm for any other iso sequences. */
+        ASSERT( ipISO <= ipHE_LIKE );
 
         t->Emis->dampXvel = (realnum)( 1.f / t->Hi->lifetime / PI4 / t->EnergyWN );
@@ -1324 +1323 @@
         }
         return;
+}
+
+/* calculate radiative lifetime of an individual iso state */
+STATIC double iso_state_lifetime( long ipISO, long nelem, long n, long l )
+{
+        /* >>refer hydro lifetimes      Horbatsch, M. W., Horbatsch, M. and Hessels, E. A. 2005, JPhysB, 38, 1765 */
+
+        double epsilon, tau, t0, eps2;
+        /* mass of electron */
+        double m = ELECTRON_MASS;
+        /* nuclear mass */
+        double M = (double)dense.AtomicWeight[nelem] * ATOMIC_MASS_UNIT;
+        double mu = (m*M)/(M+m);
+        long z = 1;
+        long Z = nelem + 1 - ipISO;
+        
+        DEBUG_ENTRY( "iso_state_lifetime()" );
+
+        eps2 = 1. - ( l*l + l + 8./47. - (l+1.)/69./n ) / POW2( (double)n );
+
+        epsilon = sqrt( eps2 );
+
+        t0 = 3. * H_BAR * pow( (double)n, 5.) / 
+                ( 2. * POW4( (double)( z * Z ) ) * pow( FINE_STRUCTURE, 5. ) * mu * POW2( SPEEDLIGHT ) ) *
+                POW2( (m + M)/(Z*m + z*M) );
+
+        tau = t0 * ( 1. - eps2 ) / 
+                ( 1. + 19./88.*( (1./eps2 - 1.) * log( 1. - eps2 ) + 1. - 
+                0.5 * eps2 - 0.025 * eps2 * eps2 ) );
+
+        ASSERT( tau > 0. );
+
+        return tau;
 }
 

trunk/source/optimize_func.cpp

@@ -411 +411 @@
                         }
 
-                        fprintf( ioQQQ, " Optimal command: %s\n", 
+                        fprintf( ioQQQ, " Varied command: %s\n", 
                           input.chCardSav[np] );
                 }

trunk/source/parse.h

@@ -117 +117 @@
 */
 void ParseAtomCO(char *chCard );
+
+/**ParseGrid parse the grid command 
+\param *chCard
+*/
+void ParseGrid(char *chCard);
 
 /**ParseInit parse the init command */

trunk/source/parse_commands.cpp

@@ -937 +937 @@
                 else if( strcmp(chKey4,"GRID") == 0 )
                 {
-                        /* option to optimize the model by varying certain parameters
+                        /* option to run grid of models by varying certain parameters
                          * in reads2 */
-                        ParseOptimize(chCard);
+                        ParseGrid(chCard);
                 }
 

trunk/source/parse_optimize.cpp

@@ -246 +246 @@
         }
 
-        else if( nMatch("GRID",chCard) || nMatch("XSPE",chCard) )
-        {
-                /* RP fake optimizer to run a grid of calculations, also accepts
-                 * keyword XSPEC */
-                strcpy( optimize.chOptRtn, "XSPE" );
-                grid.lgGrid = true;
-
-                /* 06 aug 22, change to accept three parameters: lower and upper limit and number of points. */
-                /* scan off range for the previously selected variable */
-                if( optimize.nparm > 0 )
-                {
-                        realnum highvalue;
-                        long numPoints=0;
-
-                        ASSERT( optimize.nparm < LIMPAR );
-
-                        i = 5;
-                        optimize.varang[optimize.nparm-1][0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
-                        optimize.varang[optimize.nparm-1][1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
-                        grid.paramIncrements[optimize.nparm-1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
-
-                        if( grid.paramIncrements[optimize.nparm-1] <= 0. )
-                        {
-                                fprintf( ioQQQ," The increment (third parameter) must be a positive number.\n" );
-                                fprintf( ioQQQ," Sorry.\n" );
-                                puts( "[Stop in ParseOptimize]" );
-                                cdEXIT( EXIT_FAILURE );
-                        }
-
-                        if( grid.paramIncrements[optimize.nparm-1] > ( optimize.varang[optimize.nparm-1][1] - optimize.varang[optimize.nparm-1][0] ) )
-                        {
-                                fprintf( ioQQQ," The increment (third parameter) must not be greater than the difference between the limits (first and second parameters.\n" );
-                                fprintf( ioQQQ," Sorry.\n" );
-                                puts( "[Stop in ParseOptimize]" );
-                                cdEXIT( EXIT_FAILURE );
-                        }
-
-                        if( lgEOL )
-                        {
-                                fprintf( ioQQQ," This command has changed since the definition given in Porter et al. 2006, PASP, 118, 920.\n" );
-                                fprintf( ioQQQ," The grid command now requires three parameters: lower limit, upper limit, and increment.\n" );
-                                fprintf( ioQQQ," The keywords RANGE and STEPS are no longer necessary.\n" );
-                                fprintf( ioQQQ," Sorry.\n" );
-                                puts( "[Stop in ParseOptimize]" );
-                                cdEXIT( EXIT_FAILURE );
-                        }
-                        else
-                        {
-                                ++optimize.nRangeSet;
-                        }
-
-                        /* Swap if second range parameter is less than the first. */
-                        if( optimize.varang[optimize.nparm-1][1] < optimize.varang[optimize.nparm-1][0] )
-                        {
-                                realnum temp = optimize.varang[optimize.nparm-1][0];
-                                optimize.varang[optimize.nparm-1][0] = optimize.varang[optimize.nparm-1][1];
-                                optimize.varang[optimize.nparm-1][1] = temp;
-                        }
-
-                        ASSERT( optimize.varang[optimize.nparm-1][1] - optimize.varang[optimize.nparm-1][0] > 0. );
-
-                        {
-                                realnum integer, tempNumPoints = (optimize.varang[optimize.nparm-1][1] - 
-                                        optimize.varang[optimize.nparm-1][0])/grid.paramIncrements[optimize.nparm-1];
-                                realnum fraction = modf( tempNumPoints, &integer );
-
-                                if( fraction > 0.99f || fraction < 0.01f )
-                                {
-                                        optimize.varang[optimize.nparm-1][1] += 0.01f * grid.paramIncrements[optimize.nparm-1];
-                                        fprintf( ioQQQ,"\n NOTE The range is very nearly (or exactly) an integer multiple of the increment.\n" );
-                                        fprintf( ioQQQ," NOTE The upper limit has been increased by a small fraction of the increment to ensure that\n" );
-                                        fprintf( ioQQQ," NOTE the specified upper limit is included in the grid.\n\n" );
-                                }
-                        }
-
-                        numPoints = 0;
-                        highvalue = optimize.varang[optimize.nparm-1][0];
-                        /* find the number of parameter values by incrementing from the lower limit up to the upper limit. */
-                        while( highvalue <= optimize.varang[optimize.nparm-1][1] )
-                        {
-                                ++numPoints;
-                                highvalue += grid.paramIncrements[optimize.nparm-1];
-                        }
-
-                        grid.numParamValues[optimize.nparm-1] = numPoints;
-
-                        if( grid.numParamValues[optimize.nparm-1] < 2 )
-                        {
-                                fprintf( ioQQQ, " NOTE must have at least two grid points\n" );
-                        }
-                        else if( grid.numParamValues[optimize.nparm-1] > 20 )
-                        {
-                                fprintf( ioQQQ, " NOTE There are %li grid points.  Are you "
-                                        "sure you want that many?\n" ,
-                                        grid.numParamValues[optimize.nparm-1]);
-                        }
-
-                        grid.numParamValues[optimize.nparm-1] = MAX2( 2, grid.numParamValues[optimize.nparm-1] );
-                }
-        }
-
         else
         {

trunk/source/punch_fits.cpp

@@ -184 +184 @@
 
         bytesAdded = 0;
+
+        fixit(); /* bitpix is wrong when realnum is double? */
 
         bytesAdded += addKeyword_txt( "SIMPLE"  , "T",                                  "file does conform to FITS standard", 1 );

trunk/tsuite/auto/dynamics_wind.in

@@ -37 +37 @@
 //
 // >>chng 03 dec 08, from 12 to 15, chng ots, zones
-assert itrzn < 15
+// >>chng 08 jan 23, from 15 to 21, unknown cause 
+assert itrzn < 21
 // 
 c dynamics_wind.in