Changeset 2036

Timestamp:

05/10/08 11:14:00 (6 months ago)

Author:

rporter

Message:

Add option to limit energy range in XSPEC option. Add line emission to diffuse atable option. In two comments, punch_line.c -> punch_line.cpp. Modify ism_grid to include the new range option on XSPEC output.

Location:

trunk

Files:

• source/cddrive.h (modified) (3 diffs)
• source/cdspec.cpp (modified) (3 diffs)
• source/dense_fabden.cpp (modified) (1 diff)
• source/grid.h (modified) (1 diff)
• source/grid_xspec.cpp (modified) (7 diffs)
• source/init_coreload.cpp (modified) (1 diff)
• source/parse_punch.cpp (modified) (2 diffs)
• source/punch_fits.cpp (modified) (4 diffs)
• source/punch_line.cpp (modified) (1 diff)
• tsuite/auto/ism_grid.in (modified) (1 diff)

trunk/source/cddrive.h

@@ -536 +536 @@
          3                      the reflected continuum, same units
         
-         4                      diffuse continuous emission outward direction
-         5                      diffuse continuous emission, reflected
+         4                      diffuse emission outward direction
+         5                      diffuse emission, reflected
         
          6                      collisional+recombination lines, outward
@@ -547 +547 @@
          
         \param nEnergy the number of cells + 1
+        \param ipLoEnergy
+        \param ipHiEnergy
         \param ReturnedSpectrum[] the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts
  */ 
@@ -552 +554 @@
         int Option ,
     long int nEnergy ,
+        long int ipLoEnergy,
+        long int ipHiEnergy,
     realnum ReturnedSpectrum[] );
 

trunk/source/cdspec.cpp

@@ -318 +318 @@
          * 3    the reflected continuum, same units
          *
-         * 4    diffuse continuous emission outward direction
-         * 5    diffuse continuous emission, reflected
-         *
-         * 6    collisional+recombination lines, outward
-         * 7    collisional+recombination lines, reflected
+         * 4    diffuse emission, lines + continuum, outward
+         * 5    diffuse emission, lines + continuum, reflected
+         *
+         * 6    diffuse continuous emission outward direction
+         * 7    diffuse continuous emission, reflected
          *
          * 8    total transmitted, lines and continuum
@@ -333 +333 @@
         int nOption ,
 
-        /* the number of cells + 1*/
-        long int nEnergy ,
+        /* the number of cells */
+        long int nEnergy,
+
+        /* the index of the lowest and highest energy bounds to use. */
+        long ipLoEnergy,
+        long ipHiEnergy,
 
         /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
@@ -340 +344 @@
 
 {
-        /* this pointer will bet set to one of the cloudy continuum arrays */
-        realnum *cont , 
-                refac;
-        long int ncell , j;
-
-        /* flag saying whether we will need to free cont at the end */
-        bool lgFREE;
+        realnum refac;
 
         DEBUG_ENTRY( "cdSPEC2()" );
 
+        ASSERT( ipLoEnergy >= 0 );
+        ASSERT( ipLoEnergy < ipHiEnergy );
+        ASSERT( ipHiEnergy < rfield.nupper );
+        ASSERT( nEnergy == (ipHiEnergy-ipLoEnergy+1) );
+        ASSERT( nEnergy >= 2 );
+
         ASSERT( nOption <= NUM_OUTPUT_TYPES );
 
-        if( nOption == 1 )
-        {
-                /* this is the incident continuum, col 2 of punch continuum command */
-                cont = rfield.flux_total_incident;
-                lgFREE = false;
-        }
-        else if( nOption == 2 )
-        {
-                /* the attenuated transmitted continuum, no diffuse emission,
-                 * col 3 of punch continuum command */
-                cont = rfield.flux;
-                lgFREE = false;
-        }
-        else if( nOption == 3 )
-        {
-                /* reflected incident continuum, col 6 of punch continuum command */
-                lgFREE = false;
-                cont = rfield.ConRefIncid;
-        }
-        else if( nOption == 4 )
-        {
-                /* diffuse continuous emission outward direction  */
-                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
-                /* need to free the vector once done */
-                lgFREE = true;
-                refac = (realnum)radius.r1r0sq*geometry.covgeo;
-                for( j=0; j<rfield.nflux; ++j)
-                {
-                        cont[j] = rfield.ConEmitOut[j]*refac;
-                }
-        }
-        else if( nOption == 5 )
-        {
-                /* reflected diffuse continuous emission */
-                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
-                /* need to free the vector once done */
-                lgFREE = true;
-                for( j=0; j<rfield.nflux; ++j)
-                {
-                        cont[j] = rfield.ConEmitReflec[j];
-                }
-        }
-        else if( nOption == 6 )
-        {
-                /* all outward lines,   */
-                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
-                /* need to free the vector once done */
-                lgFREE = true;
-                /* correct back to inner radius */
-                refac = (realnum)radius.r1r0sq*geometry.covgeo;
-                for( j=0; j<rfield.nflux; ++j)
-                {
-                        cont[j] = rfield.outlin[j]*refac;
-                }
-        }
-        else if( nOption == 7 )
-        {
-                /* all reflected lines */
-                if( geometry.lgSphere )
-                {
-                        refac = 0.;
-                }
-                else
-                {
-                        refac = 1.;
-                }
-
-                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
-                /* need to free the vector once done */
-                lgFREE = true;
-                for( j=0; j<rfield.nflux; ++j)
-                {
-                        /* units of lines here are to cancel out with tricks applied to continuum cells
-                         * when finally extracted below */
-                        cont[j] = rfield.reflin[j]*refac;
-                }
-        }
-        else if( nOption == 8 )
-        {
-                /* total transmitted continuum */
-                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
-                /* need to free the vector once done */
-                lgFREE = true;
-                /* correct back to inner radius */
-                refac = (realnum)radius.r1r0sq*geometry.covgeo;
-                for( j=0; j<rfield.nflux; ++j)
-                {
-                        /* units of lines here are to cancel out with tricks applied to continuum cells
-                         * when finally extracted below */
-                        cont[j] = (rfield.ConEmitOut[j]+ rfield.outlin[j])*refac
-                                          + rfield.flux[j]*(realnum)radius.r1r0sq;
-                }
-        }
-        else if( nOption == 9 )
-        {
-                /* total reflected continuum */
-                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
-                /* need to free the vector once done */
-                lgFREE = true;
-                for( j=0; j<rfield.nflux; ++j)
-                {
-                        /* units of lines here are to cancel out with tricks applied to continuum cells
-                         * when finally extracted below */
-                        cont[j] = ((rfield.ConRefIncid[j]+rfield.ConEmitReflec[j]) +
+        for( long i = 0; i < nEnergy; i++ )
+        {
+                long j = ipLoEnergy + i;
+
+                if( j >= rfield.nflux )
+                {
+                        ReturnedSpectrum[i] = SMALLFLOAT;
+                        continue;
+                }
+
+                if( nOption == 1 )
+                {
+                        /* this is the incident continuum, col 2 of punch continuum command */
+                        ReturnedSpectrum[i] = rfield.flux_total_incident[j];
+                }
+                else if( nOption == 2 )
+                {
+                        /* the attenuated transmitted continuum, no diffuse emission,
+                         * col 3 of punch continuum command */
+                        ReturnedSpectrum[i] = rfield.flux[j];
+                }
+                else if( nOption == 3 )
+                {
+                        /* reflected incident continuum, col 6 of punch continuum command */
+                        ReturnedSpectrum[i] = rfield.ConRefIncid[j];
+                }
+                else if( nOption == 4 )
+                {
+                        /* correct back to inner radius */
+                        refac = (realnum)radius.r1r0sq*geometry.covgeo;
+                        ReturnedSpectrum[i] = (rfield.outlin[j]+rfield.ConEmitOut[j])*refac;
+                }
+                else if( nOption == 5 )
+                {
+                        /* all reflected diffuse emission */
+                        if( geometry.lgSphere )
+                        {
+                                refac = 0.;
+                        }
+                        else
+                        {
+                                refac = 1.;
+                        }
+
+                        ReturnedSpectrum[i] = (rfield.reflin[j]+rfield.ConEmitReflec[j])*refac;
+                }
+                else if( nOption == 6 )
+                {
+                        /* correct back to inner radius */
+                        refac = (realnum)radius.r1r0sq*geometry.covgeo;
+                        ReturnedSpectrum[i] = rfield.outlin[j]*refac;
+                }
+                else if( nOption == 7 )
+                {
+                        /* all reflected line emission */
+                        if( geometry.lgSphere )
+                        {
+                                refac = 0.;
+                        }
+                        else
+                        {
+                                refac = 1.;
+                        }
+
+                        ReturnedSpectrum[i] = rfield.reflin[j]*refac;
+                }
+                else if( nOption == 8 )
+                {
+                        /* total transmitted continuum */
+                        /* correct back to inner radius */
+                        refac = (realnum)radius.r1r0sq*geometry.covgeo;
+                        ReturnedSpectrum[i] = (rfield.ConEmitOut[j]+ rfield.outlin[j])*refac
+                                + rfield.flux[j]*(realnum)radius.r1r0sq;
+                }
+                else if( nOption == 9 )
+                {
+                        /* total reflected continuum */
+                        ReturnedSpectrum[i] = ((rfield.ConRefIncid[j]+rfield.ConEmitReflec[j]) +
                                 rfield.reflin[j]);
                 }
-        }
-        else if( nOption == 10 )
-        {
-                /* this is exp(-tau) */
-                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
-                /* need to free the vector once done */
-                lgFREE = true;
-                for( j=0; j<nEnergy; ++j)
-                {
+                else if( nOption == 10 )
+                {
+                        /* this is exp(-tau) */
                         /* This is the TOTAL attenuation in both the continuum and lines.  
                          * Jon Miller discovered that the line attenuation was missing in 07.02 */
-                        cont[j] = opac.ExpmTau[j]*rfield.trans_coef_total[j];
-                }
-        }
-        else
-        {
-                fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
-                cdEXIT(EXIT_FAILURE);
-        }
-
-        /* now generate the continua */
-        for( ncell = 0; ncell < nEnergy; ++ncell )
-        {
-        if( ncell >= rfield.nflux )
-                {
-                        ReturnedSpectrum[ncell] = SMALLFLOAT;
+                        ReturnedSpectrum[i] = opac.ExpmTau[j]*rfield.trans_coef_total[j];
                 }
                 else
                 {
-                        ReturnedSpectrum[ncell] = cont[ncell];
-                }
-                ASSERT( ReturnedSpectrum[ncell] >=0.f );
-        }
-
-        /* need to free the vector once done if this flag was set */
-        if( lgFREE )
-        {
-                free(cont);
-        }
+                        fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
+                        cdEXIT(EXIT_FAILURE);
+                }
+
+                ASSERT( ReturnedSpectrum[i] >=0.f );
+        }
+
         return;
 }

trunk/source/dense_fabden.cpp

@@ -10 +10 @@
   double depth)
 {
-        double a , b;
-        a = depth;
-        b = radius;
-
-        fprintf(ioQQQ,"PROBLEM DISASTER The default version of dense_fabden is "
-                "still in dense_fabden.cpp - edit and replace this file to use your version.\n");
-        /* return val is example - must return density */
-        return(a*b);
-
-        /* these are two examples of using this routine */
-#       if 0
-        double fabden_v;
-        double z , z0 , density , rad_midplane;
-        DEBUG_ENTRY( "dense_fabden()" );
-
-        /* this routine is called when the DLAW command is given,
-         * and must return the hydrogen density at the current radius or depth
-         * RADIUS is the radius, the distance from the center of symmetry.
-         * DEPTH is the depth into the cloud, from the illuminated face
-         * both are in cm
-         *
-         * radius, depth, and the array DensityLaw are double precision, although
-         * dense_fabden itself is not
-         *
-         * this is one way to generate a density
-         * dense_fabden = radius * depth
-         *
-         * this is how to use the parameters on the dlaw command
-         * dense_fabden = DensityLaw(1) + radius * DensityLaw(2)
-         *
-         * following must be removed if this sub is to be used */
-        /*fprintf( ioQQQ, " the old version of dense_fabden must be deleted, and a new one put in place.  Sorry.\n" );*/
-
-        if( depth > radius )
-                TotalInsanity();
-
-        rad_midplane = radius * cos( dense.DensityLaw[0] );
-        z = radius * sin( dense.DensityLaw[0] );
-        z0 = 0.047 * pow( rad_midplane/AU , 1.25 ) * AU;
-        /* density in gm / cm3 */
-        density = 1.4e-9 * pow( rad_midplane/AU , -11./4. ) * sexp( POW2(z/z0) );
-        /* convert density to cm-3 */
-        fabden_v = density/(ATOMIC_MASS_UNIT * 1.2);
-
-        fprintf(ioQQQ, "bug dense_fabden zone %.2f radius %e density %e \n",
-                fnzone , radius , fabden_v );
-        return( fabden_v );
-#       endif
-
-#if 0
-
-        /* update stromgren radius to most recent value */
+        double fabden_v = pow(10.,dense.DensityLaw[0]);
         if( rfield.lgUSphON )
-                dense.DensityLaw[2] = rfield.rstrom/PARSEC;
-
-        /* this is the density law used in the Wen & O'Dell, 1995, ApJ 438, 784 paper */
-        powexp = MIN2((radius/parsec-dense.DensityLaw[2])/dense.DensityLaw[1],dense.DensityLaw[3]);
-        fabden_v = pow(10.,dense.DensityLaw[0])*exp(powexp);
-
-        fabden_v = dense.DensityLaw[0];
-        /* just here to stop compilers from flagging unused vars */
-        temp = radius + depth; 
-        if( fabden_v == 0. )
-        {
-                cdEXIT(EXIT_FAILURE);
-        }
-        else
-        {
-                /* when this routine is used the following branch is the correct exit */
-                /* following should return correct value of the density at this position, 
-                 * temp is there to trick lint */
-                return( fabden_v*temp );
-        }
-#endif
+                fabden_v *= pow(radius/rfield.rstrom,dense.DensityLaw[1]);
+        return fabden_v; 
 }

trunk/source/grid.h

@@ -57 +57 @@
         bool lgOutputTypeOn[NUM_OUTPUT_TYPES];
 
+        long ipLoEnergy, ipHiEnergy;
+        realnum LoEnergy_keV, HiEnergy_keV;
+
 } grid;
 

trunk/source/grid_xspec.cpp

@@ -9 +9 @@
 #include "optimize.h"
 #include "cddrive.h"
+#include "continuum.h"
 #include "rfield.h"
 #include "grid.h"
+#include "ipoint.h"
 #include "called.h"
+#include "physconst.h"
 #include "prt.h"
 
@@ -17 +20 @@
 void gridXspec(realnum xc[], long int nInterpVars)
 {
-        long int i, j, k;
-        double averageChi2;
+        long int i;
 
         DEBUG_ENTRY( "gridXspec()" );
@@ -79 +81 @@
                 grid.paramRange[i][5] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f;
 
-                for( j=0; j<grid.numParamValues[i]; j++ )
+                for( long j=0; j<grid.numParamValues[i]; j++ )
                 {
                         grid.paramData[i][j] = xc[i]+grid.paramIncrements[i]*j;
@@ -90 +92 @@
         }
 
-        /* >>chng 06 aug 23, this logic now allows non-square parameter spaces.  */
         for( i=0; i< grid.totNumModels; i++ )
         {
+                long j;
+                double averageChi2;
                 realnum variableVector[LIMPAR];
 
@@ -105 +108 @@
                          * second "volume" is product of grid.numParamValues[2]*grid.numParamValues[3]*....grid.numParamValues[n]
                          * last "volume" is unity.  */
-                        for( k=j+1; k<nInterpVars; k++ )
+                        for( long k=j+1; k<nInterpVars; k++ )
                         {
                                 volumeOtherDimensions *= grid.numParamValues[k];
@@ -210 +213 @@
         {
                 long i1, i2;
-                grid.numEnergies = rfield.nupper-2;
+
+                if( grid.LoEnergy_keV == 0. )
+                        grid.ipLoEnergy = 0;
+                else
+                        grid.ipLoEnergy = ipoint( grid.LoEnergy_keV * 1000. / EVRYD );
+
+                if( grid.HiEnergy_keV == 0. || grid.HiEnergy_keV >= continuum.filbnd[continuum.nrange] )
+                        grid.ipHiEnergy = rfield.nflux - 1;
+                else
+                        grid.ipHiEnergy = ipoint( grid.HiEnergy_keV * 1000. / EVRYD );
+
+                grid.numEnergies = grid.ipHiEnergy - grid.ipLoEnergy + 1;
                 grid.Energies = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numEnergies) );
                 grid.Spectra = (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(NUM_OUTPUT_TYPES) );
@@ -241 +255 @@
                  * at this point nOptimiz has already been incremented for first model */
                 if( grid.lgOutputTypeOn[i] )
-                        cdSPEC2( i, grid.numEnergies, grid.Spectra[i][optimize.nOptimiz]);
+                        cdSPEC2( i, grid.numEnergies, grid.ipLoEnergy, grid.ipHiEnergy, grid.Spectra[i][optimize.nOptimiz]);
         }
         return;

trunk/source/init_coreload.cpp

@@ -306 +306 @@
         grid.lgGridDone = false;
         grid.lgStrictRepeat = false;
+        /* these are energy range... if not changed with command, 0. says just use energy limits of mesh */
+        grid.LoEnergy_keV = 0.;
+        grid.HiEnergy_keV = 0.;
+
 
         PunchFilesInit();

trunk/source/parse_punch.cpp

@@ -1619 +1619 @@
                         strcpy( punch.chPunch[punch.npunch], "LINR" );
                         /* punch some details needed for line radiative transfer 
-                         * routine in punch_line.c */
+                         * routine in punch_line.cpp */
                         Punch_Line_RT(punch.ipPnunit[punch.npunch],"READ");
                 }
@@ -2143 +2143 @@
                 {
                         punch.lgPunLstIter[punch.npunch] = true;
+                }
+
+                ipFFmt = 5;
+                /* range option - important since so much data */
+                if( nMatch("RANGE",chCard) ) 
+                {
+                        /* get lower and upper range, must be in keV */
+                        punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
+                        punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
+                        if( lgEOL )
+                        {
+                                fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in keV.\nSorry.\n");
+                                cdEXIT(EXIT_FAILURE);
+                        }
+                        if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] )
+                        {
+                                fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
+                                cdEXIT(EXIT_FAILURE);
+                        }
+
+                        grid.LoEnergy_keV = punch.punarg[punch.npunch][0];
+                        grid.HiEnergy_keV = punch.punarg[punch.npunch][1];
+                }
+                else
+                {
+                        /* these mean full energy range */
+                        punch.punarg[punch.npunch][0] = 0;
+                        punch.punarg[punch.npunch][1] = 0;
                 }
 

trunk/source/punch_fits.cpp

@@ -77 +77 @@
                                                                   long totNumModels, long numEnergies, long nintparm, long naddparm  );
 STATIC void punchFITS_GenericHeader( long numEnergies );
-STATIC void punchFITS_GenericData( long numEnergies );
+STATIC void punchFITS_GenericData( long numEnergies, long ipLoEnergy, long ipHiEnergy );
 STATIC void writeCloudyDetails( void );
 STATIC long addComment( const char *CommentToAdd );
@@ -128 +128 @@
                 punchFITS_PrimaryHeader( false );
                 punchFITS_GenericHeader( rfield.nflux - 1 );
-                punchFITS_GenericData( rfield.nflux -1 );
+                punchFITS_GenericData( rfield.nflux -1, 0, rfield.nflux -2 );
         }
         /* These are specially designed XSPEC outputs. */
@@ -673 +673 @@
 }
 
-STATIC void punchFITS_GenericData( long numEnergies )
+STATIC void punchFITS_GenericData( long numEnergies, long ipLoEnergy, long ipHiEnergy )
 {
         long i;
@@ -683 +683 @@
         TransmittedSpectrum = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(numEnergies) );
 
-        cdSPEC2( 8, numEnergies, TransmittedSpectrum );
+        cdSPEC2( 8, numEnergies, ipLoEnergy, ipHiEnergy, TransmittedSpectrum );
 
         /* Now add the energies data */

trunk/source/punch_line.cpp

@@ -343 +343 @@
                         if(nLine>=LIMLINE )
                         {
-                                fprintf(ioQQQ," PUNCH RT has too many lines - increase LIMLINE in punch_line.c\n");
+                                fprintf(ioQQQ," PUNCH RT has too many lines - increase LIMLINE in punch_line.cpp\n");
                                 cdEXIT(EXIT_FAILURE);
                         }

trunk/tsuite/auto/ism_grid.in

@@ -39 +39 @@
 punch heating "ism_grid.het" 
 punch coolign "ism_grid.col" 
-punch xspec mtable "ism_grid.fit" 
+punch xspec mtable "ism_grid.fit" range 0.1 3 keV
 punch temperature "ism_grid.tem" 
 c