Changeset 2079 for branches/newmole/source

Timestamp:

05/17/08 09:51:43 (8 months ago)

Author:

rjrw

Message:

Merged from trunk r2022:2078

Location:

branches/newmole/source

Files:

• atmdat_chianti.cpp (modified) (4 diffs)
• atmdat_lamda.cpp (modified) (2 diffs)
• cddrive.h (modified) (3 diffs)
• cdspec.cpp (modified) (3 diffs)
• cont_createmesh.cpp (modified) (1 diff)
• cpu.cpp (modified) (5 diffs)
• cpu.h (modified) (2 diffs)
• date.h (modified) (1 diff)
• dense_fabden.cpp (modified) (1 diff)
• grains.cpp (modified) (1 diff)
• grid.h (modified) (1 diff)
• grid_xspec.cpp (modified) (8 diffs)
• init_coreload.cpp (modified) (1 diff)
• init_defaults_preparse.cpp (modified) (2 diffs)
• ion_trim.cpp (modified) (2 diffs)
• iso_collide.cpp (modified) (3 diffs)
• iso_error.cpp (modified) (1 diff)
• iso_level.cpp (modified) (1 diff)
• iso_radiative_recomb.cpp (modified) (1 diff)
• nemala2.cpp (modified) (3 diffs)
• parse_caseb.cpp (modified) (2 diffs)
• parse_commands.cpp (modified) (2 diffs)
• parse_punch.cpp (modified) (3 diffs)
• punch_do.cpp (modified) (1 diff)
• punch_fits.cpp (modified) (4 diffs)
• punch_line.cpp (modified) (1 diff)
• radius_increment.cpp (modified) (1 diff)
• rt_diffuse.cpp (modified) (1 diff)
• rt_tau_init.cpp (modified) (2 diffs)
• sanity_check.cpp (modified) (1 diff)
• service.cpp (modified) (5 diffs)
• stars.cpp (modified) (15 diffs)
• version.h (modified) (2 diffs)
• zero.cpp (modified) (2 diffs)

branches/newmole/source/atmdat_chianti.cpp

@@ -18 +18 @@
         /* type is set to 0 for non chianti and 1 for chianti*/
         long int i,j,nMolLevs,intCollTran,intcollindex,ipLo,ipHi;
-        FILE *atmolLevDATA , *atmolTraDATA=NULL,*atmolEleColDATA=NULL;
+        FILE *atmolLevDATA , *atmolTraDATA=NULL,*atmolEleColDATA=NULL,*atmolProColDATA=NULL;
         realnum  fstatwt,fenergyK,fenergyWN,fWLAng,fenergy,feinsteina;
         double fScalingParam,fEnergyDiff,fGF,*xs,*spl,*spl2;
@@ -40 +40 @@
         long int *intNewIndex=NULL,*intOldIndex=NULL;
         int interror;
-        bool lgEneLevOrder;
+        bool lgProtonData=false,lgEneLevOrder;
 
         /*This is to identify where the fields start */
@@ -99 +99 @@
         if( trace.lgTrace )
                 fprintf( ioQQQ," moldat_readin opening %s:",chProColFilename);
+        /*We will set a flag here to indicate if the proton collision strengths are available */
+        if( ( atmolProColDATA = fopen( chProColFilename , "r" ) ) != NULL )
+        {
+                lgProtonData = true;
+        }
+        else
+        {
+                lgProtonData = false;
+        }
 
         nMolLevs = 0;
@@ -146 +155 @@
         for( intcollindex = 0; intcollindex<NUM_COLLIDERS; intcollindex++ )
         {
-                CollRatesArray[intNS][intcollindex] = (double**)MALLOC((nMolLevs)*sizeof(double*));
+                CollRatesArray[intNS][intcollindex] = NULL;
+        }
+        /*Allocating space just for the electron*/
+        CollRatesArray[intNS][0] = (double**)MALLOC((nMolLevs)*sizeof(double*));
+        for( ipHi = 0; ipHi<nMolLevs; ipHi++ )
+        {
+                CollRatesArray[intNS][0][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double));
+                for( ipLo = 0; ipLo<nMolLevs; ipLo++)
+                {
+                        CollRatesArray[intNS][0][ipHi][ipLo] = 0.0;
+                }
+        }
+        /*Allocating space for the proton*/
+        if(lgProtonData)
+        {
+                CollRatesArray[intNS][1] = (double**)MALLOC((nMolLevs)*sizeof(double*));
                 for( ipHi = 0; ipHi<nMolLevs; ipHi++ )
                 {
-                        CollRatesArray[intNS][intcollindex][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double));
+                        CollRatesArray[intNS][1][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double));
                         for( ipLo = 0; ipLo<nMolLevs; ipLo++)
                         {
-                                CollRatesArray[intNS][intcollindex][ipHi][ipLo] = 0.0;
-                        }
-                }
-        }
-
+                                CollRatesArray[intNS][1][ipHi][ipLo] = 0.0;
+                        }
+                }
+        }
         /*Rewind the energy levels files*/
         if( fseek( atmolLevDATA , 0 , SEEK_SET ) != 0 )

branches/newmole/source/atmdat_lamda.cpp

@@ -113 +113 @@
         for( intcollindex = 0; intcollindex <NUM_COLLIDERS; intcollindex++ )
         {
-                CollRatesArray[intNS][intcollindex] = (double**)MALLOC((unsigned long)(ngMolLevs) * sizeof(double*));
-                for( ipHi = 0; ipHi<ngMolLevs; ipHi++ )
-                {
-                        CollRatesArray[intNS][intcollindex][ipHi] = (double*)MALLOC((unsigned long)(ngMolLevs) * sizeof(double));
-                        for( ipLo = 0; ipLo<ngMolLevs; ipLo++ )
-                        {
-                                CollRatesArray[intNS][intcollindex][ipHi][ipLo] = 0.0;
-                        }
-                }
-        }
+                CollRatesArray[intNS][intcollindex] = NULL;
+        }
+
 
         if(read_whole_line( chLine , (int)sizeof(chLine) , atmolLevDATA ) == NULL )
@@ -342 +335 @@
                         TotalInsanity();
                 }
+                /*This is where we allocate memory if the collider exists*/
+                /*Needed to take care of the he collisions*/
+                CollRatesArray[intNS][intCollIndex] = (double**)MALLOC((unsigned long)(ngMolLevs) * sizeof(double*));
+                for( ipHi = 0; ipHi<ngMolLevs; ipHi++ )
+                {
+                        CollRatesArray[intNS][intCollIndex][ipHi] = (double*)MALLOC((unsigned long)(ngMolLevs) * sizeof(double));
+                        for( ipLo = 0; ipLo<ngMolLevs; ipLo++ )
+                        {
+                                CollRatesArray[intNS][intCollIndex][ipHi][ipLo] = 0.0;
+                        }
+                }
                 if(read_whole_line( chLine , (int)sizeof(chLine)  , atmolLevDATA ) == NULL )
                 {

branches/newmole/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[] );
 

branches/newmole/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;
 }

branches/newmole/source/cont_createmesh.cpp

@@ -457 +457 @@
          * 
          * TeLowestFineOpacity of 1e4 K is temperature were line width is 
-         * evaluationed.  Tests were done using the stop temperature in its place
+         * evaluated.  Tests were done using the stop temperature in its place
          * Te below 1e4 K made fine opacity grid huge 
          * do not let temp get higher than 1e4 either - code run with stop temp 10 set
          * stop temp of 1e10K and assert thrown at line 204 of cont_createpointers.c 
          * simply use 1e4 K as a characteristic temperature */
+        /** \todo       1       set temp of 1e4K will be too coarse a line for PDRs where
+         * H2 line overlap is very important */
         double TeLowestFineOpacity = 1e4;
         rfield.fine_opac_velocity_width = 
                 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*TeLowestFineOpacity/
-                /* TurbVel/5 will divide line into 5 cells when turbulence-dominated */
-                dense.AtomicWeight[rfield.fine_opac_nelem] + POW2(DoppVel.TurbVel/5.)) / 
+                dense.AtomicWeight[rfield.fine_opac_nelem] ) / 
                 /* we want fine_opac_nresolv continuum elements across this line
                  * default is 1, changed with SET FINE CONTINUUM command */

branches/newmole/source/cpu.cpp

@@ -359 +359 @@
 
         /* >>chng 01 oct 09, disable gradual underflow on ultrasparc whith g++
-         * (only needed for versions prior to 3.1, see Note 1).
+         * (only needed for versions < 3.1 or >= 4.3.0, see Note 1).
          *
          * compile this module with:
-         *     g++ [ other options... ] -I<include-dir> -DHAVE_SUNMATH -c setfpenv.c
+         *     g++ [ other options... ] -I<include-dir> -DHAVE_SUNMATH -c cpu.cpp
          * link the program with:
          *     g++ -L<library-dir> -o cloudy.exe *.o -lsunmath
@@ -379 +379 @@
          * g++ [ other options... ] -funsafe-math-optimizations -c *.c
          * g++ [ other options... ] -funsafe-math-optimizations -o cloudy.exe *.o
+         *
+         * Starting with g++ 4.3.0 the -funsafe-math-optimizations option can no longer be
+         * used as it implicitly enables -fno-trapping-math, which is unsafe for Cloudy
+         * because we do trap floating point exceptions.
          *
          * Note 2: Don't use nonstandard_arithmetic() with CC (the SunWorks/Forte compiler);
@@ -439 +443 @@
         string FileName( fname );
         vector<string>::size_type begin, end;
-        bool lgAbort = ( scheme == AS_DATA_ONLY || scheme == AS_DATA_LOCAL ||
+        bool lgAbort = ( scheme == AS_DATA_ONLY || scheme == AS_DATA_OPTIONAL || scheme == AS_DATA_LOCAL ||
                          scheme == AS_LOCAL_DATA || scheme == AS_LOCAL_ONLY );
 
@@ -446 +450 @@
         case AS_DATA_ONLY:
         case AS_DATA_ONLY_TRY:
+        case AS_DATA_OPTIONAL:
                 begin = 1;
                 end = cpu.chSearchPath.size()-1;
@@ -480 +485 @@
         if( handle == NULL && lgAbort )
         {
-                // failed on very first open -> probably path is not correct
-                fprintf( ioQQQ, "\nPROBLEM DISASTER I could not open the data file %s\n\n", fname );
-                fprintf( ioQQQ, "Although there may be other reasons you have received this error,\n");
-                fprintf( ioQQQ, "the most likely are that the path has not been properly set\n");
-                fprintf( ioQQQ, "or that the path points to an old version of the data.\n\n");
-                fprintf( ioQQQ, "Please have a look at the file path.h in the source directory\n");
-                fprintf( ioQQQ, "to check how the variable CLOUDY_DATA_PATH is set - \n");
-                fprintf( ioQQQ, "it should give the location of the data files I need.\n");
-                fprintf( ioQQQ, "These are the files in the data download from the web site.\n\n");
-                fprintf( ioQQQ, "Recompile the code with the correct data path set in path.h\n");
-                fprintf( ioQQQ, "or use the shell command \"set CLOUDY_DATA_PATH=path\" to set the\n");
-                fprintf( ioQQQ, "path from a command prompt.\n\n");
-                cpu.printDataPath();
-
-                fprintf( ioQQQ, " These are all the paths I tried:\n" );
-                for( vector<string>::size_type i=begin; i < end; ++i )
-                {
-                        if( cpu.chSearchPath[i].length() > 0 )
-                                fprintf( ioQQQ, "   ==%s==\n", cpu.chSearchPath[i].c_str() );
-                        else
-                                fprintf( ioQQQ, "   ==<local directory>==\n" );
+                if( scheme == AS_DATA_OPTIONAL )
+                        // presence is optional -> make warning less scary...
+                        fprintf( ioQQQ, "\nI could not open the data file %s\n\n", fname );
+                else
+                        fprintf( ioQQQ, "\nPROBLEM DISASTER I could not open the data file %s\n\n", fname );
+                if( cpu.nFileDone == 0 || scheme == AS_DATA_ONLY )
+                {
+                        // failed on very first open -> most likely path is not correct
+                        // failed on AS_DATA_ONLY -> CLOUDY_DATA_PATH may point to obsolete data dir
+                        fprintf( ioQQQ, "Although there may be other reasons you have received this error,\n");
+                        fprintf( ioQQQ, "the most likely are that the path has not been properly set\n");
+                        fprintf( ioQQQ, "or that the path points to an old version of the data.\n\n");
+                        fprintf( ioQQQ, "Please have a look at the file path.h in the source directory\n");
+                        fprintf( ioQQQ, "to check how the variable CLOUDY_DATA_PATH is set - \n");
+                        fprintf( ioQQQ, "it should give the location of the data files I need.\n");
+                        fprintf( ioQQQ, "These are the files in the data download from the web site.\n\n");
+                        fprintf( ioQQQ, "Recompile the code with the correct data path set in path.h\n");
+                        fprintf( ioQQQ, "or use the shell command \"export CLOUDY_DATA_PATH=path\" to set the\n");
+                        fprintf( ioQQQ, "path from a bash command prompt.\n\n");
+                        cpu.printDataPath();
+                }
+                else
+                {
+                        // failed on search including local directory -> most likely the file name
+                        // was mistyped on a compile command, or Cloudy is run in the wrong directory
+                        // if scheme == AS_DATA_OPTIONAL, this most likely is a stellar grid that is not installed.