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