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