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