cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "physconst.h" 00005 #include "optimize.h" 00006 #include "continuum.h" 00007 #include "called.h" 00008 #include "rfield.h" 00009 #include "stars.h" 00010 /*lint -e785 too few initializers */ 00011 /*lint -e801 use of go to depreciated */ 00012 00014 static const int NSB99 = 1250; 00016 static const int MNTS = 200; 00017 00019 static const int NRAUCH = 19951; 00021 static const int NMODS_HCA = 66; 00023 static const int NMODS_HNI = 51; 00025 static const int NMODS_PG1159 = 71; 00027 static const int NMODS_HYDR = 100; 00029 static const int NMODS_HELIUM = 81; 00031 static const int NMODS_HpHE = 117; 00032 00033 /* set to 1 to turn on debug print statements in these routines */ 00034 #define DEBUGPRT 0 00035 00036 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; } 00037 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; } 00038 00039 static const bool lgSILENT = false; 00040 static const bool lgVERBOSE = true; 00041 00042 static const bool lgLINEAR = false; 00043 static const bool lgTAKELOG = true; 00044 00045 typedef enum { 00046 IS_UNDEFINED, IS_FIRST, IS_SECOND 00047 } IntStage; 00048 00050 typedef struct 00051 { 00052 double par[MDIM]; 00053 int modid; 00054 char chGrid; 00055 } mpp; 00056 00066 /* this is the structure of the binary atmosphere file (VERSION 20060612): 00067 * 00068 * ============================ 00069 * * int32 VERSION * 00070 * * int32 MDIM * 00071 * * int32 MNAM * 00072 * * int32 ndim * 00073 * * int32 npar * 00074 * * int32 nmods * 00075 * * int32 ngrid * 00076 * * uint32 nOffset * 00077 * * uint32 nBlocksize * 00078 * * char names[MDIM][MNAM+1] * 00079 * * mpp telg[nmods] * 00080 * * realnum anu[ngrid] * 00081 * * realnum mod1[ngrid] * 00082 * * ... * 00083 * * realnum modn[ngrid] * 00084 * ============================ 00085 * 00086 * nOffset == 7*sizeof(int32) + 2*sizeof(uint32) + MDIM*(MNAM+1)*sizeof(char) + nmods*sizeof(mpp) 00087 * nBlocksize == ngrid*size(realnum) */ 00088 00090 typedef struct 00091 { 00093 string name; 00095 access_scheme scheme; 00097 FILE *ioIN; 00100 const char *ident; 00102 const char *command; 00104 IntMode imode; 00106 int32 ndim; 00108 int32 npar; 00110 int32 nmods; 00112 int32 ngrid; 00114 uint32 nOffset; 00116 uint32 nBlocksize; 00119 mpp *telg; /* telg[nmods] */ 00121 double **val; /* val[ndim][nval[n]] */ 00123 long *nval; /* nval[ndim] */ 00130 long *jlo; /* jlo(nval[0],...,nval[ndim-1]) */ 00131 long *jhi; /* jhi(nval[0],...,nval[ndim-1]) */ 00133 char names[MDIM][MNAM+1]; 00135 long *trackLen; /* trackLen[nTracks] */ 00137 long nTracks; 00139 long *jval; 00140 } stellar_grid; 00141 00142 /* internal routines */ 00143 STATIC bool lgCompileAtmosphereCoStar(const char[],const char[],const realnum[],long,process_counter&); 00144 STATIC void InterpolateGridCoStar(const stellar_grid*,const double[],double*,double*); 00145 STATIC void FindHCoStar(const stellar_grid*,long,double,long,realnum*,long*,long*); 00146 STATIC void FindVCoStar(const stellar_grid*,double,realnum*,long[]); 00147 STATIC void CoStarListModels(const stellar_grid*); 00148 STATIC int RauchInitializeSub(const char[],const char[],const mpp[],long,long, 00149 long,const double[],int); 00150 STATIC bool lgCompileAtmosphere(const char[],const char[],const realnum[],long,process_counter&); 00151 STATIC void InitGrid(stellar_grid*,bool); 00152 STATIC bool lgValidBinFile(const char*,process_counter&,access_scheme); 00153 STATIC bool lgValidAsciiFile(const char*,access_scheme); 00154 STATIC void InitGridCoStar(stellar_grid*); 00155 STATIC void CheckVal(const stellar_grid*,double[],long*,long*); 00156 STATIC void InterpolateRectGrid(const stellar_grid*,const double[],double*,double*); 00157 STATIC void FreeGrid(stellar_grid*); 00158 STATIC void InterpolateModel(const stellar_grid*,const double[],double[],const long[], 00159 const long[],long[],long,realnum[],IntStage); 00160 STATIC void InterpolateModelCoStar(const stellar_grid*,const double[],double[],const long[], 00161 const long[],long[],long,long,realnum[]); 00162 STATIC void GetModel(const stellar_grid*,long,realnum[],bool,bool); 00163 STATIC void SetLimits(const stellar_grid*,double,const long[],const long[],const long[], 00164 const realnum[],double*,double*); 00165 STATIC void SetLimitsSub(const stellar_grid*,double,const long[],const long[],long[],long, 00166 double*,double*); 00167 STATIC void InitIndexArrays(stellar_grid*,bool); 00168 STATIC void FillJ(const stellar_grid*,long[],double[],long,bool); 00169 STATIC long JIndex(const stellar_grid*,const long[]); 00170 STATIC void SearchModel(const mpp[],long,const double[],long,long*,long*); 00171 STATIC void FindIndex(const double[],long,double,long*,long*,bool*); 00172 STATIC bool lgFileReadable(const char*, process_counter&,access_scheme); 00173 STATIC void ValidateGrid(const stellar_grid*,double); 00174 STATIC bool lgValidModel(const realnum[],const realnum[],double,double); 00175 STATIC void RebinAtmosphere(long,const realnum[],const realnum[],realnum[],long,const realnum[]); 00176 STATIC realnum RebinSingleCell(realnum,realnum,const realnum[],const realnum[],const realnum[],long); 00177 STATIC long RebinFind(const realnum[],long,realnum); 00178 00179 00180 /* the version number for the ascii/binary atmosphere files */ 00181 static const long int VERSION_ASCII = 20060612L; 00182 /* binary files are incompatible when floats are converted to doubles */ 00183 #ifdef FLT_IS_DBL 00184 static const long int VERSION_BIN = 20071004L; 00185 #else 00186 static const long int VERSION_BIN = 20060612L; 00187 #endif 00188 00190 void AtmospheresAvail( void ) 00191 { 00192 DEBUG_ENTRY( "AtmospheresAvail()" ); 00193 00194 /* This routine makes a list of all the stellar atmosphere grids that are valid, 00195 * giving the parameters for use in the input script as well. It is simply a long 00196 * list of if-statements, so if any grid is added to Cloudy, it should be added in 00197 * this routine as well. 00198 * 00199 * NB NB NB -- test this routine regularly to see if the list is still complete! */ 00200 00201 fprintf( ioQQQ, "\n I will now list all stellar atmosphere grids that are ready to be used (if any).\n" ); 00202 fprintf( ioQQQ, " User-defined stellar atmosphere grids will not be included in this list.\n\n" ); 00203 00204 process_counter dum; 00205 00206 /* we always look in the data directory regardless of where we are, 00207 * it would be very confusing to the user if we did otherwise... */ 00208 access_scheme as = AS_DATA_ONLY_TRY; 00209 00210 if( lgValidBinFile( "atlas_fp10k2.mod", dum, as ) ) 00211 fprintf( ioQQQ, " table star atlas Z+1.0 <Teff> [ <log(g)> ]\n" ); 00212 if( lgValidBinFile( "atlas_fp05k2.mod", dum, as ) ) 00213 fprintf( ioQQQ, " table star atlas Z+0.5 <Teff> [ <log(g)> ]\n" ); 00214 if( lgValidBinFile( "atlas_fp03k2.mod", dum, as ) ) 00215 fprintf( ioQQQ, " table star atlas Z+0.3 <Teff> [ <log(g)> ]\n" ); 00216 if( lgValidBinFile( "atlas_fp02k2.mod", dum, as ) ) 00217 fprintf( ioQQQ, " table star atlas Z+0.2 <Teff> [ <log(g)> ]\n" ); 00218 if( lgValidBinFile( "atlas_fp01k2.mod", dum, as ) ) 00219 fprintf( ioQQQ, " table star atlas Z+0.1 <Teff> [ <log(g)> ]\n" ); 00220 if( lgValidBinFile( "atlas_fp00k2.mod", dum, as ) ) 00221 fprintf( ioQQQ, " table star atlas Z+0.0 <Teff> [ <log(g)> ]\n" ); 00222 if( lgValidBinFile( "atlas_fm01k2.mod", dum, as ) ) 00223 fprintf( ioQQQ, " table star atlas Z-0.1 <Teff> [ <log(g)> ]\n" ); 00224 if( lgValidBinFile( "atlas_fm02k2.mod", dum, as ) ) 00225 fprintf( ioQQQ, " table star atlas Z-0.2 <Teff> [ <log(g)> ]\n" ); 00226 if( lgValidBinFile( "atlas_fm03k2.mod", dum, as ) ) 00227 fprintf( ioQQQ, " table star atlas Z-0.3 <Teff> [ <log(g)> ]\n" ); 00228 if( lgValidBinFile( "atlas_fm05k2.mod", dum, as ) ) 00229 fprintf( ioQQQ, " table star atlas Z-0.5 <Teff> [ <log(g)> ]\n" ); 00230 if( lgValidBinFile( "atlas_fm10k2.mod", dum, as ) ) 00231 fprintf( ioQQQ, " table star atlas Z-1.0 <Teff> [ <log(g)> ]\n" ); 00232 if( lgValidBinFile( "atlas_fm15k2.mod", dum, as ) ) 00233 fprintf( ioQQQ, " table star atlas Z-1.5 <Teff> [ <log(g)> ]\n" ); 00234 if( lgValidBinFile( "atlas_fm20k2.mod", dum, as ) ) 00235 fprintf( ioQQQ, " table star atlas Z-2.0 <Teff> [ <log(g)> ]\n" ); 00236 if( lgValidBinFile( "atlas_fm25k2.mod", dum, as ) ) 00237 fprintf( ioQQQ, " table star atlas Z-2.5 <Teff> [ <log(g)> ]\n" ); 00238 if( lgValidBinFile( "atlas_fm30k2.mod", dum, as ) ) 00239 fprintf( ioQQQ, " table star atlas Z-3.0 <Teff> [ <log(g)> ]\n" ); 00240 if( lgValidBinFile( "atlas_fm35k2.mod", dum, as ) ) 00241 fprintf( ioQQQ, " table star atlas Z-3.5 <Teff> [ <log(g)> ]\n" ); 00242 if( lgValidBinFile( "atlas_fm40k2.mod", dum, as ) ) 00243 fprintf( ioQQQ, " table star atlas Z-4.0 <Teff> [ <log(g)> ]\n" ); 00244 if( lgValidBinFile( "atlas_fm45k2.mod", dum, as ) ) 00245 fprintf( ioQQQ, " table star atlas Z-4.5 <Teff> [ <log(g)> ]\n" ); 00246 if( lgValidBinFile( "atlas_fm50k2.mod", dum, as ) ) 00247 fprintf( ioQQQ, " table star atlas Z-5.0 <Teff> [ <log(g)> ]\n" ); 00248 00249 if( lgValidBinFile( "atlas_fp05k2_odfnew.mod", dum, as ) ) 00250 fprintf( ioQQQ, " table star atlas odfnew Z+0.5 <Teff> [ <log(g)> ]\n" ); 00251 if( lgValidBinFile( "atlas_fp02k2_odfnew.mod", dum, as ) ) 00252 fprintf( ioQQQ, " table star atlas odfnew Z+0.2 <Teff> [ <log(g)> ]\n" ); 00253 if( lgValidBinFile( "atlas_fp00k2_odfnew.mod", dum, as ) ) 00254 fprintf( ioQQQ, " table star atlas odfnew Z+0.0 <Teff> [ <log(g)> ]\n" ); 00255 if( lgValidBinFile( "atlas_fm05k2_odfnew.mod", dum, as ) ) 00256 fprintf( ioQQQ, " table star atlas odfnew Z-0.5 <Teff> [ <log(g)> ]\n" ); 00257 if( lgValidBinFile( "atlas_fm10k2_odfnew.mod", dum, as ) ) 00258 fprintf( ioQQQ, " table star atlas odfnew Z-1.0 <Teff> [ <log(g)> ]\n" ); 00259 if( lgValidBinFile( "atlas_fm15k2_odfnew.mod", dum, as ) ) 00260 fprintf( ioQQQ, " table star atlas odfnew Z-1.5 <Teff> [ <log(g)> ]\n" ); 00261 if( lgValidBinFile( "atlas_fm20k2_odfnew.mod", dum, as ) ) 00262 fprintf( ioQQQ, " table star atlas odfnew Z-2.0 <Teff> [ <log(g)> ]\n" ); 00263 if( lgValidBinFile( "atlas_fm25k2_odfnew.mod", dum, as ) ) 00264 fprintf( ioQQQ, " table star atlas odfnew Z-2.5 <Teff> [ <log(g)> ]\n" ); 00265 00266 if( lgValidBinFile( "atlas_3d.mod", dum, as ) ) 00267 fprintf( ioQQQ, " table star atlas 3-dim <Teff> <log(g)> <log(Z)>\n" ); 00268 00269 if( lgValidBinFile( "atlas_3d_odfnew.mod", dum, as ) ) 00270 fprintf( ioQQQ, " table star atlas odfnew 3-dim <Teff> <log(g)> <log(Z)>\n" ); 00271 00272 if( lgValidBinFile( "Sc1_costar_solar.mod", dum, as ) ) 00273 fprintf( ioQQQ, " table star costar solar (see Hazy for parameters)\n" ); 00274 if( lgValidBinFile( "Sc1_costar_halo.mod", dum, as ) ) 00275 fprintf( ioQQQ, " table star costar halo (see Hazy for parameters)\n" ); 00276 00277 if( lgValidBinFile( "kurucz79.mod", dum, as ) ) 00278 fprintf( ioQQQ, " table star kurucz79 <Teff>\n" ); 00279 00280 if( lgValidBinFile( "mihalas.mod", dum, as ) ) 00281 fprintf( ioQQQ, " table star mihalas <Teff>\n" ); 00282 00283 if( lgValidBinFile( "rauch_h-ca_solar.mod", dum, as ) ) 00284 fprintf( ioQQQ, " table star rauch H-Ca solar <Teff> [ <log(g)> ]\n" ); 00285 if( lgValidBinFile( "rauch_h-ca_halo.mod", dum, as ) ) 00286 fprintf( ioQQQ, " table star rauch H-Ca halo <Teff> [ <log(g)> ]\n" ); 00287 if( lgValidBinFile( "rauch_h-ca_3d.mod", dum, as ) ) 00288 fprintf( ioQQQ, " table star rauch H-Ca 3-dim <Teff> <log(g)> <log(Z)>\n" ); 00289 00290 if( lgValidBinFile( "rauch_h-ni_solar.mod", dum, as ) ) 00291 fprintf( ioQQQ, " table star rauch H-Ni solar <Teff> [ <log(g)> ]\n" ); 00292 if( lgValidBinFile( "rauch_h-ni_halo.mod", dum, as ) ) 00293 fprintf( ioQQQ, " table star rauch H-Ni halo <Teff> [ <log(g)> ]\n" ); 00294 if( lgValidBinFile( "rauch_h-ni_3d.mod", dum, as ) ) 00295 fprintf( ioQQQ, " table star rauch H-Ni 3-dim <Teff> <log(g)> <log(Z)>\n" ); 00296 00297 if( lgValidBinFile( "rauch_pg1159.mod", dum, as ) ) 00298 fprintf( ioQQQ, " table star rauch pg1159 <Teff> [ <log(g)> ]\n" ); 00299 00300 if( lgValidBinFile( "rauch_hydr.mod", dum, as ) ) 00301 fprintf( ioQQQ, " table star rauch hydrogen <Teff> [ <log(g)> ]\n" ); 00302 00303 if( lgValidBinFile( "rauch_helium.mod", dum, as ) ) 00304 fprintf( ioQQQ, " table star rauch helium <Teff> [ <log(g)> ]\n" ); 00305 00306 if( lgValidBinFile( "rauch_h+he_3d.mod", dum, as ) ) 00307 fprintf( ioQQQ, " table star rauch H+He <Teff> <log(g)> <frac(He)>\n" ); 00308 00309 if( lgValidBinFile( "starburst99.mod", dum, as ) ) 00310 fprintf( ioQQQ, " table star \"starburst99.mod\" <age>\n" ); 00311 00312 if( lgValidBinFile( "bstar2006_p03.mod", dum, as ) ) 00313 fprintf( ioQQQ, " table star tlusty Bstar Z+0.3 <Teff> [ <log(g)> ]\n" ); 00314 if( lgValidBinFile( "bstar2006_p00.mod", dum, as ) ) 00315 fprintf( ioQQQ, " table star tlusty Bstar Z+0.0 <Teff> [ <log(g)> ]\n" ); 00316 if( lgValidBinFile( "bstar2006_m03.mod", dum, as ) ) 00317 fprintf( ioQQQ, " table star tlusty Bstar Z-0.3 <Teff> [ <log(g)> ]\n" ); 00318 if( lgValidBinFile( "bstar2006_m07.mod", dum, as ) ) 00319 fprintf( ioQQQ, " table star tlusty Bstar Z-0.7 <Teff> [ <log(g)> ]\n" ); 00320 if( lgValidBinFile( "bstar2006_m10.mod", dum, as ) ) 00321 fprintf( ioQQQ, " table star tlusty Bstar Z-1.0 <Teff> [ <log(g)> ]\n" ); 00322 if( lgValidBinFile( "bstar2006_m99.mod", dum, as ) ) 00323 fprintf( ioQQQ, " table star tlusty Bstar Z-inf <Teff> [ <log(g)> ]\n" ); 00324 00325 if( lgValidBinFile( "bstar2006_3d.mod", dum, as ) ) 00326 fprintf( ioQQQ, " table star tlusty Bstar 3-dim <Teff> <log(g)> <log(Z)>\n" ); 00327 00328 if( lgValidBinFile( "ostar2002_p03.mod", dum, as ) ) 00329 fprintf( ioQQQ, " table star tlusty Ostar Z+0.3 <Teff> [ <log(g)> ]\n" ); 00330 if( lgValidBinFile( "ostar2002_p00.mod", dum, as ) ) 00331 fprintf( ioQQQ, " table star tlusty Ostar Z+0.0 <Teff> [ <log(g)> ]\n" ); 00332 if( lgValidBinFile( "ostar2002_m03.mod", dum, as ) ) 00333 fprintf( ioQQQ, " table star tlusty Ostar Z-0.3 <Teff> [ <log(g)> ]\n" ); 00334 if( lgValidBinFile( "ostar2002_m07.mod", dum, as ) ) 00335 fprintf( ioQQQ, " table star tlusty Ostar Z-0.7 <Teff> [ <log(g)> ]\n" ); 00336 if( lgValidBinFile( "ostar2002_m10.mod", dum, as ) ) 00337 fprintf( ioQQQ, " table star tlusty Ostar Z-1.0 <Teff> [ <log(g)> ]\n" ); 00338 if( lgValidBinFile( "ostar2002_m15.mod", dum, as ) ) 00339 fprintf( ioQQQ, " table star tlusty Ostar Z-1.5 <Teff> [ <log(g)> ]\n" ); 00340 if( lgValidBinFile( "ostar2002_m17.mod", dum, as ) ) 00341 fprintf( ioQQQ, " table star tlusty Ostar Z-1.7 <Teff> [ <log(g)> ]\n" ); 00342 if( lgValidBinFile( "ostar2002_m20.mod", dum, as ) ) 00343 fprintf( ioQQQ, " table star tlusty Ostar Z-2.0 <Teff> [ <log(g)> ]\n" ); 00344 if( lgValidBinFile( "ostar2002_m30.mod", dum, as ) ) 00345 fprintf( ioQQQ, " table star tlusty Ostar Z-3.0 <Teff> [ <log(g)> ]\n" ); 00346 if( lgValidBinFile( "ostar2002_m99.mod", dum, as ) ) 00347 fprintf( ioQQQ, " table star tlusty Ostar Z-inf <Teff> [ <log(g)> ]\n" ); 00348 00349 if( lgValidBinFile( "ostar2002_3d.mod", dum, as ) ) 00350 fprintf( ioQQQ, " table star tlusty Ostar 3-dim <Teff> <log(g)> <log(Z)>\n" ); 00351 00352 if( lgValidBinFile( "kwerner.mod", dum, as ) ) 00353 fprintf( ioQQQ, " table star werner <Teff> [ <log(g)> ]\n" ); 00354 00355 if( lgValidBinFile( "wmbasic.mod", dum, as ) ) 00356 fprintf( ioQQQ, " table star wmbasic <Teff> <log(g)> <log(Z)>\n" ); 00357 return; 00358 } 00359 00360 /* AtlasCompile rebin Kurucz stellar models to match energy grid of code */ 00361 /* >>chng 05 nov 16, added return value to indicate success (0) or failure (1) */ 00362 int AtlasCompile(process_counter& pc) 00363 { 00364 /* these contain frequencies for the major absorption edges */ 00365 realnum Edges[3]; 00366 00367 bool lgFail = false; 00368 00369 DEBUG_ENTRY( "AtlasCompile()" ); 00370 00371 /* This is a program to re-bin the Kurucz stellar models spectrum to match the 00372 * CLOUDY grid. For wavelengths shorter than supplied in the Kurucz files, 00373 * the flux will be set to zero. At long wavelengths a Rayleigh-Jeans 00374 * extrapolation will be used. */ 00375 00376 /* This version uses power-law interpolation between the points of the stellar 00377 * model.*/ 00378 00379 fprintf( ioQQQ, " AtlasCompile on the job.\n" ); 00380 00381 /* define the major absorption edges that require special attention during rebinning 00382 * 00383 * NB the frequencies should be chosen here such that they are somewhere in between 00384 * the two frequency points that straddle the edge in the atmosphere model, the 00385 * software in RebinAtmosphere will seek out the exact values of those two points 00386 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at 00387 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A). 00388 * 00389 * NB beware not to choose edges too close to one another (i.e. on the order of the 00390 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides 00391 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will 00392 * almost certainly lead to erroneous behaviour in RebinAtmosphere */ 00393 Edges[0] = (realnum)(RYDLAM/911.76); 00394 Edges[1] = (realnum)(RYDLAM/504.26); 00395 Edges[2] = (realnum)(RYDLAM/227.84); 00396 00397 access_scheme as = AS_LOCAL_ONLY_TRY; 00398 00399 /* >>chng 05 nov 19, add support for non-solar metalicities as well as odfnew models, PvH */ 00400 if( lgFileReadable( "atlas_fp10k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp10k2.mod", pc, as ) ) 00401 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp10k2.ascii", "atlas_fp10k2.mod", Edges, 3L, pc ); 00402 if( lgFileReadable( "atlas_fp05k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp05k2.mod", pc, as ) ) 00403 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp05k2.ascii", "atlas_fp05k2.mod", Edges, 3L, pc ); 00404 if( lgFileReadable( "atlas_fp03k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp03k2.mod", pc, as ) ) 00405 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp03k2.ascii", "atlas_fp03k2.mod", Edges, 3L, pc ); 00406 if( lgFileReadable( "atlas_fp02k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp02k2.mod", pc, as ) ) 00407 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp02k2.ascii", "atlas_fp02k2.mod", Edges, 3L, pc ); 00408 if( lgFileReadable( "atlas_fp01k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp01k2.mod", pc, as ) ) 00409 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp01k2.ascii", "atlas_fp01k2.mod", Edges, 3L, pc ); 00410 if( lgFileReadable( "atlas_fp00k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp00k2.mod", pc, as ) ) 00411 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp00k2.ascii", "atlas_fp00k2.mod", Edges, 3L, pc ); 00412 if( lgFileReadable( "atlas_fm01k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm01k2.mod", pc, as ) ) 00413 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm01k2.ascii", "atlas_fm01k2.mod", Edges, 3L, pc ); 00414 if( lgFileReadable( "atlas_fm02k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm02k2.mod", pc, as ) ) 00415 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm02k2.ascii", "atlas_fm02k2.mod", Edges, 3L, pc ); 00416 if( lgFileReadable( "atlas_fm03k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm03k2.mod", pc, as ) ) 00417 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm03k2.ascii", "atlas_fm03k2.mod", Edges, 3L, pc ); 00418 if( lgFileReadable( "atlas_fm05k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm05k2.mod", pc, as ) ) 00419 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm05k2.ascii", "atlas_fm05k2.mod", Edges, 3L, pc ); 00420 if( lgFileReadable( "atlas_fm10k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm10k2.mod", pc, as ) ) 00421 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm10k2.ascii", "atlas_fm10k2.mod", Edges, 3L, pc ); 00422 if( lgFileReadable( "atlas_fm15k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm15k2.mod", pc, as ) ) 00423 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm15k2.ascii", "atlas_fm15k2.mod", Edges, 3L, pc ); 00424 if( lgFileReadable( "atlas_fm20k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm20k2.mod", pc, as ) ) 00425 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm20k2.ascii", "atlas_fm20k2.mod", Edges, 3L, pc ); 00426 if( lgFileReadable( "atlas_fm25k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm25k2.mod", pc, as ) ) 00427 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm25k2.ascii", "atlas_fm25k2.mod", Edges, 3L, pc ); 00428 if( lgFileReadable( "atlas_fm30k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm30k2.mod", pc, as ) ) 00429 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm30k2.ascii", "atlas_fm30k2.mod", Edges, 3L, pc ); 00430 if( lgFileReadable( "atlas_fm35k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm35k2.mod", pc, as ) ) 00431 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm35k2.ascii", "atlas_fm35k2.mod", Edges, 3L, pc ); 00432 if( lgFileReadable( "atlas_fm40k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm40k2.mod", pc, as ) ) 00433 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm40k2.ascii", "atlas_fm40k2.mod", Edges, 3L, pc ); 00434 if( lgFileReadable( "atlas_fm45k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm45k2.mod", pc, as ) ) 00435 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm45k2.ascii", "atlas_fm45k2.mod", Edges, 3L, pc ); 00436 if( lgFileReadable( "atlas_fm50k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm50k2.mod", pc, as ) ) 00437 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm50k2.ascii", "atlas_fm50k2.mod", Edges, 3L, pc ); 00438 00439 if( lgFileReadable( "atlas_fp05k2_odfnew.ascii", pc, as ) && 00440 !lgValidBinFile( "atlas_fp05k2_odfnew.mod", pc, as ) ) 00441 00442 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp05k2_odfnew.ascii", "atlas_fp05k2_odfnew.mod", 00443 Edges, 3L, pc ); 00444 if( lgFileReadable( "atlas_fp02k2_odfnew.ascii", pc, as ) && 00445 !lgValidBinFile( "atlas_fp02k2_odfnew.mod", pc, as ) ) 00446 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp02k2_odfnew.ascii", "atlas_fp02k2_odfnew.mod", 00447 Edges, 3L, pc ); 00448 if( lgFileReadable( "atlas_fp00k2_odfnew.ascii", pc, as ) && 00449 !lgValidBinFile( "atlas_fp00k2_odfnew.mod", pc, as ) ) 00450 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp00k2_odfnew.ascii", "atlas_fp00k2_odfnew.mod", 00451 Edges, 3L, pc ); 00452 if( lgFileReadable( "atlas_fm05k2_odfnew.ascii", pc, as ) && 00453 !lgValidBinFile( "atlas_fm05k2_odfnew.mod", pc, as ) ) 00454 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm05k2_odfnew.ascii", "atlas_fm05k2_odfnew.mod", 00455 Edges, 3L, pc ); 00456 if( lgFileReadable( "atlas_fm10k2_odfnew.ascii", pc, as ) && 00457 !lgValidBinFile( "atlas_fm10k2_odfnew.mod", pc, as ) ) 00458 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm10k2_odfnew.ascii", "atlas_fm10k2_odfnew.mod", 00459 Edges, 3L, pc ); 00460 if( lgFileReadable( "atlas_fm15k2_odfnew.ascii", pc, as ) && 00461 !lgValidBinFile( "atlas_fm15k2_odfnew.mod", pc, as ) ) 00462 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm15k2_odfnew.ascii", "atlas_fm15k2_odfnew.mod", 00463 Edges, 3L, pc ); 00464 if( lgFileReadable( "atlas_fm20k2_odfnew.ascii", pc, as ) && 00465 !lgValidBinFile( "atlas_fm20k2_odfnew.mod", pc, as ) ) 00466 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm20k2_odfnew.ascii", "atlas_fm20k2_odfnew.mod", 00467 Edges, 3L, pc ); 00468 if( lgFileReadable( "atlas_fm25k2_odfnew.ascii", pc, as ) && 00469 !lgValidBinFile( "atlas_fm25k2_odfnew.mod", pc, as ) ) 00470 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm25k2_odfnew.ascii", "atlas_fm25k2_odfnew.mod", 00471 Edges, 3L, pc ); 00472 00473 if( lgFileReadable( "atlas_3d.ascii", pc, as ) && !lgValidBinFile( "atlas_3d.mod", pc, as ) ) 00474 lgFail = lgFail || lgCompileAtmosphere( "atlas_3d.ascii", "atlas_3d.mod", Edges, 3L, pc ); 00475 00476 if( lgFileReadable( "atlas_3d_odfnew.ascii", pc, as ) && 00477 !lgValidBinFile( "atlas_3d_odfnew.mod", pc, as ) ) 00478 lgFail = lgFail || lgCompileAtmosphere( "atlas_3d_odfnew.ascii", "atlas_3d_odfnew.mod", Edges, 3L, pc ); 00479 return lgFail; 00480 } 00481 00482 /* AtlasInterpolate read in and interpolate on Kurucz grid of atmospheres, originally by K Volk */ 00483 long AtlasInterpolate(double val[], /* val[nval] */ 00484 long *nval, 00485 long *ndim, 00486 const char *chMetalicity, 00487 const char *chODFNew, 00488 bool lgList, 00489 double *Tlow, 00490 double *Thigh) 00491 { 00492 char chIdent[13]; 00493 stellar_grid grid; 00494 00495 DEBUG_ENTRY( "AtlasInterpolate()" ); 00496 00497 grid.name = "atlas_"; 00498 if( *ndim == 3 ) 00499 grid.name += "3d"; 00500 else 00501 { 00502 grid.name += "f"; 00503 grid.name += chMetalicity; 00504 grid.name += "k2"; 00505 } 00506 grid.name += chODFNew; 00507 grid.name += ".mod"; 00508 grid.scheme = AS_DATA_OPTIONAL; 00509 /* identification of this atmosphere set, used in 00510 * the Cloudy output, *must* be 12 characters long */ 00511 if( *ndim == 3 ) 00512 { 00513 strcpy( chIdent, "3-dim" ); 00514 } 00515 else 00516 { 00517 strcpy( chIdent, "Z " ); 00518 strcat( chIdent, chMetalicity ); 00519 } 00520 strcat( chIdent, ( strlen(chODFNew) == 0 ? " Kurucz" : " ODFNew" ) ); 00521 grid.ident = chIdent; 00522 /* the Cloudy command needed to recompile the binary model file */ 00523 grid.command = "COMPILE STARS"; 00524 00525 InitGrid( &grid, lgList ); 00526 00527 CheckVal( &grid, val, nval, ndim ); 00528 00529 /* Note on the interpolation (solar abundance grid): 26 October 2000 (Peter van Hoof) 00530 * 00531 * I computed the effective temperature for a random sample of interpolated 00532 * atmospheres by integrating the flux as shown above and compared the results 00533 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC. 00534 * 00535 * I found that the average discrepancy was: 00536 * 00537 * DELTA = -0.10% +/- 0.06% (sample size 5000) 00538 * 00539 * The most extreme discrepancies were 00540 * -0.30% <= DELTA <= 0.21% 00541 * 00542 * The most negative discrepancies were for Teff = 36 - 39 kK, log(g) = 4.5 - 5 00543 * The most positive discrepancies were for Teff = 3.5 - 4.0 kK, log(g) = 0 - 1 00544 * 00545 * The interpolation in the ATLAS grid is clearly very accurate */ 00546 00547 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 00548 00549 FreeGrid( &grid ); 00550 return rfield.nupper; 00551 } 00552 00553 /* CoStarCompile rebin costar stellar models to match energy grid of code*/ 00554 int CoStarCompile(process_counter& pc) 00555 { 00556 realnum Edges[3]; 00557 bool lgFail = false; 00558 00559 DEBUG_ENTRY( "CoStarCompile()" ); 00560 00561 fprintf( ioQQQ, " CoStarCompile on the job.\n" ); 00562 00563 /* define the major absorption edges that require special attention during rebinning 00564 * 00565 * NB the frequencies should be chosen here such that they are somewhere in between 00566 * the two frequency points that straddle the edge in the atmosphere model, the 00567 * software in RebinAtmosphere will seek out the exact values of those two points 00568 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at 00569 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A). 00570 * 00571 * NB beware not to choose edges too close to one another (i.e. on the order of the 00572 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides 00573 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will 00574 * almost certainly lead to erroneous behaviour in RebinAtmosphere */ 00575 Edges[0] = (realnum)(RYDLAM/911.76); 00576 Edges[1] = (realnum)(RYDLAM/504.26); 00577 Edges[2] = (realnum)(RYDLAM/227.84); 00578 00579 access_scheme as = AS_LOCAL_ONLY_TRY; 00580 00581 if( lgFileReadable( "Sc1_costar_z020_lb.fluxes", pc, as ) && !lgValidBinFile( "Sc1_costar_solar.mod", pc, as ) ) 00582 lgFail = lgFail || lgCompileAtmosphereCoStar( "Sc1_costar_z020_lb.fluxes", "Sc1_costar_solar.mod", 00583 Edges, 3L, pc ); 00584 if( lgFileReadable( "Sc1_costar_z004_lb.fluxes", pc, as ) && !lgValidBinFile( "Sc1_costar_halo.mod", pc, as ) ) 00585 lgFail = lgFail || lgCompileAtmosphereCoStar( "Sc1_costar_z004_lb.fluxes", "Sc1_costar_halo.mod", 00586 Edges, 3L, pc ); 00587 return lgFail; 00588 } 00589 00590 /* CoStarInterpolate read in and interpolate on CoStar grid of atmospheres */ 00591 long CoStarInterpolate(double val[], /* requested model parameters */ 00592 long *nval, 00593 long *ndim, 00594 IntMode imode, /* which interpolation mode is requested */ 00595 bool lgHalo, /* flag indicating whether solar (==0) or halo (==1) abundances */ 00596 bool lgList, 00597 double *val0_lo, 00598 double *val0_hi) 00599 { 00600 stellar_grid grid; 00601 00602 DEBUG_ENTRY( "CoStarInterpolate()" ); 00603 00604 grid.name = ( lgHalo ? "Sc1_costar_halo.mod" : "Sc1_costar_solar.mod" ); 00605 grid.scheme = AS_DATA_OPTIONAL; 00606 /* identification of this atmosphere set, used in 00607 * the Cloudy output, *must* be 12 characters long */ 00608 grid.ident = " costar"; 00609 /* the Cloudy command needed to recompile the binary model file */ 00610 grid.command = "COMPILE STARS"; 00611 00612 /* listing the models in the grid is implemented in CoStarListModels() */ 00613 InitGrid( &grid, false ); 00614 /* now sort the models according to track */ 00615 InitGridCoStar( &grid ); 00616 /* override default interpolation mode */ 00617 grid.imode = imode; 00618 00619 if( lgList ) 00620 { 00621 CoStarListModels( &grid ); 00622 cdEXIT(EXIT_SUCCESS); 00623 } 00624 00625 CheckVal( &grid, val, nval, ndim ); 00626 00627 /* Note on the interpolation: 26 October 2000 (Peter van Hoof) 00628 * 00629 * I computed the effective temperature for a random sample of interpolated 00630 * atmospheres by integrating the flux as shown above and compared the results 00631 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC. 00632 * 00633 * I found that the average discrepancy was: 00634 * 00635 * DELTA = -1.16% +/- 0.69% (SOLAR models, sample size 4590) 00636 * DELTA = -1.17% +/- 0.70% (HALO models, sample size 4828) 00637 * 00638 * The most extreme discrepancies for the SOLAR models were 00639 * -3.18% <= DELTA <= -0.16% 00640 * 00641 * The most negative discrepancies were for Teff = 35 kK, log(g) = 3.5 00642 * The least negative discrepancies were for Teff = 50 kK, log(g) = 4.1 00643 * 00644 * The most extreme discrepancies for the HALO models were 00645 * -2.90% <= DELTA <= -0.13% 00646 * 00647 * The most negative discrepancies were for Teff = 35 kK, log(g) = 3.5 00648 * The least negative discrepancies were for Teff = 50 kK, log(g) = 4.1 00649 * 00650 * Since Cloudy checks the scaling elsewhere there is no need to re-scale 00651 * things here, but this inaccuracy should be kept in mind since it could 00652 * indicate problems with the flux distribution */ 00653 00654 InterpolateGridCoStar( &grid, val, val0_lo, val0_hi ); 00655 00656 FreeGrid( &grid ); 00657 return rfield.nupper; 00658 } 00659 00660 /* GridCompile rebin user supplied stellar models to match energy grid of code */ 00661 bool GridCompile(const char *InName) 00662 { 00663 bool lgFail = false; 00664 realnum Edges[1]; 00665 string OutName( InName ); 00666 stellar_grid grid; 00667 00668 DEBUG_ENTRY( "GridCompile()" ); 00669 00670 fprintf( ioQQQ, " GridCompile on the job.\n" ); 00671 00672 // replace filename extension with ".mod" 00673 string::size_type ptr = OutName.find( '.' ); 00674 ASSERT( ptr != string::npos ); 00675 OutName.replace( ptr, string::npos, ".mod" ); 00676 00677 process_counter dum; 00678 lgFail = lgCompileAtmosphere( InName, OutName.c_str(), Edges, 0L, dum ); 00679 00680 if( !lgFail ) 00681 { 00682 /* the file must be local */ 00683 grid.name = OutName; 00684 grid.scheme = AS_LOCAL_ONLY; 00685 grid.ident = "bogus ident."; 00686 grid.command = "bogus command."; 00687 00688 InitGrid( &grid, false ); 00689 00690 /* check whether the models in the grid have the correct effective temperature */ 00691 00692 fprintf( ioQQQ, " GridCompile: checking effective temperatures...\n" ); 00693 ValidateGrid( &grid, 0.02 ); 00694 00695 FreeGrid( &grid ); 00696 } 00697 00698 return lgFail; 00699 } 00700 00701 /* GridInterpolate read in and interpolate on user supplied grid of atmospheres */ 00702 long GridInterpolate(double val[], /* val[nval] */ 00703 long *nval, 00704 long *ndim, 00705 const char *FileName, 00706 bool lgList, 00707 double *Tlow, 00708 double *Thigh) 00709 { 00710 char chIdent[13]; 00711 stellar_grid grid; 00712 00713 DEBUG_ENTRY( "GridInterpolate()" ); 00714 00715 // make filename without extension 00716 string chTruncName( FileName ); 00717 string::size_type ptr = chTruncName.find( '.' ); 00718 if( ptr != string::npos ) 00719 chTruncName.replace( ptr, string::npos, "" ); 00720 00721 grid.name = FileName; 00722 grid.scheme = AS_DATA_OPTIONAL; 00723 /* identification of this atmosphere set, used in 00724 * the Cloudy output, *must* be 12 characters long */ 00725 sprintf( chIdent, "%12.12s", chTruncName.c_str() ); 00726 grid.ident = chIdent; 00727 /* the Cloudy command needed to recompile the binary model file */ 00728 string chString( "COMPILE STARS \"" + chTruncName + ".ascii\"" ); 00729 grid.command = chString.c_str(); 00730 00731 InitGrid( &grid, lgList ); 00732 00733 CheckVal( &grid, val, nval, ndim ); 00734 00735 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 00736 00737 FreeGrid( &grid ); 00738 return rfield.nupper; 00739 } 00740 00741 /* Kurucz79Compile rebin Kurucz 1979 stellar models to match energy grid of code */ 00742 int Kurucz79Compile(process_counter& pc) 00743 { 00744 realnum Edges[1]; 00745 00746 bool lgFail = false; 00747 00748 DEBUG_ENTRY( "Kurucz79Compile()" ); 00749 00750 fprintf( ioQQQ, " Kurucz79Compile on the job.\n" ); 00751 00752 /* following atmospheres LTE from Kurucz 1979, Ap.J. Sup 40, 1. and 00753 * Kurucz (1989) private communication, newer opacities */ 00754 00755 access_scheme as = AS_LOCAL_ONLY_TRY; 00756 00757 if( lgFileReadable( "kurucz79.ascii", pc, as ) && !lgValidBinFile( "kurucz79.mod", pc, as ) ) 00758 lgFail = lgCompileAtmosphere( "kurucz79.ascii", "kurucz79.mod", Edges, 0L, pc ); 00759 return lgFail; 00760 } 00761 00762 /* Kurucz79Interpolate read in and interpolate on Kurucz79 grid of atmospheres */ 00763 long Kurucz79Interpolate(double val[], /* val[nval] */ 00764 long *nval, 00765 long *ndim, 00766 bool lgList, 00767 double *Tlow, 00768 double *Thigh) 00769 { 00770 stellar_grid grid; 00771 00772 DEBUG_ENTRY( "Kurucz79Interpolate()" ); 00773 00774 grid.name = "kurucz79.mod"; 00775 grid.scheme = AS_DATA_OPTIONAL; 00776 /* identification of this atmosphere set, used in 00777 * the Cloudy output, *must* be 12 characters long */ 00778 grid.ident = " Kurucz 1979"; 00779 /* the Cloudy command needed to recompile the binary model file */ 00780 grid.command = "COMPILE STARS"; 00781 00782 InitGrid( &grid, lgList ); 00783 00784 CheckVal( &grid, val, nval, ndim ); 00785 00786 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 00787 00788 FreeGrid( &grid ); 00789 return rfield.nupper; 00790 } 00791 00792 /* MihalasCompile rebin Mihalas stellar models to match energy grid of code */ 00793 int MihalasCompile(process_counter& pc) 00794 { 00795 realnum Edges[1]; 00796 00797 bool lgFail = false; 00798 00799 DEBUG_ENTRY( "MihalasCompile()" ); 00800 00801 fprintf( ioQQQ, " MihalasCompile on the job.\n" ); 00802 00803 /* following atmospheres NLTE from Mihalas, NCAR-TN/STR-76 */ 00804 00805 access_scheme as = AS_LOCAL_ONLY_TRY; 00806 00807 if( lgFileReadable( "mihalas.ascii", pc, as ) && !lgValidBinFile( "mihalas.mod", pc, as ) ) 00808 lgFail = lgCompileAtmosphere( "mihalas.ascii", "mihalas.mod", Edges, 0L, pc ); 00809 return lgFail; 00810 } 00811 00812 /* MihalasInterpolate read in and interpolate on Mihalas grid of atmospheres */ 00813 long MihalasInterpolate(double val[], /* val[nval] */ 00814 long *nval, 00815 long *ndim, 00816 bool lgList, 00817 double *Tlow, 00818 double *Thigh) 00819 { 00820 stellar_grid grid; 00821 00822 DEBUG_ENTRY( "MihalasInterpolate()" ); 00823 00824 grid.name = "mihalas.mod"; 00825 grid.scheme = AS_DATA_OPTIONAL; 00826 /* identification of this atmosphere set, used in 00827 * the Cloudy output, *must* be 12 characters long */ 00828 grid.ident = " Mihalas"; 00829 /* the Cloudy command needed to recompile the binary model file */ 00830 grid.command = "COMPILE STARS"; 00831 00832 InitGrid( &grid, lgList ); 00833 00834 CheckVal( &grid, val, nval, ndim ); 00835 00836 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 00837 00838 FreeGrid( &grid ); 00839 return rfield.nupper; 00840 } 00841 00842 /* RauchCompile create ascii and mod files for Rauch atmospheres */ 00843 int RauchCompile(process_counter& pc) 00844 { 00845 bool lgFail = false; 00846 00847 /* these contain frequencies for the major absorption edges */ 00848 realnum Edges[3]; 00849 00850 /* Before running this program issue the following command where the Rauch 00851 * model atmosphere files are kept (0050000_50_solar_bin_0.1 and so on) 00852 * 00853 * ls *solar_bin_0.1 > rauchmods.list 00854 * 00855 * and check to see that there are 66 lines in the file. 00856 */ 00857 00858 static const mpp telg1[NMODS_HCA] = { 00859 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, 00860 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, 00861 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, 00862 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, 00863 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, 00864 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}}, 00865 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}}, 00866 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}}, 00867 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}}, 00868 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}}, 00869 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}}, 00870 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}}, 00871 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}}, 00872 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}}, 00873 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}}, 00874 {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}}, 00875 {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}}, 00876 {{400000.,8.0}},{{400000.,9.0}}, 00877 {{500000.,8.0}},{{500000.,9.0}}, 00878 {{600000.,9.0}}, 00879 {{700000.,9.0}}, 00880 {{800000.,9.0}}, 00881 {{900000.,9.0}}, 00882 {{1000000.,9.0}} 00883 }; 00884 00885 static const mpp telg2[NMODS_HNI] = { 00886 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, 00887 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, 00888 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, 00889 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, 00890 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, 00891 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}}, 00892 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}}, 00893 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}}, 00894 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}}, 00895 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}}, 00896 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}}, 00897 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}}, 00898 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}}, 00899 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}}, 00900 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}} 00901 }; 00902 00903 static const mpp telg3[NMODS_PG1159] = { 00904 {{40000.,5.0}}, {{40000.,6.0}}, {{40000.,7.0}}, {{40000.,8.0}}, {{40000.,9.0}}, 00905 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}}, 00906 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}}, 00907 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}}, 00908 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}}, 00909 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}}, 00910 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}}, 00911 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}}, 00912 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}}, 00913 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}}, 00914 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}}, 00915 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}}, 00916 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}}, 00917 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}}, 00918 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}}, 00919 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}} 00920 }; 00921 00922 static const mpp telg4[NMODS_HYDR] = { 00923 {{20000.,4.0}}, {{20000.,5.0}}, {{20000.,6.0}}, {{20000.,7.0}}, {{20000.,8.0}}, {{20000.,9.0}}, 00924 {{30000.,4.0}}, {{30000.,5.0}}, {{30000.,6.0}}, {{30000.,7.0}}, {{30000.,8.0}}, {{30000.,9.0}}, 00925 {{40000.,4.0}}, {{40000.,5.0}}, {{40000.,6.0}}, {{40000.,7.0}}, {{40000.,8.0}}, {{40000.,9.0}}, 00926 {{50000.,4.0}}, {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}}, 00927 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}}, 00928 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}}, 00929 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}}, 00930 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}}, 00931 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}}, 00932 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}}, 00933 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}}, 00934 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}}, 00935 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}}, 00936 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}}, 00937 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}}, 00938 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}}, 00939 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}}, 00940 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}}, 00941 {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}}, 00942 {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}}, 00943 {{400000.,8.0}},{{400000.,9.0}}, 00944 {{500000.,8.0}},{{500000.,9.0}}, 00945 {{600000.,9.0}}, 00946 {{700000.,9.0}}, 00947 {{800000.,9.0}}, 00948 {{900000.,9.0}}, 00949 {{1000000.,9.0}} 00950 }; 00951 00952 static const mpp telg5[NMODS_HELIUM] = { 00953 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}}, 00954 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}}, 00955 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}}, 00956 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}}, 00957 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}}, 00958 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}}, 00959 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}}, 00960 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}}, 00961 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}}, 00962 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}}, 00963 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}}, 00964 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}}, 00965 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}}, 00966 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}}, 00967 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}}, 00968 {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}}, 00969 {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}}, 00970 {{400000.,8.0}},{{400000.,9.0}}, 00971 {{500000.,8.0}},{{500000.,9.0}}, 00972 {{600000.,9.0}}, 00973 {{700000.,9.0}}, 00974 {{800000.,9.0}}, 00975 {{900000.,9.0}}, 00976 {{1000000.,9.0}} 00977 }; 00978 00979 static const mpp telg6[NMODS_HpHE] = { 00980 {{50000.,5.0}}, {{50000.,5.5}}, {{50000.,6.0}}, {{50000.,6.5}}, {{50000.,7.0}}, {{50000.,7.5}}, {{50000.,8.0}}, {{50000.,8.5}}, {{50000.,9.0}}, 00981 {{60000.,5.0}}, {{60000.,5.5}}, {{60000.,6.0}}, {{60000.,6.5}}, {{60000.,7.0}}, {{60000.,7.5}}, {{60000.,8.0}}, {{60000.,8.5}}, {{60000.,9.0}}, 00982 {{70000.,5.0}}, {{70000.,5.5}}, {{70000.,6.0}}, {{70000.,6.5}}, {{70000.,7.0}}, {{70000.,7.5}}, {{70000.,8.0}}, {{70000.,8.5}}, {{70000.,9.0}}, 00983 {{80000.,5.0}}, {{80000.,5.5}}, {{80000.,6.0}}, {{80000.,6.5}}, {{80000.,7.0}}, {{80000.,7.5}}, {{80000.,8.0}}, {{80000.,8.5}}, {{80000.,9.0}}, 00984 {{90000.,5.0}}, {{90000.,5.5}}, {{90000.,6.0}}, {{90000.,6.5}}, {{90000.,7.0}}, {{90000.,7.5}}, {{90000.,8.0}}, {{90000.,8.5}}, {{90000.,9.0}}, 00985 {{100000.,5.0}},{{100000.,5.5}},{{100000.,6.0}},{{100000.,6.5}},{{100000.,7.0}},{{100000.,7.5}},{{100000.,8.0}},{{100000.,8.5}},{{100000.,9.0}}, 00986 {{110000.,6.0}},{{110000.,6.5}},{{110000.,7.0}},{{110000.,7.5}},{{110000.,8.0}},{{110000.,8.5}},{{110000.,9.0}}, 00987 {{120000.,6.0}},{{120000.,6.5}},{{120000.,7.0}},{{120000.,7.5}},{{120000.,8.0}},{{120000.,8.5}},{{120000.,9.0}}, 00988 {{130000.,6.0}},{{130000.,6.5}},{{130000.,7.0}},{{130000.,7.5}},{{130000.,8.0}},{{130000.,8.5}},{{130000.,9.0}}, 00989 {{140000.,6.0}},{{140000.,6.5}},{{140000.,7.0}},{{140000.,7.5}},{{140000.,8.0}},{{140000.,8.5}},{{140000.,9.0}}, 00990 {{150000.,6.0}},{{150000.,6.5}},{{150000.,7.0}},{{150000.,7.5}},{{150000.,8.0}},{{150000.,8.5}},{{150000.,9.0}}, 00991 {{160000.,6.0}},{{160000.,6.5}},{{160000.,7.0}},{{160000.,7.5}},{{160000.,8.0}},{{160000.,8.5}},{{160000.,9.0}}, 00992 {{170000.,6.0}},{{170000.,6.5}},{{170000.,7.0}},{{170000.,7.5}},{{170000.,8.0}},{{170000.,8.5}},{{170000.,9.0}}, 00993 {{180000.,6.0}},{{180000.,6.5}},{{180000.,7.0}},{{180000.,7.5}},{{180000.,8.0}},{{180000.,8.5}},{{180000.,9.0}}, 00994 {{190000.,6.0}},{{190000.,6.5}},{{190000.,7.0}},{{190000.,7.5}},{{190000.,8.0}},{{190000.,8.5}},{{190000.,9.0}} 00995 }; 00996 00997 /* metalicities of the solar and halo grid */ 00998 static const double par2[2] = { 0., -1. }; 00999 01000 /* Helium fraction by mass */ 01001 static const double par3[11] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 }; 01002 01003 DEBUG_ENTRY( "RauchCompile()" ); 01004 01005 fprintf( ioQQQ, " RauchCompile on the job.\n" ); 01006 01007 process_counter dum; 01008 access_scheme as = AS_LOCAL_ONLY_TRY; 01009 01010 /* this is the H-Ca grid */ 01011 if( lgFileReadable( "0050000_50_solar_bin_0.1", dum, as ) && !lgValidAsciiFile( "rauch_h-ca_solar.ascii", as ) ) 01012 { 01013 fprintf( ioQQQ, " Creating rauch_h-ca_solar.ascii....\n" ); 01014 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_solar.ascii", "_solar_bin_0.1", 01015 telg1, NMODS_HCA, 1, 1, par2, 1 ); 01016 } 01017 01018 if( lgFileReadable( "0050000_50_halo__bin_0.1", dum, as ) && !lgValidAsciiFile( "rauch_h-ca_halo.ascii", as ) ) 01019 { 01020 fprintf( ioQQQ, " Creating rauch_h-ca_halo.ascii....\n" ); 01021 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_halo.ascii", "_halo__bin_0.1", 01022 telg1, NMODS_HCA, 1, 1, par2, 1 ); 01023 } 01024 01025 if( lgFileReadable( "0050000_50_solar_bin_0.1", dum, as ) && 01026 lgFileReadable( "0050000_50_halo__bin_0.1", dum, as ) && 01027 !lgValidAsciiFile( "rauch_h-ca_3d.ascii", as ) ) 01028 { 01029 fprintf( ioQQQ, " Creating rauch_h-ca_3d.ascii....\n" ); 01030 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_3d.ascii", "_solar_bin_0.1", 01031 telg1, NMODS_HCA, 1, 2, par2, 1 ); 01032 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_3d.ascii", "_halo__bin_0.1", 01033 telg1, NMODS_HCA, 2, 2, par2, 1 ); 01034 } 01035 01036 /* this is the H-Ni grid */ 01037 if( lgFileReadable( "0050000_50_solar_iron.bin_0.1", dum, as ) && 01038 !lgValidAsciiFile( "rauch_h-ni_solar.ascii", as ) ) 01039 { 01040 fprintf( ioQQQ, " Creating rauch_h-ni_solar.ascii....\n" ); 01041 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_solar.ascii", "_solar_iron.bin_0.1", 01042 telg2, NMODS_HNI, 1, 1, par2, 1 ); 01043 } 01044 01045 if( lgFileReadable( "0050000_50_halo__iron.bin_0.1", dum, as ) && 01046 !lgValidAsciiFile( "rauch_h-ni_halo.ascii", as ) ) 01047 { 01048 fprintf( ioQQQ, " Creating rauch_h-ni_halo.ascii....\n" ); 01049 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_halo.ascii", "_halo__iron.bin_0.1", 01050 telg2, NMODS_HNI, 1, 1, par2, 1 ); 01051 } 01052 01053 if( lgFileReadable( "0050000_50_solar_iron.bin_0.1", dum, as ) && 01054 lgFileReadable( "0050000_50_halo__iron.bin_0.1", dum, as ) && 01055 !lgValidAsciiFile( "rauch_h-ni_3d.ascii", as ) ) 01056 { 01057 fprintf( ioQQQ, " Creating rauch_h-ni_3d.ascii....\n" ); 01058 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_3d.ascii", "_solar_iron.bin_0.1", 01059 telg2, NMODS_HNI, 1, 2, par2, 1 ); 01060 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_3d.ascii", "_halo__iron.bin_0.1", 01061 telg2, NMODS_HNI, 2, 2, par2, 1 ); 01062 } 01063 01064 /* this is the hydrogen deficient PG1159 grid */ 01065 if( lgFileReadable( "0040000_5.00_33_50_02_15.bin_0.1", dum, as ) && 01066 !lgValidAsciiFile( "rauch_pg1159.ascii", as ) ) 01067 { 01068 fprintf( ioQQQ, " Creating rauch_pg1159.ascii....\n" ); 01069 lgFail = lgFail || RauchInitializeSub( "rauch_pg1159.ascii", "_33_50_02_15.bin_0.1", 01070 telg3, NMODS_PG1159, 1, 1, par2, 2 ); 01071 } 01072 01073 /* this is the pure hydrogen grid */ 01074 if( lgFileReadable( "0020000_4.00_H_00005-02000A.bin_0.1", dum, as ) && 01075 !lgValidAsciiFile( "rauch_hydr.ascii", as ) ) 01076 { 01077 fprintf( ioQQQ, " Creating rauch_hydr.ascii....\n" ); 01078 lgFail = lgFail || RauchInitializeSub( "rauch_hydr.ascii", "_H_00005-02000A.bin_0.1", 01079 telg4, NMODS_HYDR, 1, 1, par2, 2 ); 01080 } 01081 01082 /* this is the pure helium grid */ 01083 if( lgFileReadable( "0050000_5.00_He_00005-02000A.bin_0.1", dum, as ) && 01084 !lgValidAsciiFile( "rauch_helium.ascii", as ) ) 01085 { 01086 fprintf( ioQQQ, " Creating rauch_helium.ascii....\n" ); 01087 lgFail = lgFail || RauchInitializeSub( "rauch_helium.ascii", "_He_00005-02000A.bin_0.1", 01088 telg5, NMODS_HELIUM, 1, 1, par2, 2 ); 01089 } 01090 01091 /* this is the 3D grid for arbitrary H+He mixtures */ 01092 if( lgFileReadable( "0050000_5.00_H+He_1.000_0.000_00005-02000A.bin_0.1", dum, as ) && 01093 !lgValidAsciiFile( "rauch_h+he_3d.ascii", as ) ) 01094 { 01095 fprintf( ioQQQ, " Creating rauch_h+he_3d.ascii....\n" ); 01096 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_1.000_0.000_00005-02000A.bin_0.1", 01097 telg6, NMODS_HpHE, 1, 11, par3, 2 ); 01098 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.900_0.100_00005-02000A.bin_0.1", 01099 telg6, NMODS_HpHE, 2, 11, par3, 2 ); 01100 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.800_0.200_00005-02000A.bin_0.1", 01101 telg6, NMODS_HpHE, 3, 11, par3, 2 ); 01102 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.700_0.300_00005-02000A.bin_0.1", 01103 telg6, NMODS_HpHE, 4, 11, par3, 2 ); 01104 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.600_0.400_00005-02000A.bin_0.1", 01105 telg6, NMODS_HpHE, 5, 11, par3, 2 ); 01106 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.500_0.500_00005-02000A.bin_0.1", 01107 telg6, NMODS_HpHE, 6, 11, par3, 2 ); 01108 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.400_0.600_00005-02000A.bin_0.1", 01109 telg6, NMODS_HpHE, 7, 11, par3, 2 ); 01110 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.300_0.700_00005-02000A.bin_0.1", 01111 telg6, NMODS_HpHE, 8, 11, par3, 2 ); 01112 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.200_0.800_00005-02000A.bin_0.1", 01113 telg6, NMODS_HpHE, 9, 11, par3, 2 ); 01114 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.100_0.900_00005-02000A.bin_0.1", 01115 telg6, NMODS_HpHE, 10, 11, par3, 2 ); 01116 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.000_1.000_00005-02000A.bin_0.1", 01117 telg6, NMODS_HpHE, 11, 11, par3, 2 ); 01118 } 01119 01120 /* define the major absorption edges that require special attention during rebinning 01121 * 01122 * NB the frequencies should be chosen here such that they are somewhere in between 01123 * the two frequency points that straddle the edge in the atmosphere model, the 01124 * software in RebinAtmosphere will seek out the exact values of those two points 01125 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at 01126 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A). 01127 * 01128 * NB beware not to choose edges too close to one another (i.e. on the order of the 01129 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides 01130 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will 01131 * almost certainly lead to erroneous behaviour in RebinAtmosphere */ 01132 Edges[0] = 0.99946789f; 01133 Edges[1] = 1.8071406f; 01134 Edges[2] = 3.9996377f; 01135 01136 if( lgFileReadable( "rauch_h-ca_solar.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_solar.mod", pc, as ) ) 01137 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_solar.ascii", "rauch_h-ca_solar.mod",Edges,3L, pc ); 01138 if( lgFileReadable( "rauch_h-ca_halo.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_halo.mod", pc, as ) ) 01139 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_halo.ascii", "rauch_h-ca_halo.mod", Edges, 3L, pc ); 01140 if( lgFileReadable( "rauch_h-ca_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_3d.mod", pc, as ) ) 01141 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_3d.ascii", "rauch_h-ca_3d.mod", Edges, 3L, pc ); 01142 01143 if( lgFileReadable( "rauch_h-ni_solar.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_solar.mod", pc, as ) ) 01144 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_solar.ascii", "rauch_h-ni_solar.mod",Edges,3L, pc ); 01145 if( lgFileReadable( "rauch_h-ni_halo.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_halo.mod", pc, as ) ) 01146 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_halo.ascii", "rauch_h-ni_halo.mod", Edges, 3L, pc ); 01147 if( lgFileReadable( "rauch_h-ni_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_3d.mod", pc, as ) ) 01148 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_3d.ascii", "rauch_h-ni_3d.mod", Edges, 3L, pc ); 01149 01150 if( lgFileReadable( "rauch_pg1159.ascii", pc, as ) && !lgValidBinFile( "rauch_pg1159.mod", pc, as ) ) 01151 lgFail = lgFail || lgCompileAtmosphere( "rauch_pg1159.ascii", "rauch_pg1159.mod", Edges, 3L, pc ); 01152 01153 if( lgFileReadable( "rauch_hydr.ascii", pc, as ) && !lgValidBinFile( "rauch_hydr.mod", pc, as ) ) 01154 lgFail = lgFail || lgCompileAtmosphere( "rauch_hydr.ascii", "rauch_hydr.mod", Edges, 3L, pc ); 01155 01156 if( lgFileReadable( "rauch_helium.ascii", pc, as ) && !lgValidBinFile( "rauch_helium.mod", pc, as ) ) 01157 lgFail = lgFail || lgCompileAtmosphere( "rauch_helium.ascii", "rauch_helium.mod", Edges, 3L, pc ); 01158 01159 if( lgFileReadable( "rauch_h+he_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h+he_3d.mod", pc, as ) ) 01160 lgFail = lgFail || lgCompileAtmosphere( "rauch_h+he_3d.ascii", "rauch_h+he_3d.mod", Edges, 3L, pc ); 01161 return lgFail; 01162 } 01163 01164 /* RauchInterpolateHCa get one of the Rauch H-Ca model atmospheres, originally by K. Volk */ 01165 long RauchInterpolateHCa(double val[], /* val[nval] */ 01166 long *nval, 01167 long *ndim, 01168 bool lgHalo, 01169 bool lgList, 01170 double *Tlow, 01171 double *Thigh) 01172 { 01173 stellar_grid grid; 01174 01175 DEBUG_ENTRY( "RauchInterpolateHCa()" ); 01176 01177 if( *ndim == 3 ) 01178 grid.name = "rauch_h-ca_3d.mod"; 01179 else 01180 grid.name = ( lgHalo ? "rauch_h-ca_halo.mod" : "rauch_h-ca_solar.mod" ); 01181 grid.scheme = AS_DATA_OPTIONAL; 01182 /* identification of this atmosphere set, used in 01183 * the Cloudy output, *must* be 12 characters long */ 01184 grid.ident = " H-Ca Rauch"; 01185 /* the Cloudy command needed to recompile the binary model file */ 01186 grid.command = "COMPILE STARS"; 01187 01188 InitGrid( &grid, lgList ); 01189 01190 CheckVal( &grid, val, nval, ndim ); 01191 01192 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01193 01194 FreeGrid( &grid ); 01195 return rfield.nupper; 01196 } 01197 01198 /* RauchInterpolateHNi get one of the Rauch H-Ni model atmospheres */ 01199 long RauchInterpolateHNi(double val[], /* val[nval] */ 01200 long *nval, 01201 long *ndim, 01202 bool lgHalo, 01203 bool lgList, 01204 double *Tlow, 01205 double *Thigh) 01206 { 01207 stellar_grid grid; 01208 01209 DEBUG_ENTRY( "RauchInterpolateHNi()" ); 01210 01211 if( *ndim == 3 ) 01212 grid.name = "rauch_h-ni_3d.mod"; 01213 else 01214 grid.name = ( lgHalo ? "rauch_h-ni_halo.mod" : "rauch_h-ni_solar.mod" ); 01215 grid.scheme = AS_DATA_OPTIONAL; 01216 /* identification of this atmosphere set, used in 01217 * the Cloudy output, *must* be 12 characters long */ 01218 grid.ident = " H-Ni Rauch"; 01219 /* the Cloudy command needed to recompile the binary model file */ 01220 grid.command = "COMPILE STARS"; 01221 01222 InitGrid( &grid, lgList ); 01223 01224 CheckVal( &grid, val, nval, ndim ); 01225 01226 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01227 01228 FreeGrid( &grid ); 01229 return rfield.nupper; 01230 } 01231 01232 /* RauchInterpolatePG1159 get one of the Rauch PG1159 model atmospheres */ 01233 long RauchInterpolatePG1159(double val[], /* val[nval] */ 01234 long *nval, 01235 long *ndim, 01236 bool lgList, 01237 double *Tlow, 01238 double *Thigh) 01239 { 01240 stellar_grid grid; 01241 01242 DEBUG_ENTRY( "RauchInterpolatePG1159()" ); 01243 01244 grid.name = "rauch_pg1159.mod"; 01245 grid.scheme = AS_DATA_OPTIONAL; 01246 /* identification of this atmosphere set, used in 01247 * the Cloudy output, *must* be 12 characters long */ 01248 grid.ident = "PG1159 Rauch"; 01249 /* the Cloudy command needed to recompile the binary model file */ 01250 grid.command = "COMPILE STARS"; 01251 01252 InitGrid( &grid, lgList ); 01253 01254 CheckVal( &grid, val, nval, ndim ); 01255 01256 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01257 01258 FreeGrid( &grid ); 01259 return rfield.nupper; 01260 } 01261 01262 /* RauchInterpolateHydr get one of the Rauch pure hydrogen model atmospheres */ 01263 long RauchInterpolateHydr(double val[], /* val[nval] */ 01264 long *nval, 01265 long *ndim, 01266 bool lgList, 01267 double *Tlow, 01268 double *Thigh) 01269 { 01270 stellar_grid grid; 01271 01272 DEBUG_ENTRY( "RauchInterpolateHydr()" ); 01273 01274 grid.name = "rauch_hydr.mod"; 01275 grid.scheme = AS_DATA_OPTIONAL; 01276 /* identification of this atmosphere set, used in 01277 * the Cloudy output, *must* be 12 characters long */ 01278 grid.ident = " Hydr Rauch"; 01279 /* the Cloudy command needed to recompile the binary model file */ 01280 grid.command = "COMPILE STARS"; 01281 01282 InitGrid( &grid, lgList ); 01283 01284 CheckVal( &grid, val, nval, ndim ); 01285 01286 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01287 01288 FreeGrid( &grid ); 01289 return rfield.nupper; 01290 } 01291 01292 /* RauchInterpolateHelium get one of the Rauch pure helium model atmospheres */ 01293 long RauchInterpolateHelium(double val[], /* val[nval] */ 01294 long *nval, 01295 long *ndim, 01296 bool lgList, 01297 double *Tlow, 01298 double *Thigh) 01299 { 01300 stellar_grid grid; 01301 01302 DEBUG_ENTRY( "RauchInterpolateHelium()" ); 01303 01304 grid.name = "rauch_helium.mod"; 01305 grid.scheme = AS_DATA_OPTIONAL; 01306 /* identification of this atmosphere set, used in 01307 * the Cloudy output, *must* be 12 characters long */ 01308 grid.ident = "Helium Rauch"; 01309 /* the Cloudy command needed to recompile the binary model file */ 01310 grid.command = "COMPILE STARS"; 01311 01312 InitGrid( &grid, lgList ); 01313 01314 CheckVal( &grid, val, nval, ndim ); 01315 01316 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01317 01318 FreeGrid( &grid ); 01319 return rfield.nupper; 01320 } 01321 01322 /* RauchInterpolateHpHe get one of the Rauch hydrogen plus helium model atmospheres */ 01323 long RauchInterpolateHpHe(double val[], /* val[nval] */ 01324 long *nval, 01325 long *ndim, 01326 bool lgList, 01327 double *Tlow, 01328 double *Thigh) 01329 { 01330 stellar_grid grid; 01331 01332 DEBUG_ENTRY( "RauchInterpolateHpHe()" ); 01333 01334 grid.name = "rauch_h+he_3d.mod"; 01335 grid.scheme = AS_DATA_OPTIONAL; 01336 /* identification of this atmosphere set, used in 01337 * the Cloudy output, *must* be 12 characters long */ 01338 grid.ident = " H+He Rauch"; 01339 /* the Cloudy command needed to recompile the binary model file */ 01340 grid.command = "COMPILE STARS"; 01341 01342 InitGrid( &grid, lgList ); 01343 01344 CheckVal( &grid, val, nval, ndim ); 01345 01346 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01347 01348 FreeGrid( &grid ); 01349 return rfield.nupper; 01350 } 01351 01352 /* StarburstInitialize does the actual work of preparing the ascii file */ 01353 bool StarburstInitialize(const char chInName[], 01354 const char chOutName[]) 01355 { 01356 char chLine[INPUT_LINE_LENGTH]; /* used for getting input lines */ 01357 01358 bool lgHeader = true; 01359 long int i, j, nmods, ngp; 01360 01361 size_t nsb_sz = (size_t)NSB99; 01362 01363 double *wavl, *fluxes[MNTS], Age[MNTS], lwavl; 01364 01365 FILE *ioOut, /* pointer to output file we came here to create*/ 01366 *ioIn; /* pointer to input files we will read */ 01367 01368 DEBUG_ENTRY( "StarburstInitialize()" ); 01369 01370 for( i=0; i < MNTS; i++ ) 01371 fluxes[i] = NULL; 01372 01373 /* grab some space for the wavelengths and fluxes */ 01374 wavl = (double *)MALLOC( sizeof(double)*nsb_sz); 01375 01376 ioIn = open_data( chInName, "r", AS_LOCAL_ONLY ); 01377 01378 lwavl = 0.; 01379 nmods = 0; 01380 ngp = 0; 01381 01382 while( read_whole_line( chLine, INPUT_LINE_LENGTH, ioIn ) != NULL ) 01383 { 01384 if( !lgHeader ) 01385 { 01386 double cage, cwavl, cfl1; 01387 01388 /* format: age/yr wavl/Angstrom log10(flux_total) log10(flux_stellar) log10(flux_neb) */ 01389 /* we are only interested in the total flux, so we ignore the remaining numbers */ 01390 if( sscanf( chLine, " %le %le %le", &cage, &cwavl, &cfl1 ) != 3 ) 01391 { 01392 fprintf( ioQQQ, "syntax error in data of Starburst grid.\n" ); 01393 goto error; 01394 } 01395 01396 if( cwavl < lwavl ) 01397 { 01398 ++nmods; 01399 ngp = 0; 01400 01401 if( nmods >= MNTS ) 01402 { 01403 fprintf( ioQQQ, "too many time steps in Starburst grid.\n" ); 01404 fprintf( ioQQQ, "please increase MNTS and recompile.\n" ); 01405 goto error; 01406 } 01407 } 01408 01409 if( ngp == 0 ) 01410 { 01411 fluxes[nmods] = (double *)MALLOC( sizeof(double)*nsb_sz); 01412 Age[nmods] = cage; 01413 } 01414 01415 if( ngp >= (long)nsb_sz ) 01416 { 01417 /* this should only be needed when nmods == 0 */ 01418 ASSERT( nmods == 0 ); 01419 01420 nsb_sz *= 2; 01421 fluxes[0] = (double *)REALLOC(fluxes[0],sizeof(double)*nsb_sz); 01422 wavl = (double *)REALLOC(wavl,sizeof(double)*nsb_sz); 01423 } 01424 01425 if( fabs(Age[nmods]-cage) > 10.*DBL_EPSILON*Age[nmods] ) 01426 { 01427 fprintf( ioQQQ, "age error in Starburst grid.\n" ); 01428 goto error; 01429 } 01430 01431 if( nmods == 0 ) 01432 wavl[ngp] = cwavl; 01433 else 01434 { 01435 if( fabs(wavl[ngp]-cwavl) > 10.*DBL_EPSILON*wavl[ngp] ) 01436 { 01437 fprintf( ioQQQ, "wavelength error in Starburst grid.\n" ); 01438 goto error; 01439 } 01440 } 01441 01442 /* arbitrarily renormalize to flux in erg/cm^2/s/A at 1kpc */ 01443 /* constant is log10( 4*pi*(kpc/cm)^2 ) */ 01444 fluxes[nmods][ngp] = pow( 10., cfl1 - 44.077911 ); 01445 01446 lwavl = cwavl; 01447 ++ngp; 01448 } 01449 01450 if( lgHeader && strncmp( &chLine[1], "TIME [YR]", 9 ) == 0 ) 01451 lgHeader = false; 01452 } 01453 01454 if( lgHeader ) 01455 { 01456 /* this happens when the "TIME [YR]" string was not found in column 1 of the file */ 01457 fprintf( ioQQQ, "syntax error in header of Starburst grid.\n" ); 01458 goto error; 01459 } 01460 01461 ++nmods; 01462 01463 /* finished - close the unit */ 01464 fclose(ioIn); 01465 01466 ioOut = open_data( chOutName, "w", AS_LOCAL_ONLY ); 01467 01468 fprintf( ioOut, " %ld\n", VERSION_ASCII ); 01469 fprintf( ioOut, " %d\n", 1 ); 01470 fprintf( ioOut, " %d\n", 1 ); 01471 fprintf( ioOut, " Age\n" ); 01472 fprintf( ioOut, " %ld\n", nmods ); 01473 fprintf( ioOut, " %ld\n", ngp ); 01474 /* Starburst99 models give the wavelength in Angstrom */ 01475 fprintf( ioOut, " lambda\n" ); 01476 /* conversion factor to Angstrom */ 01477 fprintf( ioOut, " %.8e\n", 1. ); 01478 /* Starburst99 models give the total flux F_lambda in erg/s/A, will be renormalized at 1 kpc */ 01479 fprintf( ioOut, " F_lambda\n" ); 01480 /* this factor is irrelevant since Teff check will not be carried out */ 01481 fprintf( ioOut, " %.8e\n", 1. ); 01482 /* write out the Ages */ 01483 for( i=0; i < nmods; i++ ) 01484 { 01485 fprintf( ioOut, " %.3e", Age[i] ); 01486 if( ((i+1)%4) == 0 ) 01487 fprintf( ioOut, "\n" ); 01488 } 01489 if( (i%4) != 0 ) 01490 fprintf( ioOut, "\n" ); 01491 01492 fprintf( ioQQQ, " Writing: " ); 01493 01494 /* write out the wavelength grid */ 01495 for( j=0; j < ngp; j++ ) 01496 { 01497 fprintf( ioOut, " %.4e", wavl[j] ); 01498 /* want to have 5 numbers per line */ 01499 if( ((j+1)%5) == 0 ) 01500 fprintf( ioOut, "\n" ); 01501 } 01502 /* need to throw a newline if we did not end on an exact line */ 01503 if( (j%5) != 0 ) 01504 fprintf( ioOut, "\n" ); 01505 01506 /* print to screen to show progress */ 01507 fprintf( ioQQQ, "." ); 01508 fflush( ioQQQ ); 01509 01510 for( i=0; i < nmods; i++ ) 01511 { 01512 for( j=0; j < ngp; j++ ) 01513 { 01514 fprintf( ioOut, " %.4e", fluxes[i][j] ); 01515 /* want to have 5 numbers per line */ 01516 if( ((j+1)%5) == 0 ) 01517 fprintf( ioOut, "\n" ); 01518 } 01519 /* need to throw a newline if we did not end on an exact line */ 01520 if( (j%5) != 0 ) 01521 fprintf( ioOut, "\n" ); 01522 01523 /* print to screen to show progress */ 01524 fprintf( ioQQQ, "." ); 01525 fflush( ioQQQ ); 01526 } 01527 01528 fprintf( ioQQQ, " done.\n" ); 01529 01530 fclose(ioOut); 01531 01532 /* free the space grabbed above */ 01533 for( i=0; i < MNTS; i++ ) 01534 FREE_SAFE( fluxes[i] ); 01535 FREE_CHECK( wavl ); 01536 return false; 01537 01538 error: 01539 for( i=0; i < MNTS; i++ ) 01540 FREE_SAFE( fluxes[i] ); 01541 FREE_CHECK( wavl ); 01542 return true; 01543 } 01544 01545 /* StarburstCompile, rebin Starburst99 model output to match energy grid of code */ 01546 bool StarburstCompile(process_counter& pc) 01547 { 01548 realnum Edges[1]; 01549 01550 bool lgFail = false; 01551 01552 DEBUG_ENTRY( "StarburstCompile()" ); 01553 01554 fprintf( ioQQQ, " StarburstCompile on the job.\n" ); 01555 01556 process_counter dum; 01557 access_scheme as = AS_LOCAL_ONLY_TRY; 01558 01559 if( lgFileReadable( "starburst99.stb99", dum, as ) && !lgValidAsciiFile( "starburst99.ascii", as ) ) 01560 lgFail = lgFail || StarburstInitialize( "starburst99.stb99", "starburst99.ascii" ); 01561 if( lgFileReadable( "starburst99.ascii", pc, as ) && !lgValidBinFile( "starburst99.mod", pc, as ) ) 01562 lgFail = lgFail || lgCompileAtmosphere( "starburst99.ascii", "starburst99.mod", Edges, 0L, pc ); 01563 return lgFail; 01564 } 01565 01566 /* TlustyCompile rebin Tlusty BSTAR2006/OSTAR2002 stellar models to match energy grid of code */ 01567 int TlustyCompile(process_counter& pc) 01568 { 01569 /* these contain frequencies for the major absorption edges */ 01570 realnum Edges[1]; 01571 01572 bool lgFail = false; 01573 01574 DEBUG_ENTRY( "TlustyCompile()" ); 01575 01576 fprintf( ioQQQ, " TlustyCompile on the job.\n" ); 01577 01578 access_scheme as = AS_LOCAL_ONLY_TRY; 01579 01580 if( lgFileReadable( "bstar2006_p03.ascii", pc, as ) && !lgValidBinFile( "bstar2006_p03.mod", pc, as ) ) 01581 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_p03.ascii", "bstar2006_p03.mod", Edges, 0L, pc ); 01582 if( lgFileReadable( "bstar2006_p00.ascii", pc, as ) && !lgValidBinFile( "bstar2006_p00.mod", pc, as ) ) 01583 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_p00.ascii", "bstar2006_p00.mod", Edges, 0L, pc ); 01584 if( lgFileReadable( "bstar2006_m03.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m03.mod", pc, as ) ) 01585 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m03.ascii", "bstar2006_m03.mod", Edges, 0L, pc ); 01586 if( lgFileReadable( "bstar2006_m07.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m07.mod", pc, as ) ) 01587 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m07.ascii", "bstar2006_m07.mod", Edges, 0L, pc ); 01588 if( lgFileReadable( "bstar2006_m10.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m10.mod", pc, as ) ) 01589 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m10.ascii", "bstar2006_m10.mod", Edges, 0L, pc ); 01590 if( lgFileReadable( "bstar2006_m99.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m99.mod", pc, as ) ) 01591 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m99.ascii", "bstar2006_m99.mod", Edges, 0L, pc ); 01592 01593 if( lgFileReadable( "bstar2006_3d.ascii", pc, as ) && !lgValidBinFile( "bstar2006_3d.mod", pc, as ) ) 01594 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_3d.ascii", "bstar2006_3d.mod", Edges, 0L, pc ); 01595 01596 if( lgFileReadable( "ostar2002_p03.ascii", pc, as ) && !lgValidBinFile( "ostar2002_p03.mod", pc, as ) ) 01597 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_p03.ascii", "ostar2002_p03.mod", Edges, 0L, pc ); 01598 if( lgFileReadable( "ostar2002_p00.ascii", pc, as ) && !lgValidBinFile( "ostar2002_p00.mod", pc, as ) ) 01599 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_p00.ascii", "ostar2002_p00.mod", Edges, 0L, pc ); 01600 if( lgFileReadable( "ostar2002_m03.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m03.mod", pc, as ) ) 01601 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m03.ascii", "ostar2002_m03.mod", Edges, 0L, pc ); 01602 if( lgFileReadable( "ostar2002_m07.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m07.mod", pc, as ) ) 01603 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m07.ascii", "ostar2002_m07.mod", Edges, 0L, pc ); 01604 if( lgFileReadable( "ostar2002_m10.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m10.mod", pc, as ) ) 01605 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m10.ascii", "ostar2002_m10.mod", Edges, 0L, pc ); 01606 if( lgFileReadable( "ostar2002_m15.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m15.mod", pc, as ) ) 01607 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m15.ascii", "ostar2002_m15.mod", Edges, 0L, pc ); 01608 if( lgFileReadable( "ostar2002_m17.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m17.mod", pc, as ) ) 01609 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m17.ascii", "ostar2002_m17.mod", Edges, 0L, pc ); 01610 if( lgFileReadable( "ostar2002_m20.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m20.mod", pc, as ) ) 01611 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m20.ascii", "ostar2002_m20.mod", Edges, 0L, pc ); 01612 if( lgFileReadable( "ostar2002_m30.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m30.mod", pc, as ) ) 01613 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m30.ascii", "ostar2002_m30.mod", Edges, 0L, pc ); 01614 if( lgFileReadable( "ostar2002_m99.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m99.mod", pc, as ) ) 01615 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m99.ascii", "ostar2002_m99.mod", Edges, 0L, pc ); 01616 01617 if( lgFileReadable( "ostar2002_3d.ascii", pc, as ) && !lgValidBinFile( "ostar2002_3d.mod", pc, as ) ) 01618 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_3d.ascii", "ostar2002_3d.mod", Edges, 0L, pc ); 01619 return lgFail; 01620 } 01621 01622 /* TlustyInterpolate get one of the Tlusty BSTAR2006/OSTAR2002 model atmospheres */ 01623 long TlustyInterpolate(double val[], /* val[nval] */ 01624 long *nval, 01625 long *ndim, 01626 tl_grid tlg, 01627 const char *chMetalicity, 01628 bool lgList, 01629 double *Tlow, 01630 double *Thigh) 01631 { 01632 char chIdent[13]; 01633 stellar_grid grid; 01634 01635 DEBUG_ENTRY( "TlustyInterpolate()" ); 01636 01637 if( tlg == TL_BSTAR ) 01638 grid.name = "bstar2006_"; 01639 else if( tlg == TL_OSTAR ) 01640 grid.name = "ostar2002_"; 01641 else 01642 TotalInsanity(); 01643 if( *ndim == 3 ) 01644 grid.name += "3d"; 01645 else 01646 grid.name += chMetalicity; 01647 grid.name += ".mod"; 01648 grid.scheme = AS_DATA_OPTIONAL; 01649 /* identification of this atmosphere set, used in 01650 * the Cloudy output, *must* be 12 characters long */ 01651 if( *ndim == 3 ) 01652 { 01653 strcpy( chIdent, "3-dim" ); 01654 } 01655 else 01656 { 01657 strcpy( chIdent, "Z " ); 01658 strcat( chIdent, chMetalicity ); 01659 } 01660 if( tlg == TL_BSTAR ) 01661 strcat( chIdent, " Bstr06" ); 01662 else if( tlg == TL_OSTAR ) 01663 strcat( chIdent, " Ostr02" ); 01664 else 01665 TotalInsanity(); 01666 grid.ident = chIdent; 01667 /* the Cloudy command needed to recompile the binary model file */ 01668 grid.command = "COMPILE STARS"; 01669 01670 InitGrid( &grid, lgList ); 01671 01672 CheckVal( &grid, val, nval, ndim ); 01673 01674 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01675 01676 FreeGrid( &grid ); 01677 return rfield.nupper; 01678 } 01679 01680 /* WernerCompile rebin Werner stellar models to match energy grid of code */ 01681 /* >>chng 05 nov 16, added return value to indicate success (0) or failure (1) */ 01682 int WernerCompile(process_counter& pc) 01683 { 01684 /* these contain frequencies for the major absorption edges */ 01685 realnum Edges[3]; 01686 01687 bool lgFail = false; 01688 01689 DEBUG_ENTRY( "WernerCompile()" ); 01690 01691 fprintf( ioQQQ, " WernerCompile on the job.\n" ); 01692 01693 /* define the major absorption edges that require special attention during rebinning 01694 * 01695 * NB the frequencies should be chosen here such that they are somewhere in between 01696 * the two frequency points that straddle the edge in the atmosphere model, the 01697 * software in RebinAtmosphere will seek out the exact values of those two points 01698 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at 01699 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A). 01700 * 01701 * NB beware not to choose edges too close to one another (i.e. on the order of the 01702 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides 01703 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will 01704 * almost certainly lead to erroneous behaviour in RebinAtmosphere */ 01705 Edges[0] = 0.99946789f; 01706 Edges[1] = 1.8071406f; 01707 Edges[2] = 3.9996377f; 01708 01709 /* The "kwerner.ascii" file is a modified ascii dump of the Klaus Werner 01710 * stellar model files which he gave to me in 1992. The first set of values 01711 * is the frequency grid (in Ryd) followed by the atmosphere models in order 01712 * of increasing temperature and log(g). The following comments are already 01713 * incorporated in the modified kwerner.ascii file that is supplied with Cloudy. 01714 * 01715 * >>chng 00 oct 18, The frequency grid was slightly tweaked compared to the 01716 * original values supplied by Klaus Werner to make it monotonically increasing; 01717 * this is due to there being fluxes above and below certain wavelengths where 01718 * the opacity changes (i.e. the Lyman and Balmer limits for example) which are 01719 * assigned the same wavelength in the original Klaus Werner files. PvH 01720 * 01721 * >>chng 00 oct 20, StarEner[172] is out of sequence. As per the Klaus Werner comment, 01722 * it should be omitted. The energy grid is very dense in this region and was most likely 01723 * intended to sample an absorption line which was not included in this particular grid. 01724 * StarFlux[172] is therefore always equal to the flux in neighbouring points (at least 01725 * those with slightly smaller energies). It is therefore safe to ignore this point. PvH 01726 * 01727 * >>chng 00 oct 20, As per the same comment, StarFlux[172] is also deleted. PvH */ 01728 01729 access_scheme as = AS_LOCAL_ONLY_TRY; 01730 01731 if( lgFileReadable( "kwerner.ascii", pc, as ) && !lgValidBinFile( "kwerner.mod", pc, as ) ) 01732 lgFail = lgCompileAtmosphere( "kwerner.ascii", "kwerner.mod", Edges, 3L, pc ); 01733 return lgFail; 01734 } 01735 01736 /* WernerInterpolate read in and interpolate on Werner grid of PN atmospheres, originally by K Volk */ 01737 long WernerInterpolate(double val[], /* val[nval] */ 01738 long *nval, 01739 long *ndim, 01740 bool lgList, 01741 double *Tlow, 01742 double *Thigh) 01743 { 01744 stellar_grid grid; 01745 01746 DEBUG_ENTRY( "WernerInterpolate()" ); 01747 01748 /* This subroutine was added (28 dec 1992) to read from the set of 01749 * hot white dwarf model atmospheres from Klaus Werner at Kiel. The 01750 * values are read in (energy in Rydberg units, f_nu in cgs units) 01751 * for any of the 20 models. Each model had 513 points before rebinning. 01752 * The Rayleigh-Jeans tail was extrapolated. */ 01753 01754 grid.name = "kwerner.mod"; 01755 grid.scheme = AS_DATA_OPTIONAL; 01756 /* identification of this atmosphere set, used in 01757 * the Cloudy output, *must* be 12 characters long */ 01758 grid.ident = "Klaus Werner"; 01759 /* the Cloudy command needed to recompile the binary model file */ 01760 grid.command = "COMPILE STARS"; 01761 01762 InitGrid( &grid, lgList ); 01763 01764 CheckVal( &grid, val, nval, ndim ); 01765 01766 /* Note on the interpolation: 26 October 2000 (Peter van Hoof) 01767 * 01768 * I computed the effective temperature for a random sample of interpolated 01769 * atmospheres by integrating the flux as shown above and compared the results 01770 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC. 01771 * 01772 * I found that the average discrepancy was: 01773 * 01774 * DELTA = -0.71% +/- 0.71% (sample size 5000) 01775 * 01776 * The most extreme discrepancies were 01777 * -4.37% <= DELTA <= 0.24% 01778 * 01779 * The most negative discrepancies were for Teff = 95 kK, log(g) = 5 01780 * The most positive discrepancies were for Teff = 160 kK, log(g) = 8 01781 * 01782 * Since Cloudy checks the scaling elsewhere there is no need to re-scale 01783 * things here, but this inaccuracy should be kept in mind since it could 01784 * indicate problems with the flux distribution */ 01785 01786 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01787 01788 FreeGrid( &grid ); 01789 return rfield.nupper; 01790 } 01791 01792 /* WMBASICCompile rebin WMBASIC stellar models to match energy grid of code */ 01793 int WMBASICCompile(process_counter& pc) 01794 { 01795 /* these contain frequencies for the major absorption edges */ 01796 realnum Edges[3]; 01797 01798 bool lgFail = false; 01799 01800 DEBUG_ENTRY( "WMBASICCompile()" ); 01801 01802 fprintf( ioQQQ, " WMBASICCompile on the job.\n" ); 01803 01804 /* define the major absorption edges that require special attention during rebinning */ 01805 Edges[0] = 0.99946789f; 01806 Edges[1] = 1.8071406f; 01807 Edges[2] = 3.9996377f; 01808 01809 access_scheme as = AS_LOCAL_ONLY_TRY; 01810 01811 if( lgFileReadable( "wmbasic.ascii", pc, as ) && !lgValidBinFile( "wmbasic.mod", pc, as ) ) 01812 lgFail = lgCompileAtmosphere( "wmbasic.ascii", "wmbasic.mod", Edges, 3L, pc ); 01813 return lgFail; 01814 } 01815 01816 /* WMBASICInterpolate read in and interpolate on WMBASIC grid of hot star atmospheres */ 01817 long WMBASICInterpolate(double val[], /* val[nval] */ 01818 long *nval, 01819 long *ndim, 01820 bool lgList, 01821 double *Tlow, 01822 double *Thigh) 01823 { 01824 stellar_grid grid; 01825 01826 DEBUG_ENTRY( "WMBASICInterpolate()" ); 01827 01828 grid.name = "wmbasic.mod"; 01829 grid.scheme = AS_DATA_OPTIONAL; 01830 /* identification of this atmosphere set, used in 01831 * the Cloudy output, *must* be 12 characters long */ 01832 grid.ident = " WMBASIC"; 01833 /* the Cloudy command needed to recompile the binary model file */ 01834 grid.command = "COMPILE STARS"; 01835 01836 InitGrid( &grid, lgList ); 01837 01838 CheckVal( &grid, val, nval, ndim ); 01839 01840 InterpolateRectGrid( &grid, val, Tlow, Thigh ); 01841 01842 FreeGrid( &grid ); 01843 return rfield.nupper; 01844 } 01845 01846 /* CompileAtmosphereCoStar rebin costar stellar atmospheres to match cloudy energy grid, 01847 * called by the compile stars command */ 01848 STATIC bool lgCompileAtmosphereCoStar(const char chFNameIn[], 01849 const char chFNameOut[], 01850 const realnum Edges[], /* Edges[nedges] */ 01851 long nedges, 01852 process_counter& pc) 01853 { 01854 char chLine[INPUT_LINE_LENGTH]; 01855 char names[MDIM][MNAM+1]; 01856 int32 val[7]; 01857 uint32 uval[2]; 01858 long int i, j, nskip, nModels, nWL; 01859 01860 /* these will be malloced into large work arrays*/ 01861 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL; 01862 /* this will hold all the model parameters */ 01863 mpp *telg = NULL; 01864 01865 FILE *ioIN; /* used for input */ 01866 FILE *ioOUT; /* used for output */ 01867 01868 DEBUG_ENTRY( "CompileAtmosphereCoStar()" ); 01869 01870 /* This is a program to re-bin the costar stellar model spectra to match the 01871 * Cloudy grid. For short wavelengths I will use a power law extrapolation 01872 * of the model values (which should be falling rapidly) if needed. At long 01873 * wavelengths I will assume Rayleigh-Jeans from the last stellar model point 01874 * to extrapolate to 1 cm wavelength. */ 01875 01876 /* This version uses power-law interpolation between the points of the stellar model. */ 01877 01878 /* read the original data file obtained off the web, 01879 * open as read only */ 01880 try 01881 { 01882 ioIN = open_data( chFNameIn, "r", AS_LOCAL_ONLY ); 01883 } 01884 catch( cloudy_exit ) 01885 { 01886 goto error; 01887 } 01888 fprintf( ioQQQ, " CompileAtmosphereCoStar got %s.\n", chFNameIn ); 01889 01890 /* get first line and see how many more to skip */ 01891 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL ) 01892 { 01893 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading nskip.\n" ); 01894 goto error; 01895 } 01896 sscanf( chLine, "%li", &nskip ); 01897 01898 /* now skip the header information */ 01899 for( i=0; i < nskip; ++i ) 01900 { 01901 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL ) 01902 { 01903 fprintf( ioQQQ, " CompileAtmosphereCoStar fails skipping header.\n" ); 01904 goto error; 01905 } 01906 } 01907 01908 /* now get number of models and number of wavelengths */ 01909 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL ) 01910 { 01911 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading nModels, nWL.\n" ); 01912 goto error; 01913 } 01914 sscanf( chLine, "%li%li", &nModels, &nWL ); 01915 01916 if( nModels <= 0 || nWL <= 0 ) 01917 { 01918 fprintf( ioQQQ, " CompileAtmosphereCoStar scanned off impossible values for nModels=%li or nWL=%li\n", 01919 nModels, nWL ); 01920 goto error; 01921 } 01922 01923 /* allocate space for the stellar parameters */ 01924 telg = (mpp *)CALLOC( (size_t)nModels, sizeof(mpp) ); 01925 01926 /* get all model parameters for the atmospheres */ 01927 for( i=0; i < nModels; ++i ) 01928 { 01929 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL ) 01930 { 01931 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading model parameters.\n" ); 01932 goto error; 01933 } 01934 /* first letter on line is indicator of grid */ 01935 telg[i].chGrid = chLine[0]; 01936 /* get the model id number */ 01937 sscanf( chLine+1, "%i", &telg[i].modid ); 01938 /* get the temperature */ 01939 sscanf( chLine+23, "%lg", &telg[i].par[0] ); 01940 /* get the surface gravity */ 01941 sscanf( chLine+31, "%lg", &telg[i].par[1] ); 01942 /* get the ZAMS mass */ 01943 sscanf( chLine+7, "%lg", &telg[i].par[2] ); 01944 /* get the model age */ 01945 sscanf( chLine+15, "%lg", &telg[i].par[3] ); 01946 01947 /* the code in parse_table.cpp implicitly depends on this! */ 01948 ASSERT( telg[i].par[2] > 10. ); 01949 ASSERT( telg[i].par[3] > 10. ); 01950 01951 /* convert ZAMS masses to logarithms */ 01952 telg[i].par[2] = log10(telg[i].par[2]); 01953 } 01954 01955 /* this will be the file we create, that will be read to compute models, 01956 * open to write binary */ 01957 try 01958 { 01959 ioOUT = open_data( chFNameOut, "wb", AS_LOCAL_ONLY ); 01960 } 01961 catch( cloudy_exit ) 01962 { 01963 goto error; 01964 } 01965 01966 /* >>chng 05 dec 01, use int32 instead of long so that binary file is compatible for 32/64-bit code */ 01967 /* >>chng 05 dec 01, add several parameters to binary file, change sequence, PvH */ 01968 /* >>chng 06 jun 10, add several parameters to binary file, change sequence, PvH */ 01969 val[0] = (int32)VERSION_BIN; 01970 val[1] = (int32)MDIM; 01971 val[2] = (int32)MNAM; 01972 val[3] = (int32)2; 01973 val[4] = (int32)4; 01974 val[5] = (int32)nModels; 01975 val[6] = (int32)rfield.nupper; 01976 uval[0] = sizeof(val) + sizeof(uval) + sizeof(names) + nModels*sizeof(mpp); /* nOffset */ 01977 uval[1] = rfield.nupper*sizeof(realnum); /* nBlocksize */ 01978 01979 strncpy( names[0], "Teff\0\0", MNAM+1 ); 01980 strncpy( names[1], "log(g)", MNAM+1 ); 01981 strncpy( names[2], "log(M)", MNAM+1 ); 01982 strncpy( names[3], "Age\0\0\0", MNAM+1 ); 01983 01984 if( fwrite( val, sizeof(val), 1, ioOUT ) != 1 || 01985 fwrite( uval, sizeof(uval), 1, ioOUT ) != 1 || 01986 fwrite( names, sizeof(names), 1, ioOUT ) != 1 || 01987 /* write out the array of {Teff,log(g)} pairs */ 01988 fwrite( telg, sizeof(mpp), (size_t)nModels, ioOUT ) != (size_t)nModels || 01989 /* write out the cloudy energy grid for later sanity checks */ 01990 fwrite( rfield.AnuOrg, (size_t)uval[1], 1, ioOUT ) != 1 ) 01991 { 01992 fprintf( ioQQQ, " CompileAtmosphereCoStar failed writing header of output file.\n" ); 01993 goto error; 01994 } 01995 01996 /* MALLOC some workspace */ 01997 StarEner = (realnum *)MALLOC( sizeof(realnum)*nWL ); 01998 StarFlux = (realnum *)MALLOC( sizeof(realnum)*nWL ); 01999 CloudyFlux = (realnum *)MALLOC( (size_t)uval[1] ); 02000 02001 fprintf( ioQQQ, " Compiling: " ); 02002 02003 /* get the star data */ 02004 for( i=0; i < nModels; ++i ) 02005 { 02006 /* get number to skip */ 02007 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL ) 02008 { 02009 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading the skip to next spectrum.\n" ); 02010 goto error; 02011 } 02012 sscanf( chLine, "%li", &nskip ); 02013 02014 for( j=0; j < nskip; ++j ) 02015 { 02016 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL ) 02017 { 02018 fprintf( ioQQQ, " CompileAtmosphereCoStar fails doing the skip.\n" ); 02019 goto error; 02020 } 02021 } 02022 02023 /* now read in the wavelength and flux for this star, read in 02024 * backwards since we want to be in increasing energy order rather 02025 * than wavelength */ 02026 for( j=nWL-1; j >= 0; --j ) 02027 { 02028 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL ) 02029 { 02030 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading the spectral data.\n" ); 02031 goto error; 02032 } 02033 double help1, help2; 02034 sscanf( chLine, "%lg %lg", &help1, &help2 ); 02035 02036 /* continuum flux was log, convert to linear, also do 02037 * conversion from "astrophysical" flux to F_nu in cgs units */ 02038 StarFlux[j] = (realnum)(PI*pow(10.,help2)); 02039 /* StarEner was in Angstroms, convert to Ryd */ 02040 StarEner[j] = (realnum)(RYDLAM/help1); 02041 02042 /* sanity check */ 02043 if( j < nWL-1 ) 02044 ASSERT( StarEner[j] < StarEner[j+1] ); 02045 } 02046 02047 /* this will do the heavy lifting, and define arrays used below, 02048 * NB the lowest energy point in these grids appears to be bogus. 02049 * tell rebin about nWL-1 */ 02050 RebinAtmosphere(nWL-1, StarEner+1, StarFlux+1, CloudyFlux, nedges, Edges ); 02051 02052 /* write the continuum out as a binary file */ 02053 if( fwrite( CloudyFlux, (size_t)uval[1], 1, ioOUT ) != 1 ) 02054 { 02055 fprintf( ioQQQ, " CompileAtmosphereCoStar failed writing star flux.\n" ); 02056 goto error; 02057 } 02058 02059 fprintf( ioQQQ, "." ); 02060 fflush( ioQQQ ); 02061 } 02062 02063 fprintf( ioQQQ, " done.\n" ); 02064 02065 fclose( ioIN ); 02066 fclose( ioOUT ); 02067 02068 FREE_CHECK( telg ); 02069 FREE_CHECK( StarEner ); 02070 FREE_CHECK( StarFlux ); 02071 FREE_CHECK( CloudyFlux ); 02072 02073 fprintf( ioQQQ, "\n CompileAtmosphereCoStar completed ok.\n\n" ); 02074 ++pc.nOK; 02075 return false; 02076 02077 error: 02078 FREE_SAFE( telg ); 02079 FREE_SAFE( StarEner ); 02080 FREE_SAFE( StarFlux ); 02081 FREE_SAFE( CloudyFlux ); 02082 ++pc.nFail; 02083 return true; 02084 } 02085 02086 /* InterpolateGridCoStar read in and interpolate on costar grid of windy O atmospheres */ 02087 STATIC void InterpolateGridCoStar(const stellar_grid *grid, /* struct with all the grid parameters */ 02088 const double val[], /* val[0]: Teff for imode = 1,2; M_ZAMS for imode = 3; 02089 * age for imode = 4 */ 02090 /* val[1]: nmodid for imode = 1; log(g) for imode = 2; 02091 * age for imode = 3; M_ZAMS for imode = 4 */ 02092 double *val0_lo, 02093 double *val0_hi) 02094 { 02095 long i, j, k, nmodid, off, ptr; 02096 long *indloTr, *indhiTr, useTr[2]; 02097 long indlo[2], indhi[2], index[2]; 02098 realnum *ValTr; 02099 double lval[2], aval[4]; 02100 02101 DEBUG_ENTRY( "InterpolateGridCoStar()" ); 02102 02103 switch( grid->imode ) 02104 { 02105 case IM_COSTAR_TEFF_MODID: 02106 case IM_COSTAR_TEFF_LOGG: 02107 lval[0] = val[0]; 02108 lval[1] = val[1]; 02109 off = 0; 02110 break; 02111 case IM_COSTAR_MZAMS_AGE: 02112 lval[0] = log10(val[0]); /* use log10(M_ZAMS) internally */ 02113 lval[1] = val[1]; 02114 off = 2; 02115 break; 02116 case IM_COSTAR_AGE_MZAMS: 02117 /* swap parameters, hence mimic IM_COSTAR_MZAMS_AGE */ 02118 lval[0] = log10(val[1]); /* use log10(M_ZAMS) internally */ 02119 lval[1] = val[0]; 02120 off = 2; 02121 break; 02122 default: 02123 fprintf( ioQQQ, " InterpolateGridCoStar called with insane value for imode: %d.\n", grid->imode ); 02124 cdEXIT(EXIT_FAILURE); 02125 } 02126 02127 nmodid = (long)(lval[1]+0.5); 02128 02129 ASSERT( rfield.lgContMalloc[rfield.nspec] ); 02130 02131 /* read in the saved cloudy energy scale so we can confirm this is a good image */ 02132 GetModel( grid, -1, rfield.tNuRyd[rfield.nspec], lgSILENT, lgLINEAR ); 02133 02134 # if DEBUGPRT 02135 /* check whether the models in the grid have the correct effective temperature */ 02136 ValidateGrid( grid, 0.005 ); 02137 # endif 02138 02139 /* now allocate some temp workspace */ 02140 ValTr = (realnum *)MALLOC( sizeof(realnum)*grid->nTracks ); 02141 indloTr = (long *)MALLOC( sizeof(long)*grid->nTracks ); 02142 indhiTr = (long *)MALLOC( sizeof(long)*grid->nTracks ); 02143 02144 /* first do horizontal search, i.e. search along individual tracks */ 02145 for( j=0; j < grid->nTracks; j++ ) 02146 { 02147 if( grid->imode == IM_COSTAR_TEFF_MODID ) 02148 { 02149 if( grid->trackLen[j] >= nmodid ) { 02150 index[0] = nmodid - 1; 02151 index[1] = j; 02152 ptr = grid->jval[JIndex(grid,index)]; 02153 indloTr[j] = ptr; 02154 indhiTr[j] = ptr; 02155 ValTr[j] = (realnum)grid->telg[ptr].par[off]; 02156 } 02157 else 02158 { 02159 indloTr[j] = -2; 02160 indhiTr[j] = -2; 02161 ValTr[j] = -FLT_MAX; 02162 } 02163 } 02164 else 02165 { 02166 FindHCoStar( grid, j, lval[1], off, ValTr, indloTr, indhiTr ); 02167 } 02168 } 02169 02170 # if DEBUGPRT 02171 for( j=0; j < grid->nTracks; j++ ) 02172 { 02173 if( indloTr[j] >= 0 ) 02174 printf( "track %c: models %c%d, %c%d, val %g\n", 02175 (char)('A'+j), grid->telg[indloTr[j]].chGrid, grid->telg[indloTr[j]].modid, 02176 grid->telg[indhiTr[j]].chGrid, grid->telg[indhiTr[j]].modid, ValTr[j]); 02177 } 02178 # endif 02179 02180 /* now do vertical search, i.e. interpolate between tracks */ 02181 FindVCoStar( grid, lval[0], ValTr, useTr ); 02182 02183 /* This should only happen when InterpolateGridCoStar is called in non-optimizing mode, 02184 * when optimizing InterpolateGridCoStar should report back to optimize_func()... 02185 * The fact that FindVCoStar allows interpolation between non-adjoining tracks 02186 * should guarantee that this will not happen. */ 02187 if( useTr[0] < 0 ) 02188 { 02189 fprintf( ioQQQ, " The parameters for the requested CoStar model are out of range.\n" ); 02190 cdEXIT(EXIT_FAILURE); 02191 } 02192 02193 ASSERT( useTr[0] >= 0 && useTr[0] < grid->nTracks ); 02194 ASSERT( useTr[1] >= 0 && useTr[1] < grid->nTracks ); 02195 ASSERT( indloTr[useTr[0]] >= 0 && indloTr[useTr[0]] < (int)grid->nmods ); 02196 ASSERT( indhiTr[useTr[0]] >= 0 && indhiTr[useTr[0]] < (int)grid->nmods ); 02197 ASSERT( indloTr[useTr[1]] >= 0 && indloTr[useTr[1]] < (int)grid->nmods ); 02198 ASSERT( indhiTr[useTr[1]] >= 0 && indhiTr[useTr[1]] < (int)grid->nmods ); 02199 02200 # if DEBUGPRT 02201 printf( "interpolate between tracks %c and %c\n", (char)('A'+useTr[0]), (char)('A'+useTr[1]) ); 02202 # endif 02203 02204 indlo[0] = indloTr[useTr[0]]; 02205 indhi[0] = indhiTr[useTr[0]]; 02206 indlo[1] = indloTr[useTr[1]]; 02207 indhi[1] = indhiTr[useTr[1]]; 02208 02209 InterpolateModelCoStar( grid, lval, aval, indlo, indhi, index, 0, off, rfield.tslop[rfield.nspec] ); 02210 02211 for( i=0; i < rfield.nupper; i++ ) 02212 { 02213 rfield.tslop[rfield.nspec][i] = (realnum)pow((realnum)10.f,rfield.tslop[rfield.nspec][i]); 02214 if( rfield.tslop[rfield.nspec][i] < 1e-37 ) 02215 rfield.tslop[rfield.nspec][i] = 0.; 02216 } 02217 02218 if( false ) 02219 { 02220 FILE *ioBUG = fopen( "interpolated.txt", "w" ); 02221 for( k=0; k < rfield.nupper; ++k ) 02222 fprintf( ioBUG, "%e %e\n", rfield.tNuRyd[rfield.nspec][k], rfield.tslop[rfield.nspec][k] ); 02223 fclose( ioBUG ); 02224 } 02225 02226 /* sanity check: see whether this model has the correct effective temperature */ 02227 if( ! lgValidModel( rfield.tNuRyd[rfield.nspec], rfield.tslop[rfield.nspec], aval[0], 0.05 ) ) 02228 TotalInsanity(); 02229 02230 /* set limits for optimizer */ 02231 SetLimits( grid, lval[0], NULL, NULL, useTr, ValTr, val0_lo, val0_hi ); 02232 02233 /* now write some final info */ 02234 if( called.lgTalk ) 02235 { 02236 fprintf( ioQQQ, " * c<< FINAL: T_eff = %7.1f, ", aval[0] ); 02237 fprintf( ioQQQ, "log(g) = %4.2f, M(ZAMS) = %5.1f, age = ", aval[1], pow(10.,aval[2]) ); 02238 fprintf( ioQQQ, PrintEfmt("%8.2e",aval[3]) ); 02239 fprintf( ioQQQ, " >>> *\n" ); 02240 } 02241 02242 FREE_CHECK( indhiTr ); 02243 FREE_CHECK( indloTr ); 02244 FREE_CHECK( ValTr ); 02245 return; 02246 } 02247 02248 /* find which models to use for interpolation along a given evolutionary track */ 02249 STATIC void FindHCoStar(const stellar_grid *grid, 02250 long track, 02251 double par2, /* requested log(g) or age */ 02252 long off, /* determines which parameter to match 0 -> log(g), 2 -> age */ 02253 realnum *ValTr, /* ValTr[track]: Teff/log(M) value for interpolated model along track */ 02254 long *indloTr, /* indloTr[track]: model number for first model used in interpolation */ 02255 long *indhiTr) /* indhiTr[track]: model number for second model used in interpolation */ 02256 { 02257 long index[2], j, mod1, mod2; 02258 02259 DEBUG_ENTRY( "FindHCoStar()" ); 02260 02261 indloTr[track] = -2; 02262 indhiTr[track] = -2; 02263 ValTr[track] = -FLT_MAX; 02264 02265 index[1] = track; 02266 02267 for( j=0; j < grid->trackLen[track]; j++ ) 02268 { 02269 index[0] = j; 02270 mod1 = grid->jval[JIndex(grid,index)]; 02271 02272 /* do we have an exact match ? */ 02273 if( fabs(par2-grid->telg[mod1].par[off+1]) <= 10.*FLT_EPSILON*fabs(grid->telg[mod1].par[off+1]) ) 02274 { 02275 indloTr[track] = mod1; 02276 indhiTr[track] = mod1; 02277 ValTr[track] = (realnum)grid->telg[mod1].par[off]; 02278 return; 02279 } 02280 } 02281 02282 for( j=0; j < grid->trackLen[track]-1; j++ ) 02283 { 02284 index[0] = j; 02285 mod1 = grid->jval[JIndex(grid,index)]; 02286 index[0] = j+1; 02287 mod2 = grid->jval[JIndex(grid,index)]; 02288 02289 /* do we interpolate ? */ 02290 if( (par2 - grid->telg[mod1].par[off+1])*(par2 - grid->telg[mod2].par[off+1]) < 0. ) 02291 { 02292 double frac; 02293 02294 indloTr[track] = mod1; 02295 indhiTr[track] = mod2; 02296 frac = (par2 - grid->telg[mod2].par[off+1])/ 02297 (grid->telg[mod1].par[off+1] - grid->telg[mod2].par[off+1]); 02298 ValTr[track] = (realnum)(frac*grid->telg[mod1].par[off] + 02299 (1.-frac)*grid->telg[mod2].par[off] ); 02300 break; 02301 } 02302 } 02303 return; 02304 } 02305 02306 /* find which tracks to use for interpolation in between tracks */ 02307 STATIC void FindVCoStar(const stellar_grid *grid, 02308 double par1, /* requested Teff or ZAMS mass */ 02309 realnum *ValTr, /* internal workspace */ 02310 long useTr[]) /* useTr[0]: track number for first track to be used in interpolation 02311 * (i.e., 0 means 'A', etc.) 02312 * useTr[1]: track number for second track to be used in interpolation 02313 * NOTE: FindVCoStar raises a flag when interpolating between non-adjoining 02314 * tracks, i.e. when (useTr[1]-useTr[0]) > 1 */ 02315 { 02316 long j; 02317 02318 DEBUG_ENTRY( "FindVCoStar()" ); 02319 02320 useTr[0] = -1; 02321 useTr[1] = -1; 02322 02323 for( j=0; j < grid->nTracks; j++ ) 02324 { 02325 /* do we have an exact match ? */ 02326 if( ValTr[j] != -FLT_MAX && fabs(par1-(double)ValTr[j]) <= 10.*FLT_EPSILON*fabs(ValTr[j]) ) 02327 { 02328 useTr[0] = j; 02329 useTr[1] = j; 02330 break; 02331 } 02332 } 02333 02334 if( useTr[0] >= 0 ) 02335 { 02336 return; 02337 } 02338 02339 for( j=0; j < grid->nTracks-1; j++ ) 02340 { 02341 if( ValTr[j] != -FLT_MAX ) 02342 { 02343 long int i,j2; 02344 02345 /* find next valid track */ 02346 j2 = 0; 02347 for( i = j+1; i < grid->nTracks; i++ ) 02348 { 02349 if( ValTr[i] != -FLT_MAX ) 02350 { 02351 j2 = i; 02352 break; 02353 } 02354 } 02355 02356 /* do we interpolate ? */ 02357 if( j2 > 0 && ((realnum)par1-ValTr[j])*((realnum)par1-ValTr[j2]) < 0.f ) 02358 { 02359 useTr[0] = j; 02360 useTr[1] = j2; 02361 break; 02362 } 02363 } 02364 } 02365 02366 /* raise caution when we interpolate between non-adjoining tracks */ 02367 continuum.lgCoStarInterpolationCaution = ( useTr[1]-useTr[0] > 1 ); 02368 return; 02369 } 02370 02371 /* Make a listing of all the models in the CoStar grid */ 02372 STATIC void CoStarListModels(const stellar_grid *grid) 02373 { 02374 long index[2], maxlen, n; 02375 02376 DEBUG_ENTRY( "CoStarListModels()" ); 02377 02378 maxlen = 0; 02379 for( n=0; n < grid->nTracks; n++ ) 02380 maxlen = MAX2( maxlen, grid->trackLen[n] ); 02381 02382 fprintf( ioQQQ, "\n" ); 02383 fprintf( ioQQQ, " Track\\Index |" ); 02384 for( n = 0; n < maxlen; n++ ) 02385 fprintf( ioQQQ, " %5ld ", n+1 ); 02386 fprintf( ioQQQ, "\n" ); 02387 fprintf( ioQQQ, "--------------|" ); 02388 for( n = 0; n < maxlen; n++ ) 02389 fprintf( ioQQQ, "----------------" ); 02390 fprintf( ioQQQ, "\n" ); 02391 02392 for( index[1]=0; index[1] < grid->nTracks; ++index[1] ) 02393 { 02394 long ptr; 02395 double Teff, alogg, Mass; 02396 02397 fprintf( ioQQQ, " %c", (char)('A'+index[1]) ); 02398 index[0] = 0; 02399 ptr = grid->jval[JIndex(grid,index)]; 02400 Mass = pow(10.,grid->telg[ptr].par[2]); 02401 fprintf( ioQQQ, " (%3.0f Msol) |", Mass ); 02402 02403 for( index[0]=0; index[0] < grid->trackLen[index[1]]; ++index[0] ) 02404 { 02405 ptr = grid->jval[JIndex(grid,index)]; 02406 Teff = grid->telg[ptr].par[0]; 02407 alogg = grid->telg[ptr].par[1]; 02408 fprintf( ioQQQ, " (%6.1f,%4.2f)", Teff, alogg ); 02409 } 02410 fprintf( ioQQQ, "\n" ); 02411 } 02412 return; 02413 } 02414 02415 /* RauchInitializeSub does the actual work of preparing the ascii file */ 02416 STATIC int RauchInitializeSub(const char chFName[], 02417 const char chSuff[], 02418 const mpp telg[], 02419 long nmods, 02420 long n, 02421 long ngrids, 02422 const double par2[], /* par2[ngrids] */ 02423 int format) 02424 { 02425 char chLine[INPUT_LINE_LENGTH]; /* used for getting input lines */ 02426 02427 FILE *ioOut, /* pointer to output file we came here to create*/ 02428 *ioIn; /* pointer to input files we will read */ 02429 02430 long int i, j; 02431 02432 double *wavl, *fluxes; 02433 02434 DEBUG_ENTRY( "RauchInitializeSub()" ); 02435 02436 /* grab some space for the wavelengths and fluxes */ 02437 wavl = (double *)MALLOC( sizeof(double)*NRAUCH); 02438 fluxes = (double *)MALLOC( sizeof(double)*NRAUCH); 02439 02440 try 02441 { 02442 if( n == 1 ) 02443 ioOut = open_data( chFName, "w", AS_LOCAL_ONLY ); 02444 else 02445 ioOut = open_data( chFName, "a", AS_LOCAL_ONLY ); 02446 } 02447 catch( cloudy_exit ) 02448 { 02449 goto error; 02450 } 02451 02452 if( n == 1 ) 02453 { 02454 fprintf( ioOut, " %ld\n", VERSION_ASCII ); 02455 fprintf( ioOut, " %d\n", ( ngrids == 1 ? 2 : 3 ) ); 02456 fprintf( ioOut, " %d\n", ( ngrids == 1 ? 2 : 3 ) ); 02457 fprintf( ioOut, " Teff\n" ); 02458 fprintf( ioOut, " log(g)\n" ); 02459 if( ngrids == 2 ) 02460 fprintf( ioOut, " log(Z)\n" ); 02461 else if( ngrids == 11 ) 02462 fprintf( ioOut, " f(He)\n" ); 02463 /* NB - this is based on the assumption that each of the planes in the cubic grid is the same */ 02464 fprintf( ioOut, " %ld\n", nmods*ngrids ); 02465 fprintf( ioOut, " %d\n", NRAUCH ); 02466 /* Rauch models give the wavelength in Angstrom */ 02467 fprintf( ioOut, " lambda\n" ); 02468 /* conversion factor to Angstrom */ 02469 fprintf( ioOut, " %.8e\n", 1. ); 02470 /* Rauch models give the "Astrophysical" flux F_lambda in erg/cm^2/s/cm */ 02471 fprintf( ioOut, " F_lambda\n" ); 02472 /* the factor PI*1e-8 is needed to convert to "regular" flux in erg/cm^2/s/Angstrom */ 02473 fprintf( ioOut, " %.8e\n", PI*1.e-8 ); 02474 /* NB - this is based on the assumption that each of the planes in the cubic grid is the same */ 02475 for( j=0; j < ngrids; j++ ) 02476 { 02477 /* write out the {Teff,log(g)} grid */ 02478 for( i=0; i < nmods; i++ ) 02479 { 02480 if( ngrids == 1 ) 02481 fprintf( ioOut, " %.0f %.1f", telg[i].par[0], telg[i].par[1] ); 02482 else 02483 fprintf( ioOut, " %.0f %.1f %.1f", telg[i].par[0], telg[i].par[1], par2[j] ); 02484 if( ((i+1)%4) == 0 ) 02485 fprintf( ioOut, "\n" ); 02486 } 02487 if( (i%4) != 0 ) 02488 fprintf( ioOut, "\n" ); 02489 } 02490 02491 fprintf( ioQQQ, " Writing: " ); 02492 } 02493 02494 for( i=0; i < nmods; i++ ) 02495 { 02496 /* must create name of next stellar atmosphere */ 02497 if( format == 1 ) 02498 sprintf( chLine, "%7.7ld_%2ld", (long)(telg[i].par[0]+0.5), (long)(10.*telg[i].par[1]+0.5) ); 02499 else if( format == 2 ) 02500 sprintf( chLine, "%7.7ld_%.2f", (long)(telg[i].par[0]+0.5), telg[i].par[1] ); 02501 else 02502 { 02503 fprintf( ioQQQ, " insanity in RauchInitializeSub\n" ); 02504 ShowMe(); 02505 cdEXIT(EXIT_FAILURE); 02506 } 02507 string chFileName( chLine ); 02508 chFileName += chSuff; 02509 /* now open next stellar atmosphere for reading*/ 02510 try 02511 { 02512 ioIn = open_data( chFileName.c_str(), "r", AS_LOCAL_ONLY ); 02513 } 02514 catch( cloudy_exit ) 02515 { 02516 goto error; 02517 } 02518 02519 /* get first line */ 02520 j = 0; 02521 if( read_whole_line( chLine, (int)sizeof(chLine), ioIn ) == NULL ) 02522 { 02523 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n", 02524 i, j ); 02525 goto error; 02526 } 02527 /* >>chng 02 nov 20, now keep reading them until don't hit the * 02528 * since number of comments may change */ 02529 while( chLine[0] == '*' ) 02530 { 02531 if( read_whole_line( chLine, (int)sizeof(chLine), ioIn ) == NULL ) 02532 { 02533 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n", 02534 i, j ); 02535 goto error; 02536 } 02537 ++j; 02538 } 02539 02540 for( j=0; j < NRAUCH; j++ ) 02541 { 02542 double ttemp, wl; 02543 /* get the input line */ 02544 /* >>chng 02 nov 20, don't reread very first line image since we got it above */ 02545 if( j > 0 ) 02546 { 02547 if(read_whole_line( chLine, (int)sizeof(chLine), ioIn )==NULL ) 02548 { 02549 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n", 02550 i, j ); 02551 goto error; 02552 } 02553 } 02554 02555 /* scan off wavelength and flux)*/ 02556 if( sscanf( chLine, "%lf %le", &wl, &ttemp ) != 2 ) 02557 { 02558 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n", 02559 i, j ); 02560 goto error; 02561 } 02562 02563 if( i == 0 ) 02564 wavl[j] = wl; 02565 else 02566 { 02567 /* check if this model is on the same wavelength grid as the first */ 02568 if( fabs(wavl[j]-wl) > 10.*DBL_EPSILON*wl ) 02569 { 02570 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n", 02571 i, j ); 02572 goto error; 02573 } 02574 } 02575 fluxes[j] = ttemp; 02576 } 02577 02578 /* finished - close the unit */ 02579 fclose(ioIn); 02580 02581 /* now write to output file */ 02582 if( i == 0 && n == 1 ) 02583 { 02584 /* wavelength grid is the same for all models, so write only once */ 02585 for( j=0; j < NRAUCH; j++ ) 02586 { 02587 fprintf( ioOut, " %.4e", wavl[j] ); 02588 /* want to have 5 numbers per line */ 02589 if( ((j+1)%5) == 0 ) 02590 fprintf( ioOut, "\n" ); 02591 } 02592 /* need to throw a newline if we did not end on an exact line */ 02593 if( (j%5) != 0 ) 02594 fprintf( ioOut, "\n" ); 02595 } 02596 02597 for( j=0; j < NRAUCH; j++ ) 02598 { 02599 fprintf( ioOut, " %.4e", fluxes[j] ); 02600 /* want to have 5 numbers per line */ 02601 if( ((j+1)%5) == 0 ) 02602 fprintf( ioOut, "\n" ); 02603 } 02604 /* need to throw a newline if we did not end on an exact line */ 02605 if( (j%5) != 0 ) 02606 fprintf( ioOut, "\n" ); 02607 02608 /* print to screen to show progress */ 02609 fprintf( ioQQQ, "." ); 02610 fflush( ioQQQ ); 02611 } 02612 02613 if( n == ngrids ) 02614 fprintf( ioQQQ, " done.\n" ); 02615 02616 fclose(ioOut); 02617 02618 /* free the space grabbed above */ 02619 FREE_CHECK( fluxes ); 02620 FREE_CHECK( wavl ); 02621 return 0; 02622 02623 error: 02624 FREE_CHECK( fluxes ); 02625 FREE_CHECK( wavl ); 02626 return 1; 02627 } 02628 02629 /* lgCompileAtmosphere does the actual rebinning onto the Cloudy grid and writes the binary file */ 02630 /* >>chng 01 feb 12, added return value to indicate success (0) or failure (1) */ 02631 STATIC bool lgCompileAtmosphere(const char chFNameIn[], 02632 const char chFNameOut[], 02633 const realnum Edges[], /* Edges[nedges] */ 02634 long nedges, 02635 process_counter& pc) 02636 { 02637 FILE *ioIN; /* used for input */ 02638 FILE *ioOUT; /* used for output */ 02639 02640 char chDataType[11]; 02641 char names[MDIM][MNAM+1]; 02642 02643 bool lgFreqX, lgFreqY, lgFlip; 02644 int32 val[7]; 02645 uint32 uval[2]; 02646 long int i, imod, version, nd, ndim, npar, nmods, ngrid; 02647 02648 /* these will be malloced into large work arrays */ 02649 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL, *scratch = NULL; 02650 02651 double convert_wavl, convert_flux; 02652 02653 mpp *telg = NULL; 02654 02655 DEBUG_ENTRY( "lgCompileAtmosphere()" ); 02656 02657 try 02658 { 02659 ioIN = open_data( chFNameIn, "r", AS_LOCAL_ONLY ); 02660 } 02661 catch( cloudy_exit ) 02662 { 02663 goto error; 02664 } 02665 fprintf( ioQQQ, " lgCompileAtmosphere got %s.\n", chFNameIn ); 02666 02667 /* read version number */ 02668 if( fscanf( ioIN, "%ld", &version ) != 1 ) 02669 { 02670 fprintf( ioQQQ, " lgCompileAtmosphere failed reading VERSION.\n" ); 02671 goto error; 02672 } 02673 02674 if( version != VERSION_ASCII ) 02675 { 02676 fprintf( ioQQQ, " lgCompileAtmosphere: there is a version number mismatch in" 02677 " the ascii atmosphere file: %s.\n", chFNameIn ); 02678 fprintf( ioQQQ, " lgCompileAtmosphere: Please recreate this file or download the" 02679 " latest version following the instructions on the Cloudy website.\n" ); 02680 goto error; 02681 } 02682 02683 /* >>chng 06 jun 10, read the dimension of the grid, PvH */ 02684 if( fscanf( ioIN, "%ld", &ndim ) != 1 || ndim <= 0 || ndim > MDIM ) 02685 { 02686 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid dimension of grid.\n" ); 02687 goto error; 02688 } 02689 02690 /* >>chng 06 jun 12, read the number of model parameters, PvH */ 02691 if( fscanf( ioIN, "%ld", &npar ) != 1 || npar <= 0 || npar < ndim || npar > MDIM ) 02692 { 02693 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid no. of model parameters.\n" ); 02694 goto error; 02695 } 02696 02697 /* make sure valgrind doesn't trip on the binary write of this array */ 02698 memset( names, '\0', MDIM*(MNAM+1) ); 02699 02700 for( nd=0; nd < npar; nd++ ) 02701 { 02702 if( fscanf( ioIN, "%6s", names[nd] ) != 1 ) 02703 { 02704 fprintf( ioQQQ, " lgCompileAtmosphere failed reading parameter label.\n" ); 02705 goto error; 02706 } 02707 } 02708 if( strcmp( names[0], "Teff" ) != 0 && strcmp( names[0], "Age" ) != 0 ) 02709 { 02710 fprintf( ioQQQ, " lgCompileAtmosphere: first parameter must be \"Teff\" or \"Age\".\n" ); 02711 goto error; 02712 } 02713 if( ndim > 1 && strcmp( names[0], "Age" ) == 0 ) 02714 { 02715 fprintf( ioQQQ, " First parameter is \"Age\", but ndim is not 1.\n" ); 02716 goto error; 02717 } 02718 if( ndim >= 2 && strcmp( names[1], "log(g)" ) != 0 ) 02719 { 02720 fprintf( ioQQQ, " lgCompileAtmosphere: second parameter must be \"log(g)\".\n" ); 02721 goto error; 02722 } 02723 02724 /* >>chng 05 nov 18, read the following extra parameters from the ascii file, PvH */ 02725 if( fscanf( ioIN, "%ld", &nmods ) != 1 || nmods <= 0 ) 02726 { 02727 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid number of models.\n" ); 02728 goto error; 02729 } 02730 02731 if( fscanf( ioIN, "%ld", &ngrid ) != 1 || ngrid <= 1 ) 02732 { 02733 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid number of grid points.\n" ); 02734 goto error; 02735 } 02736 02737 /* read data type for wavelengths, allowed values are lambda, nu */ 02738 if( fscanf( ioIN, "%10s", chDataType ) != 1 ) 02739 { 02740 fprintf( ioQQQ, " lgCompileAtmosphere failed reading wavl DataType string.\n" ); 02741 goto error; 02742 } 02743 02744 if( strcmp( chDataType, "lambda" ) == 0 ) 02745 lgFreqX = false; 02746 else if( strcmp( chDataType, "nu" ) == 0 ) 02747 lgFreqX = true; 02748 else { 02749 fprintf( ioQQQ, " lgCompileAtmosphere found illegal wavl DataType: %s.\n", chDataType ); 02750 goto error; 02751 } 02752 02753 if( fscanf( ioIN, "%le", &convert_wavl ) != 1 || convert_wavl <= 0. ) 02754 { 02755 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid wavl conversion factor.\n" ); 02756 goto error; 02757 } 02758 02759 /* read data type for flux, allowed values F_lambda, H_lambda, F_nu, H_nu */ 02760 if( fscanf( ioIN, "%10s", chDataType ) != 1 ) 02761 { 02762 fprintf( ioQQQ, " lgCompileAtmosphere failed reading flux DataType string.\n" ); 02763 goto error; 02764 } 02765 02766 if( strcmp( chDataType, "F_lambda" ) == 0 || strcmp( chDataType, "H_lambda" ) == 0 ) 02767 lgFreqY = false; 02768 else if( strcmp( chDataType, "F_nu" ) == 0 || strcmp( chDataType, "H_nu" ) == 0 ) 02769 lgFreqY = true; 02770 else { 02771 fprintf( ioQQQ, " lgCompileAtmosphere found illegal flux DataType: %s.\n", chDataType ); 02772 goto error; 02773 } 02774 02775 if( fscanf( ioIN, "%le", &convert_flux ) != 1 || convert_flux <= 0. ) 02776 { 02777 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid flux conversion factor.\n" ); 02778 goto error; 02779 } 02780 02781 telg = (mpp *)CALLOC( (size_t)nmods, sizeof(mpp) ); 02782 02783 for( i=0; i < nmods; i++ ) 02784 { 02785 for( nd=0; nd < npar; nd++ ) 02786 { 02787 if( fscanf( ioIN, "%le", &telg[i].par[nd] ) != 1 ) 02788 { 02789 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid model parameter.\n" ); 02790 goto error; 02791 } 02792 } 02793 if( telg[i].par[0] <= 0. ) 02794 { 02795 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid %s.\n", names[0] ); 02796 goto error; 02797 } 02798 } 02799 02800 try 02801 { 02802 ioOUT = open_data( chFNameOut, "wb", AS_LOCAL_ONLY ); 02803 } 02804 catch( cloudy_exit ) 02805 { 02806 goto error; 02807 } 02808 02809 /* >>chng 05 nov 16, use int32 instead of long so that binary file is compatible for 32/64-bit code */ 02810 /* >>chng 05 nov 18, add several parameters to binary file, change sequence, PvH */ 02811 /* >>chng 06 jun 10, add several parameters to binary file, change sequence, PvH */ 02812 val[0] = (int32)VERSION_BIN; 02813 val[1] = (int32)MDIM; 02814 val[2] = (int32)MNAM; 02815 val[3] = (int32)ndim; 02816 val[4] = (int32)npar; 02817 val[5] = (int32)nmods; 02818 val[6] = (int32)rfield.nupper; 02819 uval[0] = sizeof(val) + sizeof(uval) + sizeof(names) + nmods*sizeof(mpp); /* nOffset */ 02820 uval[1] = rfield.nupper*sizeof(realnum); /* nBlocksize */ 02821 02822 if( fwrite( val, sizeof(val), 1, ioOUT ) != 1 || 02823 fwrite( uval, sizeof(uval), 1, ioOUT ) != 1 || 02824 fwrite( names, sizeof(names), 1, ioOUT ) != 1 || 02825 /* write out the array of {Teff,log(g)} pairs */ 02826 fwrite( telg, sizeof(mpp), (size_t)nmods, ioOUT ) != (size_t)nmods || 02827 /* write out the cloudy energy grid for later sanity checks */ 02828 fwrite( rfield.AnuOrg, (size_t)uval[1], 1, ioOUT ) != 1 ) 02829 { 02830 fprintf( ioQQQ, " lgCompileAtmosphere failed writing header of output file.\n" ); 02831 goto error; 02832 } 02833 02834 /* MALLOC some workspace */ 02835 StarEner = (realnum *)MALLOC( sizeof(realnum)*ngrid ); 02836 scratch = (realnum *)MALLOC( sizeof(realnum)*ngrid ); 02837 StarFlux = (realnum *)MALLOC( sizeof(realnum)*ngrid ); 02838 CloudyFlux = (realnum *)MALLOC( (size_t)uval[1] ); 02839 02840 /* read wavelength grid */ 02841 for( i=0; i < ngrid; i++ ) 02842 { 02843 double help; 02844 if( fscanf( ioIN, "%lg", &help ) != 1 ) 02845 { 02846 fprintf( ioQQQ, " lgCompileAtmosphere failed reading wavelength.\n" ); 02847 goto error; 02848 } 02849 /* this conversion makes sure that scratch[i] is 02850 * either wavelength in Angstrom or frequency in Hz */ 02851 scratch[i] = (realnum)(help*convert_wavl); 02852 02853 ASSERT( scratch[i] > 0.f ); 02854 } 02855 02856 lgFlip = ( !lgFreqX && scratch[0] < scratch[1] ) || ( lgFreqX && scratch[0] > scratch[1] ); 02857 02858 /* convert continuum over to increasing frequency in Ryd */ 02859 for( i=0; i < ngrid; i++ ) 02860 { 02861 /* convert scratch[i] to frequency in Ryd */ 02862 if( lgFreqX ) 02863 scratch[i] /= (realnum)FR1RYD; 02864 else 02865 scratch[i] = (realnum)(RYDLAM/scratch[i]); 02866 02867 if( lgFlip ) 02868 StarEner[ngrid-i-1] = scratch[i]; 02869 else 02870 StarEner[i] = scratch[i]; 02871 } 02872 02873 ASSERT( StarEner[0] > 0.f ); 02874 /* make sure the array is in ascending order */ 02875 for( i=1; i < ngrid; i++ ) 02876 ASSERT( StarEner[i] > StarEner[i-1] ); 02877 02878 fprintf( ioQQQ, " Compiling: " ); 02879 02880 for( imod=0; imod < nmods; imod++ ) 02881 { 02882 const realnum CONVERT_FNU = (realnum)(1.e8*SPEEDLIGHT/POW2(FR1RYD)); 02883 02884 /* now read the stellar fluxes */ 02885 for( i=0; i < ngrid; i++ ) 02886 { 02887 double help; 02888 if( fscanf( ioIN, "%lg", &help ) != 1 ) 02889 { 02890 fprintf( ioQQQ, " lgCompileAtmosphere failed reading star flux.\n" ); 02891 goto error; 02892 } 02893 /* this conversion makes sure that scratch[i] is either 02894 * F_nu in erg/cm^2/s/Hz or F_lambda in erg/cm^2/s/A */ 02895 scratch[i] = (realnum)(help*convert_flux); 02896 02897 /* this can underflow on the Wien tail */ 02898 ASSERT( scratch[i] >= 0.f ); 02899 } 02900 02901 for( i=0; i < ngrid; i++ ) 02902 { 02903 if( lgFlip ) 02904 StarFlux[ngrid-i-1] = scratch[i]; 02905 else 02906 StarFlux[i] = scratch[i]; 02907 } 02908 02909 for( i=0; i < ngrid; i++ ) 02910 { 02911 /* this converts to F_nu in erg/cm^2/s/Hz */ 02912 if( !lgFreqY ) 02913 StarFlux[i] *= CONVERT_FNU/POW2(StarEner[i]); 02914 ASSERT( StarFlux[i] >= 0.f ); 02915 } 02916 02917 /* the re-binned values are returned in the "CloudyFlux" array */ 02918 RebinAtmosphere( ngrid, StarEner, StarFlux, CloudyFlux, nedges, Edges ); 02919 02920 /* write the continuum out as a binary file */ 02921 if( fwrite( CloudyFlux, (size_t)uval[1], 1, ioOUT ) != 1 ) 02922 { 02923 fprintf( ioQQQ, " lgCompileAtmosphere failed writing star flux.\n" ); 02924 goto error; 02925 } 02926 02927 fprintf( ioQQQ, "." ); 02928 fflush( ioQQQ ); 02929 } 02930 02931 fprintf( ioQQQ, " done.\n" ); 02932 02933 fclose(ioIN); 02934 fclose(ioOUT); 02935 02936 /* now free up the memory we claimed */ 02937 FREE_CHECK( CloudyFlux ); 02938 FREE_CHECK( StarFlux ); 02939 FREE_CHECK( StarEner ); 02940 FREE_CHECK( scratch ); 02941 FREE_CHECK( telg ); 02942 02943 fprintf( ioQQQ, " lgCompileAtmosphere completed ok.\n\n" ); 02944 ++pc.nOK; 02945 return false; 02946 02947 error: 02948 FREE_SAFE( CloudyFlux ); 02949 FREE_SAFE( StarFlux ); 02950 FREE_SAFE( StarEner ); 02951 FREE_SAFE( scratch ); 02952 FREE_SAFE( telg ); 02953 ++pc.nFail; 02954 return true; 02955 } 02956 02957 STATIC void InitGrid(stellar_grid *grid, 02958 bool lgList) 02959 { 02960 long nd; 02961 int32 version, mdim, mnam; 02962 02963 DEBUG_ENTRY( "InitGrid()" ); 02964 02965 try 02966 { 02967 grid->ioIN = open_data( grid->name.c_str(), "rb", grid->scheme ); 02968 } 02969 catch( cloudy_exit ) 02970 { 02971 /* something went wrong */ 02972 /* NB NB - DO NOT CHANGE THE FOLLOWING ERROR MESSAGE! checkall.pl picks it up */ 02973 fprintf( ioQQQ, " Error: stellar atmosphere file not found.\n" ); 02974 fprintf(ioQQQ , "\n\n If the path is set then it is possible that the stellar atmosphere data files do not exist.\n"); 02975 fprintf(ioQQQ , " Have the stellar data files been downloaded and compiled with the COMPILE STARS command?\n"); 02976 fprintf(ioQQQ , " If you are simply running the test suite and do not need the stellar continua then you should simply ignore this failure\n"); 02977 cdEXIT(EXIT_FAILURE); 02978 } 02979 02980 /* >>chng 01 oct 17, add version and size to this array */ 02981 if( fread( &version, sizeof(version), 1, grid->ioIN ) != 1 || 02982 fread( &mdim, sizeof(mdim), 1, grid->ioIN ) != 1 || 02983 fread( &mnam, sizeof(mnam), 1, grid->ioIN ) != 1 || 02984 fread( &grid->ndim, sizeof(grid->ndim), 1, grid->ioIN ) != 1 || 02985 fread( &grid->npar, sizeof(grid->npar), 1, grid->ioIN ) != 1 || 02986 fread( &grid->nmods, sizeof(grid->nmods), 1, grid->ioIN ) != 1 || 02987 fread( &grid->ngrid, sizeof(grid->ngrid), 1, grid->ioIN ) != 1 || 02988 fread( &grid->nOffset, sizeof(grid->nOffset), 1, grid->ioIN ) != 1 || 02989 fread( &grid->nBlocksize, sizeof(grid->nBlocksize), 1, grid->ioIN ) != 1 ) 02990 { 02991 fprintf( ioQQQ, " InitGrid failed reading header.\n" ); 02992 cdEXIT(EXIT_FAILURE); 02993 } 02994 02995 /* do some sanity checks */ 02996 if( version != VERSION_BIN ) 02997 { 02998 fprintf( ioQQQ, " InitGrid: there is a version mismatch between" 02999 " the compiled atmospheres file I expected and the one I found.\n" ); 03000 fprintf( ioQQQ, " InitGrid: Please recompile the stellar" 03001 " atmospheres file with the command: %s.\n", grid->command ); 03002 cdEXIT(EXIT_FAILURE); 03003 } 03004 03005 if( mdim != MDIM || mnam != MNAM ) 03006 { 03007 fprintf( ioQQQ, " InitGrid: the compiled atmospheres file is produced" 03008 " with an incompatible version of Cloudy.\n" ); 03009 fprintf( ioQQQ, " InitGrid: Please recompile the stellar" 03010 " atmospheres file with the command: %s.\n", grid->command ); 03011 cdEXIT(EXIT_FAILURE); 03012 } 03013 03014 ASSERT( grid->ndim > 0 && grid->ndim <= MDIM ); 03015 ASSERT( grid->npar >= grid->ndim && grid->npar <= MDIM ); 03016 ASSERT( grid->nmods > 0 ); 03017 ASSERT( grid->ngrid > 0 ); 03018 ASSERT( grid->nOffset > 0 ); 03019 ASSERT( grid->nBlocksize > 0 ); 03020 03021 rfield.nupper = grid->ngrid; 03022 03023 if( fread( &grid->names, sizeof(grid->names), 1, grid->ioIN ) != 1 ) 03024 { 03025 fprintf( ioQQQ, " InitGrid failed reading names array.\n" ); 03026 cdEXIT(EXIT_FAILURE); 03027 } 03028 03029 grid->telg = (mpp *)MALLOC( sizeof(mpp)*grid->nmods ); 03030 grid->val = (double **)MALLOC( sizeof(double*)*grid->ndim ); 03031 for( nd=0; nd < grid->ndim; nd++ ) 03032 { 03033 grid->val[nd] = (double *)MALLOC( sizeof(double)*grid->nmods ); 03034 } 03035 grid->nval = (long *)MALLOC( sizeof(long)*grid->ndim ); 03036 03037 if( fread( grid->telg, sizeof(mpp), grid->nmods, grid->ioIN ) != (size_t)grid->nmods ) 03038 { 03039 fprintf( ioQQQ, " InitGrid failed reading model parameter block.\n" ); 03040 cdEXIT(EXIT_FAILURE); 03041 } 03042 03043 # ifdef SEEK_END 03044 /* sanity check: does the file have the correct length ? */ 03045 /* NOTE: this operation is not necessarily supported by all operating systems 03046 * but if the preprocessor symbol SEEK_END exists it is assumed to be supported */ 03047 int res = fseek( grid->ioIN, 0, SEEK_END ); 03048 if( res == 0 ) 03049 { 03050 long End = ftell( grid->ioIN ); 03051 long Expected = grid->nOffset + (grid->nmods+1)*grid->nBlocksize; 03052 if( End != Expected ) 03053 { 03054 fprintf( ioQQQ, " InitGrid: Problem performing sanity check for size of binary file.\n" ); 03055 fprintf( ioQQQ, " InitGrid: I expected to find %ld bytes, but actually found %ld bytes.\n", 03056 Expected, End ); 03057 fprintf( ioQQQ, " InitGrid: Please recompile the stellar" 03058 " atmospheres file with the command: %s.\n", grid->command ); 03059 cdEXIT(EXIT_FAILURE); 03060 } 03061 } 03062 # endif 03063 03064 InitIndexArrays( grid, lgList ); 03065 03066 /* set default interpolation mode */ 03067 grid->imode = IM_RECT_GRID; 03068 /* these are only used by CoStar grids */ 03069 grid->trackLen = NULL; 03070 grid->nTracks = 0; 03071 grid->jval = NULL; 03072 return; 03073 } 03074 03075 /* check whether a binary atmosphere exists and is up-to-date */ 03076 STATIC bool lgValidBinFile(const char *binName, process_counter& pc, access_scheme scheme) 03077 { 03078 int32 version, mdim, mnam; 03079 stellar_grid grid; 03080 03081 DEBUG_ENTRY( "lgValidBinFile()" ); 03082 03083 grid.name = binName; 03084 03085 if( (grid.ioIN = open_data( grid.name.c_str(), "rb", scheme )) == NULL ) 03086 return false; 03087 03088 if( fread( &version, sizeof(version), 1, grid.ioIN ) != 1 || 03089 fread( &mdim, sizeof(mdim), 1, grid.ioIN ) != 1 || 03090 fread( &mnam, sizeof(mnam), 1, grid.ioIN ) != 1 || 03091 fread( &grid.ndim, sizeof(grid.ndim), 1, grid.ioIN ) != 1 || 03092 fread( &grid.npar, sizeof(grid.npar), 1, grid.ioIN ) != 1 || 03093 fread( &grid.nmods, sizeof(grid.nmods), 1, grid.ioIN ) != 1 || 03094 fread( &grid.ngrid, sizeof(grid.ngrid), 1, grid.ioIN ) != 1 || 03095 fread( &grid.nOffset, sizeof(grid.nOffset), 1, grid.ioIN ) != 1 || 03096 fread( &grid.nBlocksize, sizeof(grid.nBlocksize), 1, grid.ioIN ) != 1 ) 03097 { 03098 fclose( grid.ioIN ); 03099 return false; 03100 } 03101 03102 /* do some sanity checks */ 03103 if( version != VERSION_BIN || mdim != MDIM || mnam != MNAM ) 03104 { 03105 fclose( grid.ioIN ); 03106 return false; 03107 } 03108 03109 # ifdef SEEK_END 03110 /* sanity check: does the file have the correct length ? */ 03111 /* NOTE: this operation is not necessarily supported by all operating systems 03112 * but if the preprocessor symbol SEEK_END exists it is assumed to be supported */ 03113 int res = fseek( grid.ioIN, 0, SEEK_END ); 03114 if( res == 0 ) 03115 { 03116 long End = ftell( grid.ioIN ); 03117 long Expected = grid.nOffset + (grid.nmods+1)*grid.nBlocksize; 03118 if( End != Expected ) 03119 { 03120 fclose( grid.ioIN ); 03121 return false; 03122 } 03123 } 03124 # endif 03125 03126 fclose( grid.ioIN ); 03127 ++pc.notProcessed; // the file is up-to-date -> no processing 03128 return true; 03129 } 03130 03131 /* check whether a ascii atmosphere file exists and is up-to-date */ 03132 STATIC bool lgValidAsciiFile(const char *ascName, access_scheme scheme) 03133 { 03134 long version; 03135 FILE *ioIN; 03136 03137 DEBUG_ENTRY( "lgValidAsciiFile()" ); 03138 03139 /* can we read the file? */ 03140 if( (ioIN = open_data( ascName, "r", scheme )) == NULL ) 03141 return false; 03142 03143 /* check version number */ 03144 if( fscanf( ioIN, "%ld", &version ) != 1 || version != VERSION_ASCII ) 03145 { 03146 fclose( ioIN ); 03147 return false; 03148 } 03149 03150 fclose( ioIN ); 03151 return true; 03152 } 03153 03154 /* sort CoStar models according to track and index number, store indices in grid->jval[] */ 03155 STATIC void InitGridCoStar(stellar_grid *grid) /* the grid parameters */ 03156 { 03157 char track; 03158 bool lgFound; 03159 long i, index[2]; 03160 03161 DEBUG_ENTRY( "InitGridCoStar()" ); 03162 03163 ASSERT( grid->ndim == 2 ); 03164 ASSERT( grid->jlo != NULL && grid->jhi != NULL ); 03165 03166 grid->jval = grid->jlo; 03167 FREE_CHECK( grid->jhi ); 03168 grid->jlo = grid->jhi = NULL; 03169 03170 /* invalidate contents set by InitGrid first */ 03171 memset( grid->jval, 0xff, (size_t)(grid->nval[0]*grid->nval[1]*sizeof(long)) ); 03172 03173 grid->trackLen = (long *)CALLOC( (size_t)grid->nmods, sizeof(long) ); 03174 03175 index[1] = 0; 03176 while( true ) 03177 { 03178 index[0] = 0; 03179 track = (char)('A'+index[1]); 03180 do 03181 { 03182 lgFound = false; 03183 for( i=0; i < grid->nmods; i++ ) 03184 { 03185 if( grid->telg[i].chGrid == track && grid->telg[i].modid == index[0]+1 ) 03186 { 03187 grid->jval[JIndex(grid,index)] = i; 03188 ++index[0]; 03189 lgFound = true; 03190 break; 03191 } 03192 } 03193 } 03194 while( lgFound ); 03195 03196 if( index[0] == 0 ) 03197 break; 03198 03199 grid->trackLen[index[1]] = index[0]; 03200 ++index[1]; 03201 } 03202 03203 grid->nTracks = index[1]; 03204 return; 03205 } 03206 03207 STATIC void CheckVal(const stellar_grid *grid, 03208 double val[], /* val[ndim] */ 03209 long *nval, 03210 long *ndim) 03211 { 03212 DEBUG_ENTRY( "CheckVal()" ); 03213 03214 if( *ndim == 0 ) 03215 *ndim = (long)grid->ndim; 03216 if( *ndim == 2 && *nval == 1 ) 03217 { 03218 /* default gravity is maximum gravity */ 03219 val[*nval] = grid->val[1][grid->nval[1]-1]; 03220 ++(*nval); 03221 } 03222 if( *ndim != (long)grid->ndim ) 03223 { 03224 fprintf( ioQQQ, " A %ld-dim grid was requested, but a %ld-dim grid was found.\n", 03225 *ndim, (long)grid->ndim ); 03226 cdEXIT(EXIT_FAILURE); 03227 } 03228 if( *nval < *ndim ) 03229 { 03230 fprintf( ioQQQ, " A %ld-dim grid was requested, but only %ld parameters were entered.\n", 03231 *ndim, *nval ); 03232 cdEXIT(EXIT_FAILURE); 03233 } 03234 } 03235 03236 STATIC void InterpolateRectGrid(const stellar_grid *grid, 03237 const double val[], /* val[ndim] */ 03238 double *Tlow, 03239 double *Thigh) 03240 { 03241 bool lgInvalid; 03242 long int i, 03243 *indlo, 03244 *indhi, 03245 *index, 03246 k, 03247 nd; 03248 double *aval; 03249 03250 DEBUG_ENTRY( "InterpolateRectGrid()" ); 03251 03252 /* create some space */ 03253 indlo = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) ); 03254 indhi = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) ); 03255 index = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) ); 03256 aval = (double *)MALLOC((size_t)(grid->npar*sizeof(double)) ); 03257 03258 ASSERT( rfield.lgContMalloc[rfield.nspec] ); 03259 ASSERT( grid->nBlocksize == rfield.nupper*sizeof(realnum) ); 03260 03261 /* save energy scale for check against code's in conorm (scale not yet defined when this routine called) */ 03262 GetModel( grid, -1, rfield.tNuRyd[rfield.nspec], lgSILENT, lgLINEAR ); 03263 03264 # if DEBUGPRT 03265 /* check whether the models have the correct effective temperature, for debugging only */ 03266 ValidateGrid( grid, 0.02 ); 03267 # endif 03268 03269 /* now generate pointers for models to use */ 03270 for( nd=0; nd < grid->ndim; nd++ ) 03271 { 03272 FindIndex( grid->val[nd], grid->nval[nd], val[nd], &indlo[nd], &indhi[nd], &lgInvalid ); 03273 if( lgInvalid ) 03274 { 03275 fprintf( ioQQQ, 03276 " Requested parameter %s = %.2f is not within the range %.2f to %.2f\n", 03277 grid->names[nd], val[nd], grid->val[nd][0], grid->val[nd][grid->nval[nd]-1] ); 03278 cdEXIT(EXIT_FAILURE); 03279 } 03280 } 03281 03282 InterpolateModel( grid, val, aval, indlo, indhi, index, grid->ndim, rfield.tslop[rfield.nspec], IS_UNDEFINED ); 03283 03284 /* print the parameters of the interpolated model */ 03285 if( called.lgTalk ) 03286 { 03287 if( grid->npar == 1 ) 03288 fprintf( ioQQQ, 03289 " * c<< FINAL: %6s = %13.2f" 03290 " >>> *\n", 03291 grid->names[0], aval[0] ); 03292 else if( grid->npar == 2 ) 03293 fprintf( ioQQQ, 03294 " * c<< FINAL: %6s = %10.2f" 03295 " %6s = %8.5f >>> *\n", 03296 grid->names[0], aval[0], grid->names[1], aval[1] ); 03297 else if( grid->npar == 3 ) 03298 fprintf( ioQQQ, 03299 " * c<< FINAL: %6s = %7.0f" 03300 " %6s = %5.2f %6s = %5.2f >>> *\n", 03301 grid->names[0], aval[0], grid->names[1], aval[1], 03302 grid->names[2], aval[2] ); 03303 else if( grid->npar >= 4 ) 03304 { 03305 fprintf( ioQQQ, 03306 " * c<< FINAL: %4s = %7.0f" 03307 " %6s = %4.2f %6s = %5.2f %6s = ", 03308 grid->names[0], aval[0], grid->names[1], aval[1], 03309 grid->names[2], aval[2], grid->names[3] ); 03310 fprintf( ioQQQ, PrintEfmt( "%9.2e", aval[3] ) ); 03311 fprintf( ioQQQ, " >>> *\n" ); 03312 } 03313 } 03314 03315 for( i=0; i < rfield.nupper; i++ ) 03316 { 03317 rfield.tslop[rfield.nspec][i] = (realnum)pow((realnum)10.f,rfield.tslop[rfield.nspec][i]); 03318 if( rfield.tslop[rfield.nspec][i] < 1e-37 ) 03319 rfield.tslop[rfield.nspec][i] = 0.; 03320 } 03321 03322 if( false ) 03323 { 03324 FILE *ioBUG = fopen( "interpolated.txt", "w" ); 03325 for( k=0; k < rfield.nupper; ++k ) 03326 fprintf( ioBUG, "%e %e\n", rfield.tNuRyd[rfield.nspec][k], rfield.tslop[rfield.nspec][k] ); 03327 fclose( ioBUG ); 03328 } 03329 03330 if( strcmp( grid->names[0], "Teff" ) == 0 ) 03331 { 03332 if( ! lgValidModel( rfield.tNuRyd[rfield.nspec], rfield.tslop[rfield.nspec], val[0], 0.10 ) ) 03333 TotalInsanity(); 03334 } 03335 03336 /* set limits for optimizer */ 03337 SetLimits( grid, val[0], indlo, indhi, NULL, NULL, Tlow, Thigh ); 03338 03339 FREE_CHECK( aval ); 03340 FREE_CHECK( index ); 03341 FREE_CHECK( indhi ); 03342 FREE_CHECK( indlo ); 03343 return; 03344 } 03345 03346 STATIC void FreeGrid(stellar_grid *grid) 03347 { 03348 long i; 03349 03350 DEBUG_ENTRY( "FreeGrid()" ); 03351 03352 /* this was opened/allocated in InitGrid and subsidiaries, 03353 * this should become a destructor in C++ */ 03354 fclose( grid->ioIN ); 03355 FREE_CHECK( grid->telg ); 03356 for( i = 0; i < grid->ndim; i++ ) 03357 FREE_CHECK( grid->val[i] ); 03358 FREE_CHECK( grid->val ); 03359 FREE_CHECK( grid->nval ); 03360 FREE_SAFE( grid->jlo ); 03361 FREE_SAFE( grid->jhi ); 03362 FREE_SAFE( grid->trackLen ) 03363 FREE_SAFE( grid->jval ); 03364 return; 03365 } 03366 03367 STATIC void InterpolateModel(const stellar_grid *grid, 03368 const double val[], 03369 double aval[], 03370 const long indlo[], 03371 const long indhi[], 03372 long index[], 03373 long nd, 03374 realnum flux1[], 03375 IntStage stage) 03376 { 03377 bool lgDryRun; 03378 long i, ind, j; 03379 03380 DEBUG_ENTRY( "InterpolateModel()" ); 03381 03382 --nd; 03383 03384 lgDryRun = ( flux1 == NULL ); 03385 03386 if( nd < 0 ) 03387 { 03388 long n = JIndex(grid,index); 03389 if( stage == IS_FIRST ) 03390 ind = ( grid->jlo[n] >= 0 ) ? grid->jlo[n] : grid->jhi[n]; 03391 else if( stage == IS_SECOND ) 03392 ind = ( grid->jhi[n] >= 0 ) ? grid->jhi[n] : grid->jlo[n]; 03393 else if( grid->ndim == 1 ) 03394 /* in this case grid->jlo[n] and grid->jhi[n] should be identical */ 03395 ind = grid->jlo[n]; 03396 else 03397 TotalInsanity(); 03398 03399 if( ind < 0 ) 03400 { 03401 fprintf( ioQQQ, " The requested interpolation could not be completed, sorry.\n" ); 03402 fprintf( ioQQQ, " No suitable match was found for a model with" ); 03403 for( i=0; i < grid->ndim; i++ ) 03404 fprintf( ioQQQ, " %s=%.6g ", grid->names[i], grid->val[i][index[i]] ); 03405 fprintf( ioQQQ, "\n" ); 03406 cdEXIT(EXIT_FAILURE); 03407 } 03408 03409 for( i=0; i < grid->npar; i++ ) 03410 aval[i] = grid->telg[ind].par[i]; 03411 03412 if( !lgDryRun ) 03413 { 03414 for( i=0; i < grid->ndim && called.lgTalk; i++ ) 03415 { 03416 if( fabs(grid->val[i][index[i]]-aval[i]) > 10.*DBL_EPSILON*fabs(aval[i]) ) 03417 { 03418 fprintf( ioQQQ, " No exact match was found for a model with" ); 03419 for( j=0; j < grid->ndim; j++ ) 03420 fprintf( ioQQQ, " %s=%.6g ", grid->names[j], grid->val[j][index[j]] ); 03421 fprintf( ioQQQ, "- using the following model instead:\n" ); 03422 break; 03423 } 03424 } 03425 03426 GetModel( grid, ind, flux1, lgVERBOSE, lgTAKELOG ); 03427 } 03428 } 03429 else 03430 { 03431 realnum *flux2; 03432 double *aval2; 03433 03434 # ifndef NDEBUG 03435 const realnum SECURE = 10.f*FLT_EPSILON; 03436 # endif 03437 03438 flux2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) ); 03439 aval2 = (double*)MALLOC((size_t)(grid->npar*sizeof(double)) ); 03440 03441 /* Interpolation is carried out first in Teff (i.e., nd == 0), then the 03442 * parameter with nd == 1 (log(g)), etc. One or two atmosphere models 03443 * are read depending on whether the parameter was matched exactly or not. 03444 * If needed, logarithmic interpolation is done. 03445 */ 03446 03447 if( nd == 1 ) 03448 stage = IS_FIRST; 03449 03450 index[nd] = indlo[nd]; 03451 InterpolateModel( grid, val, aval, indlo, indhi, index, nd, flux1, stage ); 03452 03453 if( nd == 1 ) 03454 stage = IS_SECOND; 03455 03456 index[nd] = indhi[nd]; 03457 InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, NULL, stage ); 03458 03459 /* the test MUST be "greater than", otherwise it will fail for aval[nd] == aval2[nd] == 0 */ 03460 if( fabs(aval2[nd]-aval[nd]) > 10.*DBL_EPSILON*fabs(aval[nd]) ) 03461 { 03462 double fr1, fr2, fc1 = 0., fc2 = 0.; 03463 03464 if( !lgDryRun ) 03465 InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, flux2, stage ); 03466 03467 fr1 = (aval2[nd]-val[nd])/(aval2[nd]-aval[nd]); 03468 /* when interpolating in log(g) it can happen that fr1 is outside the range 0 .. 1 03469 * this can be the case when the requested log(g) was not present in the grid 03470 * and it had to be approximated by another model. In this case do not extrapolate */ 03471 if( nd == 1 ) 03472 fr1 = MIN2( MAX2( fr1, 0. ), 1. ); 03473 fr2 = 1. - fr1; 03474 03475 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE ); 03476 03477 if( !lgDryRun ) 03478 { 03479 # if DEBUGPRT 03480 fprintf( ioQQQ, "interpolation nd=%ld fr1=%g\n", nd, fr1 ); 03481 # endif 03482 03483 /* special treatment for high-temperature Rauch models */ 03484 if( nd == 0 && strcmp( grid->names[nd], "Teff" ) == 0 ) 03485 { 03486 /* The following is an approximate scaling to use for the range of 03487 * temperatures above 200000 K in the H-Ca Rauch models where the 03488 * temperature steps are large and thus the interpolations are over 03489 * large ranges. For the lower temperatures I assume that there is 03490 * no need for this. 03491 * 03492 * It should be remembered that this interpolation is not exact, and 03493 * the possible error at high temperatures might be large enough to 03494 * matter. (Kevin Volk) 03495 */ 03496 fc1 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->val[nd][indlo[nd]])*4. : 0.; 03497 fc2 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->val[nd][indhi[nd]])*4. : 0.; 03498 } 03499 03500 for( i=0; i < rfield.nupper; ++i ) 03501 flux1[i] = (realnum)(fr1*(flux1[i]+fc1) + fr2*(flux2[i]+fc2)); 03502 } 03503 03504 for( i=0; i < grid->npar; i++ ) 03505 aval[i] = fr1*aval[i] + fr2*aval2[i]; 03506 } 03507 03508 FREE_CHECK( aval2 ); 03509 FREE_CHECK( flux2 ); 03510 } 03511 return; 03512 } 03513 03514 STATIC void InterpolateModelCoStar(const stellar_grid *grid, 03515 const double val[], 03516 double aval[], 03517 const long indlo[], 03518 const long indhi[], 03519 long index[], 03520 long nd, 03521 long off, 03522 realnum flux1[]) 03523 { 03524 long i, ind; 03525 03526 DEBUG_ENTRY( "InterpolateModelCoStar()" ); 03527 03528 if( nd == 2 ) 03529 { 03530 ind = ( index[1] == 0 ) ? indlo[index[0]] : indhi[index[0]]; 03531 03532 GetModel( grid, ind, flux1, lgVERBOSE, lgTAKELOG ); 03533 03534 for( i=0; i < grid->npar; i++ ) 03535 aval[i] = grid->telg[ind].par[i]; 03536 } 03537 else 03538 { 03539 bool lgSkip; 03540 # ifndef NDEBUG 03541 const realnum SECURE = 10.f*FLT_EPSILON; 03542 # endif 03543 03544 /* Interpolation is carried out first along evolutionary tracks, then 03545 * in between evolutionary tracks. Between 1 and 4 atmosphere models are read 03546 * depending on whether the parameter/track was matched exactly or not. 03547 */ 03548 03549 index[nd] = 0; 03550 InterpolateModelCoStar( grid, val, aval, indlo, indhi, index, nd+1, off, flux1 ); 03551 03552 lgSkip = ( nd == 1 ) ? ( indhi[index[0]] == indlo[index[0]] ) : 03553 ( indlo[0] == indlo[1] && indhi[0] == indhi[1] ); 03554 03555 if( ! lgSkip ) 03556 { 03557 realnum *flux2; 03558 double fr1, fr2, *aval2; 03559 03560 flux2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) ); 03561 aval2 = (double*)MALLOC((size_t)(grid->npar*sizeof(double)) ); 03562 03563 index[nd] = 1; 03564 InterpolateModelCoStar( grid, val, aval2, indlo, indhi, index, nd+1, off, flux2 ); 03565 03566 fr1 = (aval2[nd+off]-val[nd])/(aval2[nd+off]-aval[nd+off]); 03567 fr2 = 1. - fr1; 03568 03569 # if DEBUGPRT 03570 fprintf( ioQQQ, "interpolation nd=%ld fr1=%g\n", nd, fr1 ); 03571 # endif 03572 03573 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE ); 03574 03575 for( i=0; i < rfield.nupper; ++i ) 03576 flux1[i] = (realnum)(fr1*flux1[i] + fr2*flux2[i]); 03577 03578 for( i=0; i < grid->npar; i++ ) 03579 aval[i] = fr1*aval[i] + fr2*aval2[i]; 03580 03581 FREE_CHECK( aval2 ); 03582 FREE_CHECK( flux2 ); 03583 } 03584 } 03585 return; 03586 } 03587 03588 STATIC void GetModel(const stellar_grid *grid, 03589 long ind, 03590 realnum flux[], 03591 bool lgTalk, 03592 bool lgTakeLog) 03593 { 03594 long i; 03595 03596 DEBUG_ENTRY( "GetModel()" ); 03597 03598 /* add 1 to account for frequency grid that is stored in front of all the atmospheres */ 03599 ind++; 03600 03601 /* make sure ident is exactly 12 characters long, otherwise output won't fit */ 03602 ASSERT( strlen(grid->ident) == 12 ); 03603 /* ind == 0 is the frequency grid, ind == 1 .. nmods are the atmosphere models */ 03604 ASSERT( ind >= 0 && ind <= grid->nmods ); 03605 03606 /* skip over ind stars */ 03607 /* >>chng 01 oct 18, add nOffset */ 03608 if( fseek( grid->ioIN, (long)(ind*grid->nBlocksize+grid->nOffset), SEEK_SET ) != 0 ) 03609 { 03610 fprintf( ioQQQ, " Error seeking atmosphere %ld\n", ind ); 03611 cdEXIT(EXIT_FAILURE); 03612 } 03613 03614 if( fread( flux, 1, grid->nBlocksize, grid->ioIN ) != grid->nBlocksize ) 03615 { 03616 fprintf( ioQQQ, " Error trying to read atmosphere %ld\n", ind ); 03617 cdEXIT(EXIT_FAILURE); 03618 } 03619 03620 /* print the parameters of the atmosphere model */ 03621 if( called.lgTalk && lgTalk ) 03622 { 03623 /* ind-1 below since telg doesn't have the entry for the frequency grid */ 03624 if( grid->npar == 1 ) 03625 fprintf( ioQQQ, 03626 " * c<< %s model%5ld read. " 03627 " %6s = %13.2f >>> *\n", 03628 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0] ); 03629 else if( grid->npar == 2 ) 03630 fprintf( ioQQQ, 03631 " * c<< %s model%5ld read. " 03632 " %6s = %10.2f %6s = %8.5f >>> *\n", 03633 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0], 03634 grid->names[1], grid->telg[ind-1].par[1] ); 03635 else if( grid->npar == 3 ) 03636 fprintf( ioQQQ, 03637 " * c<< %s model%5ld read. " 03638 " %6s=%7.0f %6s=%5.2f %6s=%5.2f >>> *\n", 03639 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0], 03640 grid->names[1], grid->telg[ind-1].par[1], 03641 grid->names[2], grid->telg[ind-1].par[2] ); 03642 else if( grid->npar >= 4 ) 03643 { 03644 fprintf( ioQQQ, 03645 " * c< %s mdl%4ld" 03646 " %4s=%5.0f %6s=%4.2f %6s=%5.2f %6s=", 03647 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0], 03648 grid->names[1], grid->telg[ind-1].par[1], 03649 grid->names[2], grid->telg[ind-1].par[2], grid->names[3] ); 03650 fprintf( ioQQQ, PrintEfmt( "%9.2e", grid->telg[ind-1].par[3] ) ); 03651 fprintf( ioQQQ, " >> *\n" ); 03652 } 03653 } 03654 03655 if( lgTakeLog ) 03656 { 03657 /* convert to logs since we will interpolate in log flux */ 03658 for( i=0; i < rfield.nupper; ++i ) 03659 flux[i] = (realnum)log10( MAX2( 1e-37, (double)flux[i] ) ); 03660 } 03661 return; 03662 } 03663 03664 STATIC void SetLimits(const stellar_grid *grid, 03665 double val, 03666 const long indlo[], 03667 const long indhi[], 03668 const long useTr[], 03669 const realnum ValTr[], 03670 double *loLim, 03671 double *hiLim) 03672 { 03673 DEBUG_ENTRY( "SetLimits()" ); 03674 03675 if( optimize.lgVarOn ) 03676 { 03677 int ptr0,ptr1; 03678 long index[MDIM], j; 03679 const double SECURE = (1. + 20.*(double)FLT_EPSILON); 03680 03681 *loLim = +DBL_MAX; 03682 *hiLim = -DBL_MAX; 03683 03684 switch( grid->imode ) 03685 { 03686 case IM_RECT_GRID: 03687 *loLim = -DBL_MAX; 03688 *hiLim = +DBL_MAX; 03689 SetLimitsSub( grid, val, indlo, indhi, index, grid->ndim, loLim, hiLim ); 03690 break; 03691 case IM_COSTAR_TEFF_MODID: 03692 case IM_COSTAR_TEFF_LOGG: 03693 case IM_COSTAR_MZAMS_AGE: 03694 for( j=0; j < grid->nTracks; j++ ) 03695 { 03696 if( ValTr[j] != -FLT_MAX ) 03697 { 03698 /* M_ZAMS is already logarithm, Teff is linear */ 03699 double temp = ( grid->imode == IM_COSTAR_MZAMS_AGE ) ? 03700 pow(10.,(double)ValTr[j]) : ValTr[j]; 03701 *loLim = MIN2(*loLim,temp); 03702 *hiLim = MAX2(*hiLim,temp); 03703 } 03704 } 03705 break; 03706 case IM_COSTAR_AGE_MZAMS: 03707 index[0] = 0; 03708 index[1] = useTr[0]; 03709 ptr0 = grid->jval[JIndex(grid,index)]; 03710 index[1] = useTr[1]; 03711 ptr1 = grid->jval[JIndex(grid,index)]; 03712 *loLim = MAX2(grid->telg[ptr0].par[3],grid->telg[ptr1].par[3]); 03713 # if DEBUGPRT 03714 printf( "set limit 0: (models %d, %d) %f %f\n", 03715 ptr0+1, ptr1+1, grid->telg[ptr0].par[3], grid->telg[ptr1].par[3] ); 03716 # endif 03717 index[0] = grid->trackLen[useTr[0]]-1; 03718 index[1] = useTr[0]; 03719 ptr0 = grid->jval[JIndex(grid,index)]; 03720 index[0] = grid->trackLen[useTr[1]]-1; 03721 index[1] = useTr[1]; 03722 ptr1 = grid->jval[JIndex(grid,index)]; 03723 *hiLim = MIN2(grid->telg[ptr0].par[3],grid->telg[ptr1].par[3]); 03724 # if DEBUGPRT 03725 printf( "set limit 1: (models %d, %d) %f %f\n", 03726 ptr0+1, ptr1+1, grid->telg[ptr0].par[3], grid->telg[ptr1].par[3] ); 03727 # endif 03728 break; 03729 default: 03730 fprintf( ioQQQ, " SetLimits called with insane value for imode: %d.\n", grid->imode ); 03731 cdEXIT(EXIT_FAILURE); 03732 } 03733 03734 ASSERT( fabs(*loLim) < DBL_MAX && fabs(*hiLim) < DBL_MAX ); 03735 03736 /* check sanity of optimization limits */ 03737 if( *hiLim <= *loLim ) 03738 { 03739 fprintf( ioQQQ, " no room to optimize: lower limit %.4f, upper limit %.4f.\n", 03740 *loLim,*hiLim ); 03741 cdEXIT(EXIT_FAILURE); 03742 } 03743 03744 /* make a bit of room for round-off errors */ 03745 *loLim *= SECURE; 03746 *hiLim /= SECURE; 03747 03748 # if DEBUGPRT 03749 printf("set limits: %g %g\n",*loLim,*hiLim); 03750 # endif 03751 } 03752 else 03753 { 03754 *loLim = 0.; 03755 *hiLim = 0.; 03756 } 03757 return; 03758 } 03759 03760 STATIC void SetLimitsSub(const stellar_grid *grid, 03761 double val, 03762 const long indlo[], 03763 const long indhi[], 03764 long index[], 03765 long nd, 03766 double *loLim, 03767 double *hiLim) 03768 { 03769 long n; 03770 03771 DEBUG_ENTRY( "SetLimitsSub()" ); 03772 03773 --nd; 03774 03775 if( nd < 2 ) 03776 { 03777 double loLoc = +DBL_MAX; 03778 double hiLoc = -DBL_MAX; 03779 03780 index[1] = 0; 03781 for( index[0]=0; index[0] < grid->nval[0]; ++index[0] ) 03782 { 03783 /* grid->val[0][i] is the array of Teff/Age values in the grid, which 03784 * it is sorted in strict monotonically increasing order. This routine 03785 * searches for the largest range [loLoc,hiLoc] in Teff/Age such that 03786 * loLoc <= val <= hiLoc (but loLoc < hiLoc) and at least one model 03787 * exists for each Teff/Age value in this range. This assures that 03788 * interpolation is safe and the optimizer will not trip... */ 03789 n = JIndex(grid,index); 03790 if( grid->jhi[n] < 0 ) 03791 { 03792 /* there are no models with this value of Teff */ 03793 /* this value of Teff should be outside of allowed range */ 03794 if( grid->val[0][index[0]] < val ) 03795 loLoc = DBL_MAX; 03796 if( grid->val[0][index[0]] > val ) 03797 break; 03798 } 03799 else 03800 { 03801 /* there are models with this value of Teff */ 03802 /* update range to include this value of Teff */ 03803 if( grid->val[0][index[0]] <= val ) 03804 { 03805 /* remember lowest legal value of loLoc 03806 * -> only update if previous value was illegal */ 03807 if( loLoc == DBL_MAX ) 03808 loLoc = grid->val[0][index[0]]; 03809 } 03810 if( grid->val[0][index[0]] >= val ) 03811 { 03812 /* remember highest legal value of hiLoc 03813 * -> always update */ 03814 hiLoc = grid->val[0][index[0]]; 03815 } 03816 } 03817 } 03818 03819 ASSERT( fabs(loLoc) < DBL_MAX && fabs(hiLoc) < DBL_MAX && loLoc < hiLoc ); 03820 03821 *loLim = MAX2(*loLim,loLoc); 03822 *hiLim = MIN2(*hiLim,hiLoc); 03823 } 03824 else 03825 { 03826 index[nd] = indlo[nd]; 03827 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim ); 03828 03829 if( indhi[nd] != indlo[nd] ) 03830 { 03831 index[nd] = indhi[nd]; 03832 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim ); 03833 } 03834 } 03835 return; 03836 } 03837 03838 STATIC void InitIndexArrays(stellar_grid *grid, 03839 bool lgList) 03840 { 03841 long i, *index, j, jsize, nd; 03842 double *val; 03843 03844 DEBUG_ENTRY( "InitIndexArrays()" ); 03845 03846 ASSERT( grid->telg != NULL ); 03847 ASSERT( grid->nmods > 0 ); 03848 03849 jsize = 1; 03850 03851 /* this loop creates a list of all unique model parameter values in increasing order */ 03852 for( nd=0; nd < grid->ndim; nd++ ) 03853 { 03854 double pval = grid->telg[0].par[nd]; 03855 grid->val[nd][0] = pval; 03856 grid->nval[nd] = 1; 03857 03858 for( i=1; i < grid->nmods; i++ ) 03859 { 03860 bool lgOutOfRange; 03861 long i1, i2; 03862 03863 pval = grid->telg[i].par[nd]; 03864 FindIndex( grid->val[nd], grid->nval[nd], pval, &i1, &i2, &lgOutOfRange ); 03865 /* if i1 < i2, the new parameter value was not present yet and 03866 * it needs to be inserted in between i1 and i2 --> first move 03867 * all entries from i2 to grid->nval[nd]-1 one slot upward and 03868 * then insert the new value at i2; this also works correctly 03869 * if lgOutOfRange is set, hence no special check is needed */ 03870 if( i1 < i2 ) 03871 { 03872 /* val[nd] has grid->nmods entries, so cannot overflow */ 03873 for( j = grid->nval[nd]-1; j >= i2; j-- ) 03874 grid->val[nd][j+1] = grid->val[nd][j]; 03875 grid->val[nd][i2] = pval; 03876 grid->nval[nd]++; 03877 } 03878 } 03879 03880 jsize *= grid->nval[nd]; 03881 03882 # if DEBUGPRT 03883 printf( "%s[%ld]:", grid->names[nd], grid->nval[nd] ); 03884 for( i=0; i < grid->nval[nd]; i++ ) 03885 printf( " %g", grid->val[nd][i] ); 03886 printf( "\n" ); 03887 # endif 03888 } 03889 03890 index = (long *)MALLOC( sizeof(long)*grid->ndim ); 03891 val = (double *)MALLOC( sizeof(double)*grid->ndim ); 03892 03893 /* this memory will be freed in the calling function */ 03894 grid->jlo = (long *)MALLOC( sizeof(long)*jsize ); 03895 grid->jhi = (long *)MALLOC( sizeof(long)*jsize ); 03896 03897 /* set up square array of model indices; this will be used to 03898 * choose the correct models for the interpolation process */ 03899 FillJ( grid, index, val, grid->ndim, lgList ); 03900 03901 FREE_CHECK( val ); 03902 FREE_CHECK( index ); 03903 03904 if( lgList ) 03905 cdEXIT(EXIT_SUCCESS); 03906 return; 03907 } 03908 03909 STATIC void FillJ(const stellar_grid *grid, 03910 long index[], /* index[grid->ndim] */ 03911 double val[], /* val[grid->ndim] */ 03912 long nd, 03913 bool lgList) 03914 { 03915 DEBUG_ENTRY( "FillJ()" ); 03916 03917 --nd; 03918 03919 if( nd < 0 ) 03920 { 03921 long n = JIndex(grid,index); 03922 SearchModel( grid->telg, grid->nmods, val, grid->ndim, &grid->jlo[n], &grid->jhi[n] ); 03923 } 03924 else 03925 { 03926 for( index[nd]=0; index[nd] < grid->nval[nd]; index[nd]++ ) 03927 { 03928 val[nd] = grid->val[nd][index[nd]]; 03929 FillJ( grid, index, val, nd, lgList ); 03930 } 03931 } 03932 03933 if( lgList && nd == MIN2(grid->ndim-1,1) ) 03934 { 03935 long n; 03936 bool lgAge = ( strcmp( grid->names[0], "Age" ) == 0 ); 03937 03938 fprintf( ioQQQ, "\n" ); 03939 if( grid->ndim > 2 ) 03940 { 03941 fprintf( ioQQQ, "subgrid for" ); 03942 for( n = nd+1; n < grid->ndim; n++ ) 03943 fprintf( ioQQQ, " %s=%g", grid->names[n], val[n] ); 03944 fprintf( ioQQQ, ":\n\n" ); 03945 } 03946 if( grid->ndim > 1 ) 03947 { 03948 fprintf( ioQQQ, "Teff\\lg g|" ); 03949 for( n = 0; n < grid->nval[1]; n++ ) 03950 fprintf( ioQQQ, " %5.2f", grid->val[1][n] ); 03951 fprintf( ioQQQ, "\n" ); 03952 fprintf( ioQQQ, "---------|" ); 03953 for( n = 0; n < grid->nval[1]; n++ ) 03954 fprintf( ioQQQ, "------" ); 03955 } 03956 else 03957 { 03958 if( lgAge ) 03959 fprintf( ioQQQ, " Age |\n" ); 03960 else 03961 fprintf( ioQQQ, " Teff |\n" ); 03962 fprintf( ioQQQ, "---------|------" ); 03963 } 03964 fprintf( ioQQQ, "\n" ); 03965 for( index[0]=0; index[0] < grid->nval[0]; index[0]++ ) 03966 { 03967 if( lgAge ) 03968 fprintf( ioQQQ, "%8.2e |", grid->val[0][index[0]] ); 03969 else 03970 fprintf( ioQQQ, "%8.0f |", grid->val[0][index[0]] ); 03971 if( grid->ndim > 1 ) 03972 { 03973 for( index[1]=0; index[1] < grid->nval[1]; index[1]++ ) 03974 if( grid->jlo[JIndex(grid,index)] == grid->jhi[JIndex(grid,index)] && 03975 grid->jlo[JIndex(grid,index)] >= 0 ) 03976 fprintf( ioQQQ, " %5ld", grid->jlo[JIndex(grid,index)]+1 ); 03977 else 03978 fprintf( ioQQQ, " --" ); 03979 } 03980 else 03981 { 03982 fprintf( ioQQQ, " %5ld", grid->jlo[JIndex(grid,index)]+1 ); 03983 } 03984 fprintf( ioQQQ, "\n" ); 03985 } 03986 } 03987 return; 03988 } 03989 03990 STATIC long JIndex(const stellar_grid *grid, 03991 const long index[]) /* index[grid->ndim] */ 03992 { 03993 long i, ind, mul; 03994 03995 DEBUG_ENTRY( "JIndex()" ); 03996 03997 ind = 0; 03998 mul = 1; 03999 for( i=0; i < grid->ndim; i++ ) 04000 { 04001 ind += index[i]*mul; 04002 mul *= grid->nval[i]; 04003 } 04004 return ind; 04005 } 04006 04007 STATIC void SearchModel(const mpp telg[], /* telg[nmods] */ 04008 long nmods, 04009 const double val[], /* val[ndim] */ 04010 long ndim, 04011 long *index_low, 04012 long *index_high) 04013 { 04014 long i, nd; 04015 double alogg_low = -DBL_MAX, alogg_high = DBL_MAX; 04016 04017 DEBUG_ENTRY( "SearchModel()" ); 04018 04019 /* given values for the model parameters, this routine searches for 04020 * the atmosphere model that is the best match. If all parameters 04021 * can be matched simultaneously the choice is obvious, but this 04022 * cannot always be achieved (typically for high Teff, the low 04023 * log(g) models will be missing). The rule is that all parameters 04024 * except log(g) must always be matched (such a model is not always 04025 * guaranteed to exist). If all requested parameters can be matched 04026 * exactly, both index_low and index_high will point to that model. 04027 * If all parameters except log(g) can be matched exactly, it will 04028 * return the model with the lowest log(g) value larger than the 04029 * requested value in index_high, and the model with the highest 04030 * log(g) value lower than the requested value in index_low. If 04031 * either requirement cannot be fulfilled, -2 will be returned. */ 04032 04033 *index_low = *index_high = -2; 04034 for( i=0; i < nmods; i++ ) 04035 { 04036 bool lgNext = false; 04037 /* ignore models with different parameters */ 04038 for( nd=0; nd < ndim; nd++ ) 04039 { 04040 if( nd != 1 && fabs(telg[i].par[nd]-val[nd]) > 10.*DBL_EPSILON*fabs(val[nd]) ) 04041 { 04042 lgNext = true; 04043 break; 04044 } 04045 } 04046 if( lgNext ) 04047 continue; 04048 04049 /* an exact match is found */ 04050 /* this needs to be <= otherwise the test will fail for values equal to 0 */ 04051 if( ndim == 1 || fabs(telg[i].par[1]-val[1]) <= 10.*DBL_EPSILON*fabs(val[1]) ) 04052 { 04053 *index_low = i; 04054 *index_high = i; 04055 return; 04056 } 04057 /* keep a record of the highest log(g) model smaller than alogg */ 04058 if( telg[i].par[1] < val[1] && telg[i].par[1] > alogg_low ) 04059 { 04060 *index_low = i; 04061 alogg_low = telg[i].par[1]; 04062 } 04063 /* also keep a record of the lowest log(g) model greater than alogg */ 04064 if( telg[i].par[1] > val[1] && telg[i].par[1] < alogg_high ) 04065 { 04066 *index_high = i; 04067 alogg_high = telg[i].par[1]; 04068 } 04069 } 04070 return; 04071 } 04072 04073 STATIC void FindIndex(const double xval[], /* xval[NVAL] */ 04074 long NVAL, 04075 double x, 04076 long *ind1, 04077 long *ind2, 04078 bool *lgInvalid) 04079 { 04080 bool lgOutLo, lgOutHi; 04081 long i; 04082 04083 DEBUG_ENTRY( "FindIndex()" ); 04084 04085 /* this routine searches for indices ind1, ind2 such that 04086 * xval[ind1] < x < xval[ind2] 04087 * if x is equal to one of the values in xval, then 04088 * ind1 == ind2 and xval[ind1] == x 04089 * 04090 * if x is outside the range xval[0] ... xval[NVAL-1] 04091 * then lgInvalid will be set to true 04092 * 04093 * NB NB -- this routine implicitly assumes that xval is 04094 * strictly monotonically increasing! 04095 */ 04096 04097 ASSERT( NVAL > 0 ); 04098 04099 /* is x outside of range xval[0] ... xval[NVAL-1]? */ 04100 lgOutLo = ( x-xval[0] < -10.*DBL_EPSILON*fabs(xval[0]) ); 04101 lgOutHi = ( x-xval[NVAL-1] > 10.*DBL_EPSILON*fabs(xval[NVAL-1]) ); 04102 04103 if( lgOutLo || lgOutHi ) 04104 { 04105 /* pretend there are two fictitious array elements 04106 * xval[-1] = -Inf and xval[NVAL] = +Inf, 04107 * and return ind1 and ind2 accordingly. This behavior 04108 * is needed for InitIndexArrays() to work correctly */ 04109 *ind1 = lgOutLo ? -1 : NVAL-1; 04110 *ind2 = lgOutLo ? 0 : NVAL; 04111 *lgInvalid = true; 04112 return; 04113 } 04114 04115 *lgInvalid = false; 04116 04117 /* there are more efficient ways of doing this, e.g. a binary search. 04118 * However, the xval arrays typically only have 1 or 2 dozen elements, 04119 * so the overhead is negligible and the clarity of this code is preferred */ 04120 04121 /* first look for an "exact" match */ 04122 for( i=0; i < NVAL; i++ ) 04123 { 04124 /* this needs to be <= otherwise the test will fail for x == 0 */ 04125 if( fabs(xval[i]-x) <= 10.*DBL_EPSILON*fabs(xval[i]) ) 04126 { 04127 *ind1 = i; 04128 *ind2 = i; 04129 return; 04130 } 04131 } 04132 04133 /* no match was found -> bracket the x value */ 04134 for( i=0; i < NVAL-1; i++ ) 04135 { 04136 if( xval[i] < x && x < xval[i+1] ) 04137 { 04138 *ind1 = i; 04139 *ind2 = i+1; 04140 return; 04141 } 04142 } 04143 04144 /* this should never be reached ! */ 04145 fprintf( ioQQQ, " insanity in FindIndex\n" ); 04146 ShowMe(); 04147 cdEXIT(EXIT_FAILURE); 04148 } 04149 04150 STATIC bool lgFileReadable(const char *chFnam, process_counter& pc, access_scheme scheme) 04151 { 04152 DEBUG_ENTRY( "lgFileReadable()" ); 04153 04154 FILE *ioIN; 04155 04156 ioIN = open_data( chFnam, "r", scheme ); 04157 if( ioIN != NULL ) 04158 { 04159 fclose( ioIN ); 04160 ++pc.nFound; 04161 return true; 04162 } 04163 else 04164 { 04165 return false; 04166 } 04167 } 04168 04169 /*ValidateGrid: check each model in the grid to see if it has the correct Teff */ 04170 STATIC void ValidateGrid(const stellar_grid *grid, 04171 double toler) 04172 { 04173 long i, k, nd; 04174 realnum *anu, *flux; 04175 04176 DEBUG_ENTRY( "ValidateGrid()" ); 04177 04178 if( strcmp( grid->names[0], "Teff" ) != 0 ) 04179 { 04180 return; 04181 } 04182 04183 anu = (realnum *)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) ); 04184 flux = (realnum *)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) ); 04185 04186 GetModel( grid, -1, anu, lgSILENT, lgLINEAR ); 04187 04188 for( i=0; i < grid->nmods; i++ ) 04189 { 04190 fprintf( ioQQQ, "testing model %ld ", i+1 ); 04191 for( nd=0; nd < grid->npar; nd++ ) 04192 fprintf( ioQQQ, " %s %g", grid->names[nd], grid->telg[i].par[nd] ); 04193 04194 GetModel( grid, i, flux, lgSILENT, lgLINEAR ); 04195 04196 if( lgValidModel( anu, flux, grid->telg[i].par[0], toler ) ) 04197 fprintf( ioQQQ, " OK\n" ); 04198 04199 if( false ) 04200 { 04201 FILE *ioBUG = fopen( "atmosphere_dump.txt", ( i == 0 ) ? "w" : "a" ); 04202 04203 fprintf( ioBUG, "######## MODEL %ld", i+1 ); 04204 for( nd=0; nd < grid->npar; nd++ ) 04205 fprintf( ioBUG, " %s %g", grid->names[nd], grid->telg[i].par[nd] ); 04206 fprintf( ioBUG, " ####################\n" ); 04207 04208 for( k=0; k < rfield.nupper; ++k ) 04209 fprintf( ioBUG, "%e %e\n", anu[k], flux[k] ); 04210 04211 fclose( ioBUG ); 04212 } 04213 } 04214 04215 FREE_CHECK( flux ); 04216 FREE_CHECK( anu ); 04217 return; 04218 } 04219 04220 STATIC bool lgValidModel(const realnum anu[], 04221 const realnum flux[], 04222 double Teff, 04223 double toler) 04224 { 04225 bool lgPassed = true; 04226 long k; 04227 double chk, lumi; 04228 04229 DEBUG_ENTRY( "lgValidModel()" ); 04230 04231 ASSERT( Teff > 0. ); 04232 04233 lumi = 0.; 04234 /* rebinned models are in cgs F_nu units */ 04235 for( k=1; k < rfield.nupper; k++ ) 04236 lumi += (anu[k] - anu[k-1])*(flux[k] + flux[k-1])/2.; 04237 04238 /* now convert luminosity to effective temperature */ 04239 chk = pow(lumi*FR1RYD/STEFAN_BOLTZ,0.25); 04240 /* the allowed tolerance is set by the caller in toler */ 04241 if( fabs(Teff - chk) > toler*Teff ) { 04242 fprintf( ioQQQ, "\n*** WARNING, Teff discrepancy for this model, expected Teff %.2f, ", Teff); 04243 fprintf( ioQQQ, "integration yielded Teff %.2f, delta %.2f%%\n", chk, (chk/Teff-1.)*100. ); 04244 lgPassed = false; 04245 } 04246 return lgPassed; 04247 } 04248 04249 /*RebinAtmosphere: generic routine for rebinning atmospheres onto Cloudy grid */ 04250 STATIC void RebinAtmosphere(long nCont, /* the number of points in the incident continuum*/ 04251 const realnum StarEner[], /* StarEner[nCont], the freq grid for the model, in Ryd*/ 04252 const realnum StarFlux[], /* StarFlux[nCont], the original model flux */ 04253 realnum CloudyFlux[], /* CloudyFlux[NC_ELL], the model flux on the cloudy grid */ 04254 long nEdge, /* the number of bound-free continuum edges in AbsorbEdge */ 04255 const realnum AbsorbEdge[]) /* AbsorbEdge[nEdge], energies of the edges */ 04256 { 04257 bool lgDone; 04258 long int ind, 04259 j, 04260 k; 04261 /* >>chng 00 jun 02, demoted next two to realnum, PvH */ 04262 realnum BinHigh, 04263 BinLow, 04264 BinMid, 04265 BinNext, 04266 *EdgeLow=NULL, 04267 *EdgeHigh=NULL, 04268 *StarPower; 04269 04270 DEBUG_ENTRY( "RebinAtmosphere()" ); 04271 04272 if( nEdge > 0 ) 04273 { 04274 EdgeLow = (realnum*)MALLOC( sizeof(realnum)*(unsigned)nEdge ); 04275 EdgeHigh = (realnum*)MALLOC( sizeof(realnum)*(unsigned)nEdge ); 04276 } 04277 04278 /* this loop should be before the next loop, otherwise models with a 04279 * very strong He II edge (i.e. no flux beyond that edge) will fail */ 04280 for( j=0; j < nEdge; j++ ) 04281 { 04282 ind = RebinFind(StarEner,nCont,AbsorbEdge[j]); 04283 04284 /* sanity check */ 04285 ASSERT( ind >= 0 && ind+1 < nCont ); 04286 04287 EdgeLow[j] = StarEner[ind]; 04288 EdgeHigh[j] = StarEner[ind+1]; 04289 } 04290 04291 /* cut off that part of the Wien tail that evaluated to zero */ 04292 /* >> chng 05 nov 22, inverted loop, slightly faster PvH */ 04293 /*for( j=nCont-1; j >= 0; j-- )*/ 04294 for( j=0; j < nCont; j++ ) 04295 { 04296 if( StarFlux[j] == 0.f ) 04297 { 04298 nCont = j; 04299 break; 04300 } 04301 } 04302 ASSERT( nCont > 0 ); 04303 04304 StarPower = (realnum *)MALLOC( sizeof(realnum)*(unsigned)(nCont-1) ); 04305 04306 for( j=0; j < nCont-1; j++ ) 04307 { 04308 double ratio_x, ratio_y; 04309 04310 /* >>chng 05 nov 22, add sanity check to prevent invalid fp operations */ 04311 ASSERT( StarEner[j+1] > StarEner[j] ); 04312 04313 /* >>chng 06 aug 11, on some systems (e.g., macbook pro) y/x can get evaluated as y*(1/x); 04314 * this causes overflows if x is a denormalized number, hence we force a cast to double, PvH */ 04315 ratio_x = (double)StarEner[j+1]/(double)StarEner[j]; 04316 ratio_y = (double)StarFlux[j+1]/(double)StarFlux[j]; 04317 StarPower[j] = (realnum)(log(ratio_y)/log(ratio_x)); 04318 } 04319 04320 for( j=0; j < rfield.nupper; j++ ) 04321 { 04322 /* >>chng 05 nov 22, modified BinLow, BinHigh, BinNext to make boundaries match exactly, PvH */ 04323 /* BinLow is lower bound of this continuum cell */ 04324 BinLow = ( j > 0 ) ? 04325 (realnum)sqrt(rfield.anu[j-1]*rfield.anu[j]) : (realnum)sqrt(POW3(rfield.anu[0])/rfield.anu[1]); 04326 04327 /* BinHigh is upper bound of this continuum cell */ 04328 BinHigh = ( j+1 < rfield.nupper ) ? 04329 (realnum)sqrt(rfield.anu[j]*rfield.anu[j+1]) : rfield.anu[rfield.nupper-1]; 04330 04331 /* BinNext is upper bound of next continuum cell */ 04332 BinNext = ( j+2 < rfield.nupper ) ? 04333 (realnum)sqrt(rfield.anu[j+1]*rfield.anu[j+2]) : rfield.anu[rfield.nupper-1]; 04334 04335 lgDone = false; 04336 04337 /* >>chng 00 aug 14, take special care not to interpolate over major edges, 04338 * the region in between EdgeLow and EdgeHigh should be avoided, 04339 * the spectrum is extremely steep there, leading to significant roundoff error, PvH */ 04340 for( k=0; k < nEdge; k++ ) 04341 { 04342 if( BinLow < EdgeLow[k] && BinNext > EdgeHigh[k] ) 04343 { 04344 BinMid = 0.99999f*EdgeLow[k]; 04345 CloudyFlux[j] = RebinSingleCell(BinLow,BinMid,StarEner,StarFlux,StarPower,nCont); 04346 j++; 04347 04348 /* sanity check */ 04349 ASSERT( j < rfield.nupper ); 04350 04351 BinMid = 1.00001f*EdgeHigh[k]; 04352 CloudyFlux[j] = RebinSingleCell(BinMid,BinNext,StarEner,StarFlux,StarPower,nCont); 04353 lgDone = true; 04354 break; 04355 } 04356 } 04357 04358 /* default case when we are not close to an edge */ 04359 if( !lgDone ) 04360 { 04361 CloudyFlux[j] = RebinSingleCell(BinLow,BinHigh,StarEner,StarFlux,StarPower,nCont); 04362 } 04363 } 04364 04365 FREE_CHECK( StarPower ); 04366 FREE_SAFE( EdgeHigh ); 04367 FREE_SAFE( EdgeLow ); 04368 return; 04369 } 04370 04371 STATIC realnum RebinSingleCell(realnum BinLow, 04372 realnum BinHigh, 04373 const realnum StarEner[], /* StarEner[nCont] */ 04374 const realnum StarFlux[], /* StarFlux[nCont] */ 04375 const realnum StarPower[], /* StarPower[nCont-1] */ 04376 long nCont) 04377 { 04378 long int i, 04379 ipHi, 04380 ipLo; 04381 double anu, 04382 retval, 04383 widflx; 04384 double sum, 04385 v1, 04386 val, 04387 x1, 04388 x2; 04389 04390 DEBUG_ENTRY( "RebinSingleCell()" ); 04391 04392 /* >>chng 05 nov 22, use geometric mean instead of arithmetic mean, PvH */ 04393 anu = sqrt(BinLow*BinHigh); 04394 /* >>chng 05 nov 22, reduce widflx if cell sticks out above highest frequency in model, PvH */ 04395 widflx = MIN2(BinHigh,StarEner[nCont-1])-BinLow; 04396 04397 if( BinLow < StarEner[0] ) 04398 { 04399 /* this is case where Cloudy's continuum is below stellar continuum, 04400 * (at least for part of the cell), so we do Rayleigh Jeans extrapolation */ 04401 retval = (realnum)(StarFlux[0]*pow(anu/StarEner[0],2.)); 04402 } 04403 else if( BinLow > StarEner[nCont-1] ) 04404 { 04405 /* case where cloudy continuum is entirely above highest stellar point */ 04406 retval = 0.0e00; 04407 } 04408 else 04409 { 04410 /* now go through stellar continuum to find bins corresponding to 04411 * this cloudy cell, stellar continuum defined through nCont cells */ 04412 ipLo = RebinFind(StarEner,nCont,BinLow); 04413 ipHi = RebinFind(StarEner,nCont,BinHigh); 04414 /* sanity check */ 04415 ASSERT( ipLo >= 0 && ipLo < nCont-1 && ipHi >= ipLo ); 04416 04417 if( ipHi == ipLo ) 04418 { 04419 /* Do the case where the cloudy cell and its edges are between 04420 * two adjacent stellar model points: do power-law interpolation */ 04421 retval = (realnum)(StarFlux[ipLo]*pow(anu/StarEner[ipLo],(double)StarPower[ipLo])); 04422 } 04423 else 04424 { 04425 /* Do the case where the cloudy cell and its edges span two or more 04426 * stellar model cells: add segments with power-law interpolation up to 04427 * do the averaging.*/ 04428 04429 sum = 0.; 04430 04431 /* ipHi points to stellar point at high end of cloudy continuum cell, 04432 * if the Cloudy cell extends beyond the stellar grid, ipHi == nCont-1 04433 * and the MIN2(ipHi,nCont-2) prevents access beyond allocated memory 04434 * ipLo points to low end, above we asserted that 0 <= ipLo < nCont-1 */ 04435 for( i=ipLo; i <= MIN2(ipHi,nCont-2); i++ ) 04436 { 04437 double pp1 = StarPower[i] + 1.; 04438 04439 if( i == ipLo ) 04440 { 04441 x1 = BinLow; 04442 x2 = StarEner[i+1]; 04443 v1 = StarFlux[i]*pow(x1/StarEner[i],(double)StarPower[i]); 04444 /*v2 = StarFlux[i+1];*/ 04445 } 04446 04447 else if( i == ipHi ) 04448 { 04449 x2 = BinHigh; 04450 x1 = StarEner[i]; 04451 /*v2 = StarFlux[i]*pow(x2/StarEner[i],StarPower[i]);*/ 04452 v1 = StarFlux[i]; 04453 } 04454 04455 /*if( i > ipLo && i < ipHi )*/ 04456 else 04457 { 04458 x1 = StarEner[i]; 04459 x2 = StarEner[i+1]; 04460 v1 = StarFlux[i]; 04461 /*v2 = StarFlux[i+1];*/ 04462 } 04463 04464 if( fabs(pp1) < 0.001 ) 04465 { 04466 val = x1*v1*log(x2/x1); 04467 } 04468 else 04469 { 04470 val = pow(x2/x1,pp1) - 1.; 04471 val = val*x1*v1/pp1; 04472 } 04473 sum += val; 04474 } 04475 04476 retval = sum/widflx; 04477 } 04478 } 04479 return (realnum)retval; 04480 } 04481 04482 STATIC long RebinFind(const realnum array[], /* array[nArr] */ 04483 long nArr, 04484 realnum val) 04485 { 04486 long i1, 04487 i2, 04488 i3, 04489 ind = -2, 04490 sgn; 04491 04492 DEBUG_ENTRY( "RebinFind()" ); 04493 04494 /* sanity check */ 04495 ASSERT( nArr > 1 ); 04496 04497 /* return ind(val) such that array[ind] <= val <= array[ind+1], 04498 * 04499 * NB NB: this routine assumes that array[] increases monotonically ! 04500 * 04501 * the first two clauses indicate out-of-bounds conditions and 04502 * guarantee that when val1 <= val2, also ind(val1) <= ind(val2) */ 04503 04504 if( val < array[0] ) 04505 { 04506 ind = -1; 04507 } 04508 else if( val > array[nArr-1] ) 04509 { 04510 ind = nArr-1; 04511 } 04512 else 04513 { 04514 /* do a binary search for ind */ 04515 i1 = 0; 04516 i3 = nArr-1; 04517 while( i3-i1 > 1 ) 04518 { 04519 i2 = (i1+i3)/2; 04520 sgn = sign3(val-array[i2]); 04521 04522 switch(sgn) 04523 { 04524 case -1: 04525 i3 = i2; 04526 break; 04527 case 0: 04528 ind = i2; 04529 return( ind ); 04530 case 1: 04531 i1 = i2; 04532 break; 04533 } 04534 } 04535 ind = i1; 04536 } 04537 04538 /* sanity check */ 04539 ASSERT( ind > -2 ); 04540 return ind; 04541 } 04542 /*lint +e785 too few initializers */ 04543 /*lint +e801 use of go to depreciated */