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