17static const int MNTS = 200;
37#define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
38#define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
157 long,
const double[],
int);
158STATIC void RauchReadMPP(vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&);
169 const long[],
long[],
long,vector<realnum>&,
IntStage);
171 const long[],
long[],
long,
long,vector<realnum>&);
175 const realnum[],
double*,
double*);
213 fprintf(
ioQQQ,
"\n I will now list all stellar atmosphere grids that are ready to be used (if any).\n" );
214 fprintf(
ioQQQ,
" User-defined stellar atmosphere grids will not be included in this list.\n\n" );
223 fprintf(
ioQQQ,
" table star atlas Z+1.0 <Teff> [ <log(g)> ]\n" );
225 fprintf(
ioQQQ,
" table star atlas Z+0.5 <Teff> [ <log(g)> ]\n" );
227 fprintf(
ioQQQ,
" table star atlas Z+0.3 <Teff> [ <log(g)> ]\n" );
229 fprintf(
ioQQQ,
" table star atlas Z+0.2 <Teff> [ <log(g)> ]\n" );
231 fprintf(
ioQQQ,
" table star atlas Z+0.1 <Teff> [ <log(g)> ]\n" );
233 fprintf(
ioQQQ,
" table star atlas Z+0.0 <Teff> [ <log(g)> ]\n" );
235 fprintf(
ioQQQ,
" table star atlas Z-0.1 <Teff> [ <log(g)> ]\n" );
237 fprintf(
ioQQQ,
" table star atlas Z-0.2 <Teff> [ <log(g)> ]\n" );
239 fprintf(
ioQQQ,
" table star atlas Z-0.3 <Teff> [ <log(g)> ]\n" );
241 fprintf(
ioQQQ,
" table star atlas Z-0.5 <Teff> [ <log(g)> ]\n" );
243 fprintf(
ioQQQ,
" table star atlas Z-1.0 <Teff> [ <log(g)> ]\n" );
245 fprintf(
ioQQQ,
" table star atlas Z-1.5 <Teff> [ <log(g)> ]\n" );
247 fprintf(
ioQQQ,
" table star atlas Z-2.0 <Teff> [ <log(g)> ]\n" );
249 fprintf(
ioQQQ,
" table star atlas Z-2.5 <Teff> [ <log(g)> ]\n" );
251 fprintf(
ioQQQ,
" table star atlas Z-3.0 <Teff> [ <log(g)> ]\n" );
253 fprintf(
ioQQQ,
" table star atlas Z-3.5 <Teff> [ <log(g)> ]\n" );
255 fprintf(
ioQQQ,
" table star atlas Z-4.0 <Teff> [ <log(g)> ]\n" );
257 fprintf(
ioQQQ,
" table star atlas Z-4.5 <Teff> [ <log(g)> ]\n" );
259 fprintf(
ioQQQ,
" table star atlas Z-5.0 <Teff> [ <log(g)> ]\n" );
262 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.5 <Teff> [ <log(g)> ]\n" );
264 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.2 <Teff> [ <log(g)> ]\n" );
266 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.0 <Teff> [ <log(g)> ]\n" );
268 fprintf(
ioQQQ,
" table star atlas odfnew Z-0.5 <Teff> [ <log(g)> ]\n" );
270 fprintf(
ioQQQ,
" table star atlas odfnew Z-1.0 <Teff> [ <log(g)> ]\n" );
272 fprintf(
ioQQQ,
" table star atlas odfnew Z-1.5 <Teff> [ <log(g)> ]\n" );
274 fprintf(
ioQQQ,
" table star atlas odfnew Z-2.0 <Teff> [ <log(g)> ]\n" );
276 fprintf(
ioQQQ,
" table star atlas odfnew Z-2.5 <Teff> [ <log(g)> ]\n" );
279 fprintf(
ioQQQ,
" table star atlas 3-dim <Teff> <log(g)> <log(Z)>\n" );
282 fprintf(
ioQQQ,
" table star atlas odfnew 3-dim <Teff> <log(g)> <log(Z)>\n" );
285 fprintf(
ioQQQ,
" table star costar solar (see Hazy for parameters)\n" );
287 fprintf(
ioQQQ,
" table star costar halo (see Hazy for parameters)\n" );
290 fprintf(
ioQQQ,
" table star kurucz79 <Teff>\n" );
293 fprintf(
ioQQQ,
" table star mihalas <Teff>\n" );
296 fprintf(
ioQQQ,
" table star rauch H-Ca solar <Teff> [ <log(g)> ]\n" );
298 fprintf(
ioQQQ,
" table star rauch H-Ca halo <Teff> [ <log(g)> ]\n" );
300 fprintf(
ioQQQ,
" table star rauch H-Ca 3-dim <Teff> <log(g)> <log(Z)>\n" );
303 fprintf(
ioQQQ,
" table star rauch H-Ni solar <Teff> [ <log(g)> ]\n" );
305 fprintf(
ioQQQ,
" table star rauch H-Ni halo <Teff> [ <log(g)> ]\n" );
307 fprintf(
ioQQQ,
" table star rauch H-Ni 3-dim <Teff> <log(g)> <log(Z)>\n" );
310 fprintf(
ioQQQ,
" table star rauch pg1159 <Teff> [ <log(g)> ]\n" );
312 fprintf(
ioQQQ,
" table star rauch co wd <Teff>\n" );
315 fprintf(
ioQQQ,
" table star rauch hydrogen <Teff> [ <log(g)> ]\n" );
318 fprintf(
ioQQQ,
" table star rauch helium <Teff> [ <log(g)> ]\n" );
321 fprintf(
ioQQQ,
" table star rauch H+He <Teff> <log(g)> <frac(He)>\n" );
324 fprintf(
ioQQQ,
" table star \"starburst99.mod\" <age>\n" );
326 fprintf(
ioQQQ,
" table star \"starburst99_2d.mod\" <age> <Z>\n" );
329 fprintf(
ioQQQ,
" table star tlusty OBstar Z+0.3 <Teff> [ <log(g)> ]\n" );
331 fprintf(
ioQQQ,
" table star tlusty OBstar Z+0.0 <Teff> [ <log(g)> ]\n" );
333 fprintf(
ioQQQ,
" table star tlusty OBstar Z-0.3 <Teff> [ <log(g)> ]\n" );
335 fprintf(
ioQQQ,
" table star tlusty OBstar Z-0.7 <Teff> [ <log(g)> ]\n" );
337 fprintf(
ioQQQ,
" table star tlusty OBstar Z-1.0 <Teff> [ <log(g)> ]\n" );
339 fprintf(
ioQQQ,
" table star tlusty OBstar Z-inf <Teff> [ <log(g)> ]\n" );
342 fprintf(
ioQQQ,
" table star tlusty OBstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
345 fprintf(
ioQQQ,
" table star tlusty Bstar Z+0.3 <Teff> [ <log(g)> ]\n" );
347 fprintf(
ioQQQ,
" table star tlusty Bstar Z+0.0 <Teff> [ <log(g)> ]\n" );
349 fprintf(
ioQQQ,
" table star tlusty Bstar Z-0.3 <Teff> [ <log(g)> ]\n" );
351 fprintf(
ioQQQ,
" table star tlusty Bstar Z-0.7 <Teff> [ <log(g)> ]\n" );
353 fprintf(
ioQQQ,
" table star tlusty Bstar Z-1.0 <Teff> [ <log(g)> ]\n" );
355 fprintf(
ioQQQ,
" table star tlusty Bstar Z-inf <Teff> [ <log(g)> ]\n" );
358 fprintf(
ioQQQ,
" table star tlusty Bstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
361 fprintf(
ioQQQ,
" table star tlusty Ostar Z+0.3 <Teff> [ <log(g)> ]\n" );
363 fprintf(
ioQQQ,
" table star tlusty Ostar Z+0.0 <Teff> [ <log(g)> ]\n" );
365 fprintf(
ioQQQ,
" table star tlusty Ostar Z-0.3 <Teff> [ <log(g)> ]\n" );
367 fprintf(
ioQQQ,
" table star tlusty Ostar Z-0.7 <Teff> [ <log(g)> ]\n" );
369 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.0 <Teff> [ <log(g)> ]\n" );
371 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.5 <Teff> [ <log(g)> ]\n" );
373 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.7 <Teff> [ <log(g)> ]\n" );
375 fprintf(
ioQQQ,
" table star tlusty Ostar Z-2.0 <Teff> [ <log(g)> ]\n" );
377 fprintf(
ioQQQ,
" table star tlusty Ostar Z-3.0 <Teff> [ <log(g)> ]\n" );
379 fprintf(
ioQQQ,
" table star tlusty Ostar Z-inf <Teff> [ <log(g)> ]\n" );
382 fprintf(
ioQQQ,
" table star tlusty Ostar 3-dim <Teff> <log(g)> <log(Z)>\n" );
385 fprintf(
ioQQQ,
" table star werner <Teff> [ <log(g)> ]\n" );
388 fprintf(
ioQQQ,
" table star wmbasic <Teff> <log(g)> <log(Z)>\n" );
411 fprintf(
ioQQQ,
" AtlasCompile on the job.\n" );
433 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp10k2.ascii",
"atlas_fp10k2.mod", Edges, 3L, pc );
435 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp05k2.ascii",
"atlas_fp05k2.mod", Edges, 3L, pc );
437 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp03k2.ascii",
"atlas_fp03k2.mod", Edges, 3L, pc );
439 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp02k2.ascii",
"atlas_fp02k2.mod", Edges, 3L, pc );
441 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp01k2.ascii",
"atlas_fp01k2.mod", Edges, 3L, pc );
443 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp00k2.ascii",
"atlas_fp00k2.mod", Edges, 3L, pc );
445 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm01k2.ascii",
"atlas_fm01k2.mod", Edges, 3L, pc );
447 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm02k2.ascii",
"atlas_fm02k2.mod", Edges, 3L, pc );
449 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm03k2.ascii",
"atlas_fm03k2.mod", Edges, 3L, pc );
451 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm05k2.ascii",
"atlas_fm05k2.mod", Edges, 3L, pc );
453 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm10k2.ascii",
"atlas_fm10k2.mod", Edges, 3L, pc );
455 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm15k2.ascii",
"atlas_fm15k2.mod", Edges, 3L, pc );
457 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm20k2.ascii",
"atlas_fm20k2.mod", Edges, 3L, pc );
459 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm25k2.ascii",
"atlas_fm25k2.mod", Edges, 3L, pc );
461 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm30k2.ascii",
"atlas_fm30k2.mod", Edges, 3L, pc );
463 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm35k2.ascii",
"atlas_fm35k2.mod", Edges, 3L, pc );
465 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm40k2.ascii",
"atlas_fm40k2.mod", Edges, 3L, pc );
467 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm45k2.ascii",
"atlas_fm45k2.mod", Edges, 3L, pc );
469 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm50k2.ascii",
"atlas_fm50k2.mod", Edges, 3L, pc );
474 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp05k2_odfnew.ascii",
"atlas_fp05k2_odfnew.mod",
478 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp02k2_odfnew.ascii",
"atlas_fp02k2_odfnew.mod",
482 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp00k2_odfnew.ascii",
"atlas_fp00k2_odfnew.mod",
486 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm05k2_odfnew.ascii",
"atlas_fm05k2_odfnew.mod",
490 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm10k2_odfnew.ascii",
"atlas_fm10k2_odfnew.mod",
494 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm15k2_odfnew.ascii",
"atlas_fm15k2_odfnew.mod",
498 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm20k2_odfnew.ascii",
"atlas_fm20k2_odfnew.mod",
502 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm25k2_odfnew.ascii",
"atlas_fm25k2_odfnew.mod",
510 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_3d_odfnew.ascii",
"atlas_3d_odfnew.mod", Edges, 3L, pc );
518 const char *chMetalicity,
519 const char *chODFNew,
529 grid.name =
"atlas_";
535 grid.name += chMetalicity;
538 grid.name += chODFNew;
545 strcpy( chIdent,
"3-dim" );
549 strcpy( chIdent,
"Z " );
550 strcat( chIdent, chMetalicity );
552 strcat( chIdent, ( strlen(chODFNew) == 0 ?
" Kurucz" :
" ODFNew" ) );
553 grid.ident = chIdent;
555 grid.command =
"COMPILE STARS";
593 fprintf(
ioQQQ,
" CoStarCompile on the job.\n" );
636 grid.name = ( lgHalo ?
"Sc1_costar_halo.mod" :
"Sc1_costar_solar.mod" );
640 grid.ident =
" costar";
642 grid.command =
"COMPILE STARS";
697 string OutName( InName );
701 fprintf(
ioQQQ,
" GridCompile on the job.\n" );
704 string::size_type ptr = OutName.find(
'.' );
705 ASSERT( ptr != string::npos );
706 OutName.replace( ptr, string::npos,
".mod" );
718 grid.ident =
"bogus ident.";
719 grid.command =
"bogus command.";
725 if( strcmp(
grid.names[0],
"Teff" ) == 0 )
727 fprintf(
ioQQQ,
" GridCompile: checking effective temperatures...\n" );
741 const char *FileName,
752 string chTruncName( FileName );
753 string::size_type ptr = chTruncName.find(
'.' );
754 if( ptr != string::npos )
755 chTruncName.replace( ptr, string::npos,
"" );
757 grid.name = FileName;
761 sprintf( chIdent,
"%12.12s", chTruncName.c_str() );
762 grid.ident = chIdent;
764 string chString(
"COMPILE STARS \"" + chTruncName +
".ascii\"" );
765 grid.command = chString.c_str();
786 fprintf(
ioQQQ,
" Kurucz79Compile on the job.\n" );
810 grid.name =
"kurucz79.mod";
814 grid.ident =
" Kurucz 1979";
816 grid.command =
"COMPILE STARS";
837 fprintf(
ioQQQ,
" MihalasCompile on the job.\n" );
860 grid.name =
"mihalas.mod";
864 grid.ident =
" Mihalas";
866 grid.command =
"COMPILE STARS";
902 static const double par2[2] = { 0., -1. };
905 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 };
909 fprintf(
ioQQQ,
" RauchCompile on the job.\n" );
911 RauchReadMPP( telg1, telg2, telg3, telg4, telg5, telg6 );
919 fprintf(
ioQQQ,
" Creating rauch_h-ca_solar.ascii....\n" );
926 fprintf(
ioQQQ,
" Creating rauch_h-ca_halo.ascii....\n" );
935 fprintf(
ioQQQ,
" Creating rauch_h-ca_3d.ascii....\n" );
946 fprintf(
ioQQQ,
" Creating rauch_h-ni_solar.ascii....\n" );
947 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h-ni_solar.ascii",
"_solar_iron.bin_0.1",
954 fprintf(
ioQQQ,
" Creating rauch_h-ni_halo.ascii....\n" );
955 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h-ni_halo.ascii",
"_halo__iron.bin_0.1",
963 fprintf(
ioQQQ,
" Creating rauch_h-ni_3d.ascii....\n" );
971 if(
lgFileReadable(
"0040000_5.00_33_50_02_15.bin_0.1", dum, as ) &&
974 fprintf(
ioQQQ,
" Creating rauch_pg1159.ascii....\n" );
980 if(
lgFileReadable(
"0020000_4.00_H_00005-02000A.bin_0.1", dum, as ) &&
983 fprintf(
ioQQQ,
" Creating rauch_hydr.ascii....\n" );
989 if(
lgFileReadable(
"0050000_5.00_He_00005-02000A.bin_0.1", dum, as ) &&
992 fprintf(
ioQQQ,
" Creating rauch_helium.ascii....\n" );
993 lgFail = lgFail ||
RauchInitializeSub(
"rauch_helium.ascii",
"_He_00005-02000A.bin_0.1",
998 if(
lgFileReadable(
"0050000_5.00_H+He_1.000_0.000_00005-02000A.bin_0.1", dum, as ) &&
1001 fprintf(
ioQQQ,
" Creating rauch_h+he_3d.ascii....\n" );
1002 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_1.000_0.000_00005-02000A.bin_0.1",
1004 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.900_0.100_00005-02000A.bin_0.1",
1006 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.800_0.200_00005-02000A.bin_0.1",
1008 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.700_0.300_00005-02000A.bin_0.1",
1010 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.600_0.400_00005-02000A.bin_0.1",
1012 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.500_0.500_00005-02000A.bin_0.1",
1014 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.400_0.600_00005-02000A.bin_0.1",
1016 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.300_0.700_00005-02000A.bin_0.1",
1018 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.200_0.800_00005-02000A.bin_0.1",
1020 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.100_0.900_00005-02000A.bin_0.1",
1022 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.000_1.000_00005-02000A.bin_0.1",
1038 Edges[0] = 0.99946789f;
1039 Edges[1] = 1.8071406f;
1040 Edges[2] = 3.9996377f;
1043 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_solar.ascii",
"rauch_h-ca_solar.mod",Edges,3L, pc );
1045 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_halo.ascii",
"rauch_h-ca_halo.mod", Edges, 3L, pc );
1047 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_3d.ascii",
"rauch_h-ca_3d.mod", Edges, 3L, pc );
1050 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_solar.ascii",
"rauch_h-ni_solar.mod",Edges,3L, pc );
1052 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_halo.ascii",
"rauch_h-ni_halo.mod", Edges, 3L, pc );
1054 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_3d.ascii",
"rauch_h-ni_3d.mod", Edges, 3L, pc );
1057 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_pg1159.ascii",
"rauch_pg1159.mod", Edges, 3L, pc );
1059 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_cowd.ascii",
"rauch_cowd.mod", Edges, 3L, pc );
1062 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_hydr.ascii",
"rauch_hydr.mod", Edges, 3L, pc );
1065 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_helium.ascii",
"rauch_helium.mod", Edges, 3L, pc );
1068 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h+he_3d.ascii",
"rauch_h+he_3d.mod", Edges, 3L, pc );
1086 grid.name =
"rauch_h-ca_3d.mod";
1088 grid.name = ( lgHalo ?
"rauch_h-ca_halo.mod" :
"rauch_h-ca_solar.mod" );
1092 grid.ident =
" H-Ca Rauch";
1094 grid.command =
"COMPILE STARS";
1120 grid.name =
"rauch_h-ni_3d.mod";
1122 grid.name = ( lgHalo ?
"rauch_h-ni_halo.mod" :
"rauch_h-ni_solar.mod" );
1126 grid.ident =
" H-Ni Rauch";
1128 grid.command =
"COMPILE STARS";
1152 grid.name =
"rauch_pg1159.mod";
1156 grid.ident =
"PG1159 Rauch";
1158 grid.command =
"COMPILE STARS";
1182 grid.name =
"rauch_cowd.mod";
1186 grid.ident =
"C/O WD Rauch";
1188 grid.command =
"COMPILE STARS";
1212 grid.name =
"rauch_hydr.mod";
1216 grid.ident =
" Hydr Rauch";
1218 grid.command =
"COMPILE STARS";
1242 grid.name =
"rauch_helium.mod";
1246 grid.ident =
"Helium Rauch";
1248 grid.command =
"COMPILE STARS";
1272 grid.name =
"rauch_h+he_3d.mod";
1276 grid.ident =
" H+He Rauch";
1278 grid.command =
"COMPILE STARS";
1292 const char chOutName[],
1297 bool lgHeader =
true;
1298 long int i, j, nmods, ngp;
1302 double *wavl, *fluxes[
MNTS], Age[
MNTS], lwavl;
1309 for( i=0; i <
MNTS; i++ )
1313 wavl = (
double *)
MALLOC(
sizeof(
double)*nsb_sz);
1325 double cage, cwavl, cfl, cfl1, cfl2, cfl3;
1329 if( sscanf( chLine,
" %le %le %le %le %le", &cage, &cwavl, &cfl1, &cfl2, &cfl3 ) != 5 )
1331 fprintf(
ioQQQ,
"syntax error in data of Starburst grid.\n" );
1351 fprintf(
ioQQQ,
"too many time steps in Starburst grid.\n" );
1352 fprintf(
ioQQQ,
"please increase MNTS and recompile.\n" );
1359 fluxes[nmods] = (
double *)
MALLOC(
sizeof(
double)*nsb_sz);
1363 if( ngp >= (
long)nsb_sz )
1369 fluxes[0] = (
double *)
REALLOC(fluxes[0],
sizeof(
double)*nsb_sz);
1370 wavl = (
double *)
REALLOC(wavl,
sizeof(
double)*nsb_sz);
1373 if( !
fp_equal(Age[nmods],cage,10) )
1375 fprintf(
ioQQQ,
"age error in Starburst grid.\n" );
1383 if( !
fp_equal(wavl[ngp],cwavl,10) )
1385 fprintf(
ioQQQ,
"wavelength error in Starburst grid.\n" );
1392 fluxes[nmods][ngp] = pow( 10., cfl - 44.077911 );
1398 if( lgHeader && strncmp( &chLine[1],
"TIME [YR]", 9 ) == 0 )
1405 fprintf(
ioQQQ,
"syntax error in header of Starburst grid.\n" );
1417 fprintf( ioOut,
" %d\n", 1 );
1418 fprintf( ioOut,
" %d\n", 1 );
1419 fprintf( ioOut,
" Age\n" );
1420 fprintf( ioOut,
" %ld\n", nmods );
1421 fprintf( ioOut,
" %ld\n", ngp );
1423 fprintf( ioOut,
" lambda\n" );
1425 fprintf( ioOut,
" %.8e\n", 1. );
1427 fprintf( ioOut,
" F_lambda\n" );
1429 fprintf( ioOut,
" %.8e\n", 1. );
1431 for( i=0; i < nmods; i++ )
1433 fprintf( ioOut,
" %.3e", Age[i] );
1434 if( ((i+1)%4) == 0 )
1435 fprintf( ioOut,
"\n" );
1438 fprintf( ioOut,
"\n" );
1440 fprintf(
ioQQQ,
" Writing: " );
1443 for( j=0; j < ngp; j++ )
1445 fprintf( ioOut,
" %.4e", wavl[j] );
1447 if( ((j+1)%5) == 0 )
1448 fprintf( ioOut,
"\n" );
1452 fprintf( ioOut,
"\n" );
1455 fprintf(
ioQQQ,
"." );
1458 for( i=0; i < nmods; i++ )
1460 for( j=0; j < ngp; j++ )
1462 fprintf( ioOut,
" %.4e", fluxes[i][j] );
1464 if( ((j+1)%5) == 0 )
1465 fprintf( ioOut,
"\n" );
1469 fprintf( ioOut,
"\n" );
1472 fprintf(
ioQQQ,
"." );
1476 fprintf(
ioQQQ,
" done.\n" );
1481 for( i=0; i <
MNTS; i++ )
1487 for( i=0; i <
MNTS; i++ )
1498 bool lgFail =
false;
1502 fprintf(
ioQQQ,
" StarburstCompile on the job.\n" );
1510 lgFail = lgFail ||
lgCompileAtmosphere(
"starburst99.ascii",
"starburst99.mod", Edges, 0L, pc );
1513 lgFail = lgFail ||
lgCompileAtmosphere(
"starburst99_2d.ascii",
"starburst99_2d.mod", Edges, 0L, pc );
1523 bool lgFail =
false;
1527 fprintf(
ioQQQ,
" TlustyCompile on the job.\n" );
1532 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_p03.ascii",
"obstar_merged_p03.mod", Edges, 0L, pc );
1534 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_p00.ascii",
"obstar_merged_p00.mod", Edges, 0L, pc );
1536 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m03.ascii",
"obstar_merged_m03.mod", Edges, 0L, pc );
1538 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m07.ascii",
"obstar_merged_m07.mod", Edges, 0L, pc );
1540 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m10.ascii",
"obstar_merged_m10.mod", Edges, 0L, pc );
1542 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m99.ascii",
"obstar_merged_m99.mod", Edges, 0L, pc );
1545 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_3d.ascii",
"obstar_merged_3d.mod", Edges, 0L, pc );
1548 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_p03.ascii",
"bstar2006_p03.mod", Edges, 0L, pc );
1550 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_p00.ascii",
"bstar2006_p00.mod", Edges, 0L, pc );
1552 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m03.ascii",
"bstar2006_m03.mod", Edges, 0L, pc );
1554 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m07.ascii",
"bstar2006_m07.mod", Edges, 0L, pc );
1556 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m10.ascii",
"bstar2006_m10.mod", Edges, 0L, pc );
1558 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m99.ascii",
"bstar2006_m99.mod", Edges, 0L, pc );
1561 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_3d.ascii",
"bstar2006_3d.mod", Edges, 0L, pc );
1564 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_p03.ascii",
"ostar2002_p03.mod", Edges, 0L, pc );
1566 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_p00.ascii",
"ostar2002_p00.mod", Edges, 0L, pc );
1568 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m03.ascii",
"ostar2002_m03.mod", Edges, 0L, pc );
1570 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m07.ascii",
"ostar2002_m07.mod", Edges, 0L, pc );
1572 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m10.ascii",
"ostar2002_m10.mod", Edges, 0L, pc );
1574 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m15.ascii",
"ostar2002_m15.mod", Edges, 0L, pc );
1576 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m17.ascii",
"ostar2002_m17.mod", Edges, 0L, pc );
1578 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m20.ascii",
"ostar2002_m20.mod", Edges, 0L, pc );
1580 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m30.ascii",
"ostar2002_m30.mod", Edges, 0L, pc );
1582 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m99.ascii",
"ostar2002_m99.mod", Edges, 0L, pc );
1585 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_3d.ascii",
"ostar2002_3d.mod", Edges, 0L, pc );
1594 const char *chMetalicity,
1605 grid.name =
"obstar_merged_";
1607 grid.name =
"bstar2006_";
1609 grid.name =
"ostar2002_";
1615 grid.name += chMetalicity;
1616 grid.name +=
".mod";
1622 strcpy( chIdent,
"3-dim" );
1626 strcpy( chIdent,
"Z " );
1627 strcat( chIdent, chMetalicity );
1630 strcat( chIdent,
" OBstar" );
1632 strcat( chIdent,
" Bstr06" );
1634 strcat( chIdent,
" Ostr02" );
1637 grid.ident = chIdent;
1639 grid.command =
"COMPILE STARS";
1658 bool lgFail =
false;
1662 fprintf(
ioQQQ,
" WernerCompile on the job.\n" );
1676 Edges[0] = 0.99946789f;
1677 Edges[1] = 1.8071406f;
1678 Edges[2] = 3.9996377f;
1725 grid.name =
"kwerner.mod";
1729 grid.ident =
"Klaus Werner";
1731 grid.command =
"COMPILE STARS";
1769 bool lgFail =
false;
1773 fprintf(
ioQQQ,
" WMBASICCompile on the job.\n" );
1776 Edges[0] = 0.99946789f;
1777 Edges[1] = 1.8071406f;
1778 Edges[2] = 3.9996377f;
1799 grid.name =
"wmbasic.mod";
1803 grid.ident =
" WMBASIC";
1805 grid.command =
"COMPILE STARS";
1820 const char chFNameOut[],
1831 long int i, j, nskip, nModels, nWL;
1834 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL;
1840 vector<realnum> SaveAnu(
rfield.nupper);
1862 fprintf(
ioQQQ,
" CompileAtmosphereCoStar got %s.\n", chFNameIn );
1867 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading nskip.\n" );
1870 sscanf( chLine,
"%li", &nskip );
1873 for( i=0; i < nskip; ++i )
1877 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails skipping header.\n" );
1885 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading nModels, nWL.\n" );
1888 sscanf( chLine,
"%li%li", &nModels, &nWL );
1890 if( nModels <= 0 || nWL <= 0 )
1892 fprintf(
ioQQQ,
" CompileAtmosphereCoStar scanned off impossible values for nModels=%li or nWL=%li\n",
1901 for( i=0; i < nModels; ++i )
1905 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading model parameters.\n" );
1909 telg[i].
chGrid = chLine[0];
1911 sscanf( chLine+1,
"%i", &telg[i].modid );
1913 sscanf( chLine+23,
"%lg", &telg[i].par[0] );
1915 sscanf( chLine+31,
"%lg", &telg[i].par[1] );
1917 sscanf( chLine+7,
"%lg", &telg[i].par[2] );
1919 sscanf( chLine+15,
"%lg", &telg[i].par[3] );
1922 ASSERT( telg[i].par[2] > 10. );
1923 ASSERT( telg[i].par[3] > 10. );
1926 telg[i].
par[2] = log10(telg[i].par[2]);
1941 val[1] = (int32)
MDIM;
1942 val[2] = (int32)
MNAM;
1945 val[5] = (int32)nModels;
1946 val[6] = (int32)
rfield.nupper;
1947 uval[0] =
sizeof(val) +
sizeof(uval) +
sizeof(dval) +
sizeof(md5sum) +
1948 sizeof(names) + nModels*
sizeof(
mpp);
1950 dval[0] = double(
rfield.emm);
1951 dval[1] = double(
rfield.egamry);
1952 dval[2] = double(
continuum.ResolutionScaleFactor);
1954 strncpy( md5sum,
continuum.mesh_md5sum.c_str(),
sizeof(md5sum) );
1956 strncpy( names[0],
"Teff\0\0",
MNAM+1 );
1957 strncpy( names[1],
"log(g)",
MNAM+1 );
1958 strncpy( names[2],
"log(M)",
MNAM+1 );
1959 strncpy( names[3],
"Age\0\0\0",
MNAM+1 );
1961 for (
long i=0; i<
rfield.nupper; ++i)
1963 if( fwrite( val,
sizeof(val), 1, ioOUT ) != 1 ||
1964 fwrite( uval,
sizeof(uval), 1, ioOUT ) != 1 ||
1966 fwrite( dval,
sizeof(dval), 1, ioOUT ) != 1 ||
1968 fwrite( md5sum,
sizeof(md5sum), 1, ioOUT ) != 1 ||
1969 fwrite( names,
sizeof(names), 1, ioOUT ) != 1 ||
1971 fwrite( telg,
sizeof(
mpp), (
size_t)nModels, ioOUT ) != (
size_t)nModels ||
1973 fwrite(
get_ptr(SaveAnu), (
size_t)uval[1], 1, ioOUT ) != 1 )
1975 fprintf(
ioQQQ,
" CompileAtmosphereCoStar failed writing header of output file.\n" );
1984 fprintf(
ioQQQ,
" Compiling: " );
1987 for( i=0; i < nModels; ++i )
1992 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading the skip to next spectrum.\n" );
1995 sscanf( chLine,
"%li", &nskip );
1997 for( j=0; j < nskip; ++j )
2001 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails doing the skip.\n" );
2009 for( j=nWL-1; j >= 0; --j )
2013 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading the spectral data.\n" );
2016 double help1, help2;
2017 sscanf( chLine,
"%lg %lg", &help1, &help2 );
2021 StarFlux[j] = (
realnum)(
PI*pow(10.,help2));
2027 ASSERT( StarEner[j] < StarEner[j+1] );
2033 RebinAtmosphere(nWL-1, StarEner+1, StarFlux+1, CloudyFlux, nedges, Edges );
2036 if( fwrite( CloudyFlux, (
size_t)uval[1], 1, ioOUT ) != 1 )
2038 fprintf(
ioQQQ,
" CompileAtmosphereCoStar failed writing star flux.\n" );
2042 fprintf(
ioQQQ,
"." );
2046 fprintf(
ioQQQ,
" done.\n" );
2056 fprintf(
ioQQQ,
"\n CompileAtmosphereCoStar completed ok.\n\n" );
2078 long i, j, k, nmodid, off, ptr;
2079 long *indloTr, *indhiTr, useTr[2];
2080 long indlo[2], indhi[2], index[2];
2082 double lval[2], aval[4];
2086 switch(
grid->imode )
2095 lval[0] = log10(val[0]);
2101 lval[0] = log10(val[1]);
2106 fprintf(
ioQQQ,
" InterpolateGridCoStar called with insane value for imode: %d.\n",
grid->imode );
2110 nmodid = (long)(lval[1]+0.5);
2124 indloTr = (
long *)
MALLOC(
sizeof(
long)*
grid->nTracks );
2125 indhiTr = (
long *)
MALLOC(
sizeof(
long)*
grid->nTracks );
2128 for( j=0; j <
grid->nTracks; j++ )
2132 if(
grid->trackLen[j] >= nmodid ) {
2133 index[0] = nmodid - 1;
2144 ValTr[j] = -FLT_MAX;
2154 for( j=0; j <
grid->nTracks; j++ )
2156 if( indloTr[j] >= 0 )
2157 printf(
"track %c: models %c%d, %c%d, val %g\n",
2158 (
char)(
'A'+j),
grid->telg[indloTr[j]].chGrid,
grid->telg[indloTr[j]].modid,
2159 grid->telg[indhiTr[j]].chGrid,
grid->telg[indhiTr[j]].modid, ValTr[j]);
2172 fprintf(
ioQQQ,
" The parameters for the requested CoStar model are out of range.\n" );
2176 ASSERT( useTr[0] >= 0 && useTr[0] <
grid->nTracks );
2177 ASSERT( useTr[1] >= 0 && useTr[1] <
grid->nTracks );
2178 ASSERT( indloTr[useTr[0]] >= 0 && indloTr[useTr[0]] < (
int)
grid->nmods );
2179 ASSERT( indhiTr[useTr[0]] >= 0 && indhiTr[useTr[0]] < (
int)
grid->nmods );
2180 ASSERT( indloTr[useTr[1]] >= 0 && indloTr[useTr[1]] < (
int)
grid->nmods );
2181 ASSERT( indhiTr[useTr[1]] >= 0 && indhiTr[useTr[1]] < (
int)
grid->nmods );
2184 printf(
"interpolate between tracks %c and %c\n", (
char)(
'A'+useTr[0]), (
char)(
'A'+useTr[1]) );
2187 indlo[0] = indloTr[useTr[0]];
2188 indhi[0] = indhiTr[useTr[0]];
2189 indlo[1] = indloTr[useTr[1]];
2190 indhi[1] = indhiTr[useTr[1]];
2194 for( i=0; i <
rfield.nupper; i++ )
2203 FILE *ioBUG = fopen(
"interpolated.txt",
"w" );
2204 for( k=0; k <
rfield.nupper; ++k )
2214 SetLimits(
grid, lval[0], NULL, NULL, useTr, ValTr, val0_lo, val0_hi );
2219 fprintf(
ioQQQ,
" * c<< FINAL: T_eff = %7.1f, ", aval[0] );
2220 fprintf(
ioQQQ,
"log(g) = %4.2f, M(ZAMS) = %5.1f, age = ", aval[1], pow(10.,aval[2]) );
2222 fprintf(
ioQQQ,
" >>> *\n" );
2240 long index[2], j, mod1, mod2;
2244 indloTr[track] = -2;
2245 indhiTr[track] = -2;
2246 ValTr[track] = -FLT_MAX;
2250 for( j=0; j <
grid->trackLen[track]; j++ )
2256 if( fabs(par2-
grid->telg[mod1].par[off+1]) <= 10.*FLT_EPSILON*fabs(
grid->telg[mod1].par[off+1]) )
2258 indloTr[track] = mod1;
2259 indhiTr[track] = mod1;
2260 ValTr[track] = (
realnum)
grid->telg[mod1].par[off];
2265 for( j=0; j <
grid->trackLen[track]-1; j++ )
2273 if( (par2 -
grid->telg[mod1].par[off+1])*(par2 -
grid->telg[mod2].par[off+1]) < 0. )
2277 indloTr[track] = mod1;
2278 indhiTr[track] = mod2;
2279 frac = (par2 -
grid->telg[mod2].par[off+1])/
2280 (
grid->telg[mod1].par[off+1] -
grid->telg[mod2].par[off+1]);
2281 ValTr[track] = (
realnum)(frac*
grid->telg[mod1].par[off] +
2282 (1.-frac)*
grid->telg[mod2].par[off] );
2306 for( j=0; j <
grid->nTracks; j++ )
2309 if( ValTr[j] != -FLT_MAX && fabs(par1-(
double)ValTr[j]) <= 10.*FLT_EPSILON*fabs(ValTr[j]) )
2322 for( j=0; j <
grid->nTracks-1; j++ )
2324 if( ValTr[j] != -FLT_MAX )
2330 for( i = j+1; i <
grid->nTracks; i++ )
2332 if( ValTr[i] != -FLT_MAX )
2340 if( j2 > 0 && ((
realnum)par1-ValTr[j])*((
realnum)par1-ValTr[j2]) < 0.f )
2350 continuum.lgCoStarInterpolationCaution = ( useTr[1]-useTr[0] > 1 );
2357 long index[2], maxlen, n;
2362 for( n=0; n <
grid->nTracks; n++ )
2363 maxlen =
MAX2( maxlen,
grid->trackLen[n] );
2365 fprintf(
ioQQQ,
"\n" );
2366 fprintf(
ioQQQ,
" Track\\Index |" );
2367 for( n = 0; n < maxlen; n++ )
2368 fprintf(
ioQQQ,
" %5ld ", n+1 );
2369 fprintf(
ioQQQ,
"\n" );
2370 fprintf(
ioQQQ,
"--------------|" );
2371 for( n = 0; n < maxlen; n++ )
2372 fprintf(
ioQQQ,
"----------------" );
2373 fprintf(
ioQQQ,
"\n" );
2375 for( index[1]=0; index[1] <
grid->nTracks; ++index[1] )
2378 double Teff, alogg, Mass;
2380 fprintf(
ioQQQ,
" %c", (
char)(
'A'+index[1]) );
2383 Mass = pow(10.,
grid->telg[ptr].par[2]);
2384 fprintf(
ioQQQ,
" (%3.0f Msol) |", Mass );
2386 for( index[0]=0; index[0] <
grid->trackLen[index[1]]; ++index[0] )
2389 Teff =
grid->telg[ptr].par[0];
2390 alogg =
grid->telg[ptr].par[1];
2391 fprintf(
ioQQQ,
" (%6.1f,%4.2f)", Teff, alogg );
2393 fprintf(
ioQQQ,
"\n" );
2400 const char chSuff[],
2401 const vector<mpp>& telg,
2405 const double par2[],
2415 double *wavl, *fluxes;
2438 fprintf( ioOut,
" %d\n", ( ngrids == 1 ? 2 : 3 ) );
2439 fprintf( ioOut,
" %d\n", ( ngrids == 1 ? 2 : 3 ) );
2440 fprintf( ioOut,
" Teff\n" );
2441 fprintf( ioOut,
" log(g)\n" );
2443 fprintf( ioOut,
" log(Z)\n" );
2444 else if( ngrids == 11 )
2445 fprintf( ioOut,
" f(He)\n" );
2447 fprintf( ioOut,
" %ld\n", nmods*ngrids );
2448 fprintf( ioOut,
" %d\n",
NRAUCH );
2450 fprintf( ioOut,
" lambda\n" );
2452 fprintf( ioOut,
" %.8e\n", 1. );
2454 fprintf( ioOut,
" F_lambda\n" );
2456 fprintf( ioOut,
" %.8e\n",
PI*1.e-8 );
2458 for( j=0; j < ngrids; j++ )
2461 for( i=0; i < nmods; i++ )
2464 fprintf( ioOut,
" %.0f %.1f", telg[i].par[0], telg[i].par[1] );
2466 fprintf( ioOut,
" %.0f %.1f %.1f", telg[i].par[0], telg[i].par[1], par2[j] );
2467 if( ((i+1)%4) == 0 )
2468 fprintf( ioOut,
"\n" );
2471 fprintf( ioOut,
"\n" );
2474 fprintf(
ioQQQ,
" Writing: " );
2477 for( i=0; i < nmods; i++ )
2481 sprintf( chLine,
"%7.7ld_%2ld", (
long)(telg[i].par[0]+0.5), (
long)(10.*telg[i].par[1]+0.5) );
2482 else if( format == 2 )
2483 sprintf( chLine,
"%7.7ld_%.2f", (
long)(telg[i].par[0]+0.5), telg[i].par[1] );
2486 fprintf(
ioQQQ,
" insanity in RauchInitializeSub\n" );
2490 string chFileName( chLine );
2491 chFileName += chSuff;
2506 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2512 while( chLine[0] ==
'*' )
2516 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2523 for( j=0; j <
NRAUCH; j++ )
2532 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2539 if( sscanf( chLine,
"%lf %le", &wl, &ttemp ) != 2 )
2541 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2553 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2565 if( i == 0 && n == 1 )
2568 for( j=0; j <
NRAUCH; j++ )
2570 fprintf( ioOut,
" %.4e", wavl[j] );
2572 if( ((j+1)%5) == 0 )
2573 fprintf( ioOut,
"\n" );
2577 fprintf( ioOut,
"\n" );
2580 for( j=0; j <
NRAUCH; j++ )
2582 fprintf( ioOut,
" %.4e", fluxes[j] );
2584 if( ((j+1)%5) == 0 )
2585 fprintf( ioOut,
"\n" );
2589 fprintf( ioOut,
"\n" );
2592 fprintf(
ioQQQ,
"." );
2597 fprintf(
ioQQQ,
" done.\n" );
2621 const char fnam[] =
"rauch_models.dat";
2628 istringstream iss( line );
2632 fprintf(
ioQQQ,
" RauchReadMPP: the version of %s is not the current version.\n", fnam );
2633 fprintf(
ioQQQ,
" Please obtain the current version from the Cloudy web site.\n" );
2634 fprintf(
ioQQQ,
" I expected to find version %ld and got %ld instead.\n",
2640 unsigned long ndata;
2641 istringstream iss2( line );
2643 ASSERT( ndata == telg1.size() );
2646 getline( ioDATA, line );
2648 for(
unsigned long i=0; i < ndata; ++i )
2649 ioDATA >> telg1[i].par[0] >> telg1[i].par[1];
2650 getline( ioDATA, line );
2653 istringstream iss3( line );
2655 ASSERT( ndata == telg2.size() );
2656 getline( ioDATA, line );
2658 for(
unsigned long i=0; i < ndata; ++i )
2659 ioDATA >> telg2[i].par[0] >> telg2[i].par[1];
2660 getline( ioDATA, line );
2663 istringstream iss4( line );
2665 ASSERT( ndata == telg3.size() );
2666 getline( ioDATA, line );
2668 for(
unsigned long i=0; i < ndata; ++i )
2669 ioDATA >> telg3[i].par[0] >> telg3[i].par[1];
2670 getline( ioDATA, line );
2673 istringstream iss5( line );
2675 ASSERT( ndata == telg4.size() );
2676 getline( ioDATA, line );
2678 for(
unsigned long i=0; i < ndata; ++i )
2679 ioDATA >> telg4[i].par[0] >> telg4[i].par[1];
2680 getline( ioDATA, line );
2683 istringstream iss6( line );
2685 ASSERT( ndata == telg5.size() );
2686 getline( ioDATA, line );
2688 for(
unsigned long i=0; i < ndata; ++i )
2689 ioDATA >> telg5[i].par[0] >> telg5[i].par[1];
2690 getline( ioDATA, line );
2693 istringstream iss7( line );
2695 ASSERT( ndata == telg6.size() );
2696 getline( ioDATA, line );
2698 for(
unsigned long i=0; i < ndata; ++i )
2699 ioDATA >> telg6[i].par[0] >> telg6[i].par[1];
2700 getline( ioDATA, line );
2703 istringstream iss8( line );
2715 getline( ioDATA, line );
2717 while( line[0] ==
'#' );
2724 const char chFNameOut[],
2732 char chDataType[11];
2735 bool lgFreqX, lgFreqY, lgFlip;
2740 long int i, imod, version, nd, ndim, npar, nmods, ngrid;
2743 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL, *scratch = NULL;
2744 vector<realnum> SaveAnu(
rfield.nupper);
2746 double convert_wavl, convert_flux;
2760 fprintf(
ioQQQ,
" lgCompileAtmosphere got %s.\n", chFNameIn );
2763 if( fscanf( ioIN,
"%ld", &version ) != 1 )
2765 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading VERSION.\n" );
2771 fprintf(
ioQQQ,
" lgCompileAtmosphere: there is a version number mismatch in"
2772 " the ascii atmosphere file: %s.\n", chFNameIn );
2773 fprintf(
ioQQQ,
" lgCompileAtmosphere: Please recreate this file or download the"
2774 " latest version following the instructions on the Cloudy website.\n" );
2779 if( fscanf( ioIN,
"%ld", &ndim ) != 1 || ndim <= 0 || ndim >
MDIM )
2781 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid dimension of grid.\n" );
2786 if( fscanf( ioIN,
"%ld", &npar ) != 1 || npar <= 0 || npar < ndim || npar >
MDIM )
2788 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid no. of model parameters.\n" );
2793 memset( names,
'\0',
MDIM*(
MNAM+1) );
2795 for( nd=0; nd < npar; nd++ )
2797 if( fscanf( ioIN,
"%6s", names[nd] ) != 1 )
2799 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading parameter label.\n" );
2805 if( fscanf( ioIN,
"%ld", &nmods ) != 1 || nmods <= 0 )
2807 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid number of models.\n" );
2811 if( fscanf( ioIN,
"%ld", &ngrid ) != 1 || ngrid <= 1 )
2813 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid number of grid points.\n" );
2818 if( fscanf( ioIN,
"%10s", chDataType ) != 1 )
2820 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading wavl DataType string.\n" );
2824 if( strcmp( chDataType,
"lambda" ) == 0 )
2826 else if( strcmp( chDataType,
"nu" ) == 0 )
2829 fprintf(
ioQQQ,
" lgCompileAtmosphere found illegal wavl DataType: %s.\n", chDataType );
2833 if( fscanf( ioIN,
"%le", &convert_wavl ) != 1 || convert_wavl <= 0. )
2835 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid wavl conversion factor.\n" );
2840 if( fscanf( ioIN,
"%10s", chDataType ) != 1 )
2842 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading flux DataType string.\n" );
2846 if( strcmp( chDataType,
"F_lambda" ) == 0 || strcmp( chDataType,
"H_lambda" ) == 0 )
2848 else if( strcmp( chDataType,
"F_nu" ) == 0 || strcmp( chDataType,
"H_nu" ) == 0 )
2851 fprintf(
ioQQQ,
" lgCompileAtmosphere found illegal flux DataType: %s.\n", chDataType );
2855 if( fscanf( ioIN,
"%le", &convert_flux ) != 1 || convert_flux <= 0. )
2857 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid flux conversion factor.\n" );
2863 for( i=0; i < nmods; i++ )
2865 for( nd=0; nd < npar; nd++ )
2867 if( fscanf( ioIN,
"%le", &telg[i].par[nd] ) != 1 )
2869 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid model parameter.\n" );
2885 val[1] = (int32)
MDIM;
2886 val[2] = (int32)
MNAM;
2887 val[3] = (int32)ndim;
2888 val[4] = (int32)npar;
2889 val[5] = (int32)nmods;
2890 val[6] = (int32)
rfield.nupper;
2891 uval[0] =
sizeof(val) +
sizeof(uval) +
sizeof(dval) +
sizeof(md5sum) +
2892 sizeof(names) + nmods*
sizeof(
mpp);
2894 dval[0] = double(
rfield.emm);
2895 dval[1] = double(
rfield.egamry);
2896 dval[2] = double(
continuum.ResolutionScaleFactor);
2898 strncpy( md5sum,
continuum.mesh_md5sum.c_str(),
sizeof(md5sum) );
2900 for (
long i=0; i<
rfield.nupper; ++i)
2903 if( fwrite( val,
sizeof(val), 1, ioOUT ) != 1 ||
2904 fwrite( uval,
sizeof(uval), 1, ioOUT ) != 1 ||
2906 fwrite( dval,
sizeof(dval), 1, ioOUT ) != 1 ||
2908 fwrite( md5sum,
sizeof(md5sum), 1, ioOUT ) != 1 ||
2909 fwrite( names,
sizeof(names), 1, ioOUT ) != 1 ||
2911 fwrite( telg,
sizeof(
mpp), (
size_t)nmods, ioOUT ) != (
size_t)nmods ||
2913 fwrite(
get_ptr(SaveAnu), (
size_t)uval[1], 1, ioOUT ) != 1 )
2915 fprintf(
ioQQQ,
" lgCompileAtmosphere failed writing header of output file.\n" );
2926 for( i=0; i < ngrid; i++ )
2929 if( fscanf( ioIN,
"%lg", &help ) != 1 )
2931 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading wavelength.\n" );
2936 scratch[i] = (
realnum)(help*convert_wavl);
2938 if( scratch[i] <= 0.f )
2940 fprintf(
ioQQQ,
" PROBLEM: a non-positive %s was found, value: %e\n",
2941 lgFreqX ?
"frequency" :
"wavelength", help );
2946 lgFlip = ( !lgFreqX && scratch[0] < scratch[1] ) || ( lgFreqX && scratch[0] > scratch[1] );
2949 for( i=0; i < ngrid; i++ )
2958 StarEner[ngrid-i-1] = scratch[i];
2960 StarEner[i] = scratch[i];
2963 ASSERT( StarEner[0] > 0.f );
2965 for( i=1; i < ngrid; i++ )
2967 if( StarEner[i] <= StarEner[i-1] )
2969 fprintf(
ioQQQ,
" PROBLEM: the %s grid is not strictly monotonically increasing/decreasing\n",
2970 lgFreqX ?
"frequency" :
"wavelength" );
2975 fprintf(
ioQQQ,
" Compiling: " );
2977 for( imod=0; imod < nmods; imod++ )
2982 for( i=0; i < ngrid; i++ )
2985 if( fscanf( ioIN,
"%lg", &help ) != 1 )
2987 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading star flux.\n" );
2992 scratch[i] = (
realnum)(help*convert_flux);
2995 if( scratch[i] < 0.f )
2997 fprintf(
ioQQQ,
"\n PROBLEM: a negative flux was found, model number %ld, value: %e\n",
3003 for( i=0; i < ngrid; i++ )
3006 StarFlux[ngrid-i-1] = scratch[i];
3008 StarFlux[i] = scratch[i];
3011 for( i=0; i < ngrid; i++ )
3015 StarFlux[i] *= CONVERT_FNU/
POW2(StarEner[i]);
3016 ASSERT( StarFlux[i] >= 0.f );
3020 RebinAtmosphere( ngrid, StarEner, StarFlux, CloudyFlux, nedges, Edges );
3023 if( fwrite( CloudyFlux, (
size_t)uval[1], 1, ioOUT ) != 1 )
3025 fprintf(
ioQQQ,
" lgCompileAtmosphere failed writing star flux.\n" );
3029 fprintf(
ioQQQ,
"." );
3033 fprintf(
ioQQQ,
" done.\n" );
3045 fprintf(
ioQQQ,
" lgCompileAtmosphere completed ok.\n\n" );
3063 int32 version, mdim, mnam;
3064 double mesh_elo, mesh_ehi;
3077 fprintf(
ioQQQ,
" Error: stellar atmosphere file not found.\n" );
3078 fprintf(
ioQQQ ,
"\n\n If the path is set then it is possible that the stellar atmosphere data files do not exist.\n");
3079 fprintf(
ioQQQ ,
" Have the stellar data files been downloaded and compiled with the COMPILE STARS command?\n");
3080 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");
3085 if( fread( &version,
sizeof(version), 1,
grid->ioIN ) != 1 ||
3086 fread( &mdim,
sizeof(mdim), 1,
grid->ioIN ) != 1 ||
3087 fread( &mnam,
sizeof(mnam), 1,
grid->ioIN ) != 1 ||
3088 fread( &
grid->ndim,
sizeof(
grid->ndim), 1,
grid->ioIN ) != 1 ||
3089 fread( &
grid->npar,
sizeof(
grid->npar), 1,
grid->ioIN ) != 1 ||
3090 fread( &
grid->nmods,
sizeof(
grid->nmods), 1,
grid->ioIN ) != 1 ||
3091 fread( &
grid->ngrid,
sizeof(
grid->ngrid), 1,
grid->ioIN ) != 1 ||
3092 fread( &
grid->nOffset,
sizeof(
grid->nOffset), 1,
grid->ioIN ) != 1 ||
3093 fread( &
grid->nBlocksize,
sizeof(
grid->nBlocksize), 1,
grid->ioIN ) != 1 ||
3094 fread( &mesh_elo,
sizeof(mesh_elo), 1,
grid->ioIN ) != 1 ||
3095 fread( &mesh_ehi,
sizeof(mesh_ehi), 1,
grid->ioIN ) != 1 ||
3097 fread( md5sum,
sizeof(md5sum), 1,
grid->ioIN ) != 1 )
3099 fprintf(
ioQQQ,
" InitGrid failed reading header.\n" );
3106 fprintf(
ioQQQ,
" InitGrid: there is a version mismatch between"
3107 " the compiled atmospheres file I expected and the one I found.\n" );
3108 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3109 " atmospheres file with the command: %s.\n",
grid->command );
3115 fprintf(
ioQQQ,
" InitGrid: the compiled atmospheres file is produced"
3116 " with an incompatible version of Cloudy.\n" );
3117 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3118 " atmospheres file with the command: %s.\n",
grid->command );
3124 strncmp(
continuum.mesh_md5sum.c_str(), md5sum,
NMD5 ) != 0 )
3126 fprintf(
ioQQQ,
" InitGrid: the compiled atmospheres file is produced"
3127 " with an incompatible frequency grid.\n" );
3128 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3129 " atmospheres file with the command: %s.\n",
grid->command );
3142 if( fread( &
grid->names,
sizeof(
grid->names), 1,
grid->ioIN ) != 1 )
3144 fprintf(
ioQQQ,
" InitGrid failed reading names array.\n" );
3148 grid->lgIsTeffLoggGrid = (
grid->ndim >= 2 &&
3149 strcmp(
grid->names[0],
"Teff" ) == 0 &&
3150 strcmp(
grid->names[1],
"log(g)" ) == 0 );
3154 for( nd=0; nd <
grid->ndim; nd++ )
3162 fprintf(
ioQQQ,
" InitGrid failed reading model parameter block.\n" );
3170 int res = fseek(
grid->ioIN, 0, SEEK_END );
3173 long End = ftell(
grid->ioIN );
3174 long Expected =
grid->nOffset + (
grid->nmods+1)*
grid->nBlocksize;
3175 if( End != Expected )
3177 fprintf(
ioQQQ,
" InitGrid: Problem performing sanity check for size of binary file.\n" );
3178 fprintf(
ioQQQ,
" InitGrid: I expected to find %ld bytes, but actually found %ld bytes.\n",
3180 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3181 " atmospheres file with the command: %s.\n",
grid->command );
3192 grid->trackLen = NULL;
3201 int32 version, mdim, mnam;
3202 double mesh_elo, mesh_ehi, mesh_res_factor;
3222 grid.name = binName;
3227 if( fread( &version,
sizeof(version), 1,
grid.ioIN ) != 1 ||
3228 fread( &mdim,
sizeof(mdim), 1,
grid.ioIN ) != 1 ||
3229 fread( &mnam,
sizeof(mnam), 1,
grid.ioIN ) != 1 ||
3230 fread( &
grid.ndim,
sizeof(
grid.ndim), 1,
grid.ioIN ) != 1 ||
3231 fread( &
grid.npar,
sizeof(
grid.npar), 1,
grid.ioIN ) != 1 ||
3232 fread( &
grid.nmods,
sizeof(
grid.nmods), 1,
grid.ioIN ) != 1 ||
3233 fread( &
grid.ngrid,
sizeof(
grid.ngrid), 1,
grid.ioIN ) != 1 ||
3234 fread( &
grid.nOffset,
sizeof(
grid.nOffset), 1,
grid.ioIN ) != 1 ||
3235 fread( &
grid.nBlocksize,
sizeof(
grid.nBlocksize), 1,
grid.ioIN ) != 1 ||
3236 fread( &mesh_elo,
sizeof(mesh_elo), 1,
grid.ioIN ) != 1 ||
3237 fread( &mesh_ehi,
sizeof(mesh_ehi), 1,
grid.ioIN ) != 1 ||
3238 fread( &mesh_res_factor,
sizeof(mesh_res_factor), 1,
grid.ioIN ) != 1 ||
3239 fread( md5sum,
sizeof(md5sum), 1,
grid.ioIN ) != 1 )
3241 fclose(
grid.ioIN );
3250 strncmp(
continuum.mesh_md5sum.c_str(), md5sum,
NMD5 ) != 0 )
3252 fclose(
grid.ioIN );
3260 int res = fseek(
grid.ioIN, 0, SEEK_END );
3263 long End = ftell(
grid.ioIN );
3264 long Expected =
grid.nOffset + (
grid.nmods+1)*
grid.nBlocksize;
3265 if( End != Expected )
3267 fclose(
grid.ioIN );
3273 fclose(
grid.ioIN );
3287 if( (ioIN =
open_data( ascName,
"r", scheme )) == NULL )
3291 if( fscanf( ioIN,
"%ld", &version ) != 1 || version !=
VERSION_ASCII )
3318 memset(
grid->jval, 0xff, (
size_t)(
grid->nval[0]*
grid->nval[1]*
sizeof(
long)) );
3320 grid->trackLen = (
long *)
CALLOC( (
size_t)
grid->nmods,
sizeof(
long) );
3326 track = (char)(
'A'+index[1]);
3330 for( i=0; i <
grid->nmods; i++ )
3332 if(
grid->telg[i].chGrid == track &&
grid->telg[i].modid == index[0]+1 )
3346 grid->trackLen[index[1]] = index[0];
3350 grid->nTracks = index[1];
3362 *ndim = (long)
grid->ndim;
3363 if( *ndim == 2 && *nval == 1 &&
grid->lgIsTeffLoggGrid )
3366 val[*nval] =
grid->val[1][
grid->nval[1]-1];
3369 if( *ndim != (
long)
grid->ndim )
3371 fprintf(
ioQQQ,
" A %ld-dim grid was requested, but a %ld-dim grid was found.\n",
3372 *ndim, (
long)
grid->ndim );
3377 fprintf(
ioQQQ,
" A %ld-dim grid was requested, but only %ld parameters were entered.\n",
3400 indlo = (
long *)
MALLOC((
size_t)(
grid->ndim*
sizeof(
long)) );
3401 indhi = (
long *)
MALLOC((
size_t)(
grid->ndim*
sizeof(
long)) );
3402 index = (
long *)
MALLOC((
size_t)(
grid->ndim*
sizeof(
long)) );
3403 aval = (
double *)
MALLOC((
size_t)(
grid->npar*
sizeof(
double)) );
3417 for( nd=0; nd <
grid->ndim; nd++ )
3419 FindIndex(
grid->val[nd],
grid->nval[nd], val[nd], &indlo[nd], &indhi[nd], &lgInvalid );
3423 " Requested parameter %s = %.2f is not within the range %.2f to %.2f\n",
3424 grid->names[nd], val[nd],
grid->val[nd][0],
grid->val[nd][
grid->nval[nd]-1] );
3434 if(
grid->npar == 1 )
3436 " * c<< FINAL: %6s = %13.2f"
3438 grid->names[0], aval[0] );
3439 else if(
grid->npar == 2 )
3441 " * c<< FINAL: %6s = %10.2f"
3442 " %6s = %8.5f >>> *\n",
3443 grid->names[0], aval[0],
grid->names[1], aval[1] );
3444 else if(
grid->npar == 3 )
3446 " * c<< FINAL: %6s = %7.0f"
3447 " %6s = %5.2f %6s = %5.2f >>> *\n",
3448 grid->names[0], aval[0],
grid->names[1], aval[1],
3449 grid->names[2], aval[2] );
3450 else if(
grid->npar >= 4 )
3453 " * c<< FINAL: %4s = %7.0f"
3454 " %6s = %4.2f %6s = %5.2f %6s = ",
3455 grid->names[0], aval[0],
grid->names[1], aval[1],
3456 grid->names[2], aval[2],
grid->names[3] );
3458 fprintf(
ioQQQ,
" >>> *\n" );
3462 for( i=0; i <
rfield.nupper; i++ )
3471 FILE *ioBUG = fopen(
"interpolated.txt",
"w" );
3472 for( k=0; k <
rfield.nupper; ++k )
3477 if( strcmp(
grid->names[0],
"Teff" ) == 0 )
3484 SetLimits(
grid, val[0], indlo, indhi, NULL, NULL, Tlow, Thigh );
3501 fclose(
grid->ioIN );
3503 for( i = 0; i <
grid->ndim; i++ )
3521 vector<realnum>& flux1,
3531 lgDryRun = ( flux1.size() == 0 );
3537 ind = (
grid->jlo[n] >= 0 ) ?
grid->jlo[n] :
grid->jhi[n];
3539 ind = (
grid->jhi[n] >= 0 ) ?
grid->jhi[n] :
grid->jlo[n];
3540 else if(
grid->ndim == 1 )
3548 fprintf(
ioQQQ,
" The requested interpolation could not be completed, sorry.\n" );
3549 fprintf(
ioQQQ,
" No suitable match was found for a model with" );
3550 for( i=0; i <
grid->ndim; i++ )
3551 fprintf(
ioQQQ,
" %s=%.6g ",
grid->names[i],
grid->val[i][index[i]] );
3552 fprintf(
ioQQQ,
"\n" );
3556 for( i=0; i <
grid->npar; i++ )
3557 aval[i] =
grid->telg[ind].par[i];
3561 for( i=0; i <
grid->ndim &&
called.lgTalk; i++ )
3565 fprintf(
ioQQQ,
" No exact match was found for a model with" );
3566 for( j=0; j <
grid->ndim; j++ )
3567 fprintf(
ioQQQ,
" %s=%.6g ",
grid->names[j],
grid->val[j][index[j]] );
3568 fprintf(
ioQQQ,
"- using the following model instead:\n" );
3578 vector<realnum> flux2(
rfield.nupper);
3582 const realnum SECURE = 10.f*FLT_EPSILON;
3585 aval2 = (
double*)
MALLOC((
size_t)(
grid->npar*
sizeof(
double)) );
3596 index[nd] = indlo[nd];
3602 index[nd] = indhi[nd];
3603 vector<realnum> empty;
3606 if( !
fp_equal(aval2[nd],aval[nd],10) )
3608 double fr1, fr2, fc1 = 0., fc2 = 0.;
3613 fr1 = (aval2[nd]-val[nd])/(aval2[nd]-aval[nd]);
3621 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3626 fprintf(
ioQQQ,
"interpolation nd=%ld fr1=%g\n", nd, fr1 );
3630 if( nd == 0 && strcmp(
grid->names[nd],
"Teff" ) == 0 )
3642 fc1 = ( val[nd] > 200000. ) ? log10(val[nd]/
grid->val[nd][indlo[nd]])*4. : 0.;
3643 fc2 = ( val[nd] > 200000. ) ? log10(val[nd]/
grid->val[nd][indhi[nd]])*4. : 0.;
3646 for( i=0; i <
rfield.nupper; ++i )
3647 flux1[i] = (
realnum)(fr1*(flux1[i]+fc1) + fr2*(flux2[i]+fc2));
3650 for( i=0; i <
grid->npar; i++ )
3651 aval[i] = fr1*aval[i] + fr2*aval2[i];
3667 vector<realnum>& flux1)
3675 ind = ( index[1] == 0 ) ? indlo[index[0]] : indhi[index[0]];
3679 for( i=0; i <
grid->npar; i++ )
3680 aval[i] =
grid->telg[ind].par[i];
3686 const realnum SECURE = 10.f*FLT_EPSILON;
3697 lgSkip = ( nd == 1 ) ? ( indhi[index[0]] == indlo[index[0]] ) :
3698 ( indlo[0] == indlo[1] && indhi[0] == indhi[1] );
3702 vector<realnum> flux2(
rfield.nupper);
3703 double fr1, fr2, *aval2;
3705 aval2 = (
double*)
MALLOC((
size_t)(
grid->npar*
sizeof(
double)) );
3710 fr1 = (aval2[nd+off]-val[nd])/(aval2[nd+off]-aval[nd+off]);
3714 fprintf(
ioQQQ,
"interpolation nd=%ld fr1=%g\n", nd, fr1 );
3717 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3719 for( i=0; i <
rfield.nupper; ++i )
3720 flux1[i] = (
realnum)(fr1*flux1[i] + fr2*flux2[i]);
3722 for( i=0; i <
grid->npar; i++ )
3723 aval[i] = fr1*aval[i] + fr2*aval2[i];
3732 vector<Energy>& ener)
3743 if( fseek(
grid->ioIN, (
long)(
grid->nOffset), SEEK_SET ) != 0 )
3745 fprintf(
ioQQQ,
" Error finding atmosphere frequency bins\n");
3749 vector<realnum> data(
rfield.nupper);
3752 fprintf(
ioQQQ,
" Error reading atmosphere frequency bins\n" );
3756 for(
long i=0; i <
rfield.nupper; ++i )
3757 ener[i].set(data[i]);
3764 vector<realnum>& flux,
3778 ASSERT( ind >= 0 && ind <= grid->nmods );
3782 if( fseek(
grid->ioIN, (
long)(ind*
grid->nBlocksize+
grid->nOffset), SEEK_SET ) != 0 )
3784 fprintf(
ioQQQ,
" Error seeking atmosphere %ld\n", ind );
3790 fprintf(
ioQQQ,
" Error trying to read atmosphere %ld\n", ind );
3795 if(
called.lgTalk && lgTalk )
3798 if(
grid->npar == 1 )
3800 " * c<< %s model%5ld read. "
3801 " %6s = %13.2f >>> *\n",
3802 grid->ident, ind,
grid->names[0],
grid->telg[ind-1].par[0] );
3803 else if(
grid->npar == 2 )
3805 " * c<< %s model%5ld read. "
3806 " %6s = %10.2f %6s = %8.5f >>> *\n",
3807 grid->ident, ind,
grid->names[0],
grid->telg[ind-1].par[0],
3808 grid->names[1],
grid->telg[ind-1].par[1] );
3809 else if(
grid->npar == 3 )
3811 " * c<< %s model%5ld read. "
3812 " %6s=%7.0f %6s=%5.2f %6s=%5.2f >>> *\n",
3813 grid->ident, ind,
grid->names[0],
grid->telg[ind-1].par[0],
3814 grid->names[1],
grid->telg[ind-1].par[1],
3815 grid->names[2],
grid->telg[ind-1].par[2] );
3816 else if(
grid->npar >= 4 )
3820 " %4s=%5.0f %6s=%4.2f %6s=%5.2f %6s=",
3821 grid->ident, ind,
grid->names[0],
grid->telg[ind-1].par[0],
3822 grid->names[1],
grid->telg[ind-1].par[1],
3823 grid->names[2],
grid->telg[ind-1].par[2],
grid->names[3] );
3825 fprintf(
ioQQQ,
" >> *\n" );
3832 for( i=0; i <
rfield.nupper; ++i )
3837 volatile double help = flux[i];
3862 long index[
MDIM], j;
3863 const double SECURE = (1. + 20.*(double)FLT_EPSILON);
3868 switch(
grid->imode )
3878 for( j=0; j <
grid->nTracks; j++ )
3880 if( ValTr[j] != -FLT_MAX )
3884 pow(10.,(
double)ValTr[j]) : ValTr[j];
3885 *loLim =
MIN2(*loLim,temp);
3886 *hiLim =
MAX2(*hiLim,temp);
3892 index[1] = useTr[0];
3894 index[1] = useTr[1];
3896 *loLim =
MAX2(
grid->telg[ptr0].par[3],
grid->telg[ptr1].par[3]);
3898 printf(
"set limit 0: (models %d, %d) %f %f\n",
3899 ptr0+1, ptr1+1,
grid->telg[ptr0].par[3],
grid->telg[ptr1].par[3] );
3901 index[0] =
grid->trackLen[useTr[0]]-1;
3902 index[1] = useTr[0];
3904 index[0] =
grid->trackLen[useTr[1]]-1;
3905 index[1] = useTr[1];
3907 *hiLim =
MIN2(
grid->telg[ptr0].par[3],
grid->telg[ptr1].par[3]);
3909 printf(
"set limit 1: (models %d, %d) %f %f\n",
3910 ptr0+1, ptr1+1,
grid->telg[ptr0].par[3],
grid->telg[ptr1].par[3] );
3914 fprintf(
ioQQQ,
" SetLimits called with insane value for imode: %d.\n",
grid->imode );
3918 ASSERT( fabs(*loLim) < DBL_MAX && fabs(*hiLim) < DBL_MAX );
3921 if( *hiLim <= *loLim )
3923 fprintf(
ioQQQ,
" no room to optimize: lower limit %.4f, upper limit %.4f.\n",
3933 printf(
"set limits: %g %g\n",*loLim,*hiLim);
3961 double loLoc = +DBL_MAX;
3962 double hiLoc = -DBL_MAX;
3964 for( index[0]=0; index[0] <
grid->nval[0]; ++index[0] )
3973 if(
grid->jlo[n] < 0 &&
grid->jhi[n] < 0 )
3977 if(
grid->val[0][index[0]] < val )
3980 if(
grid->val[0][index[0]] > val )
3987 if(
grid->val[0][index[0]] <= val )
3991 if( loLoc == DBL_MAX )
3992 loLoc =
grid->val[0][index[0]];
3994 if(
grid->val[0][index[0]] >= val )
3998 hiLoc =
grid->val[0][index[0]];
4003 ASSERT( fabs(loLoc) < DBL_MAX && fabs(hiLoc) < DBL_MAX && loLoc <= hiLoc );
4005 *loLim =
MAX2(*loLim,loLoc);
4006 *hiLim =
MIN2(*hiLim,hiLoc);
4010 index[nd] = indlo[nd];
4013 if( indhi[nd] != indlo[nd] )
4015 index[nd] = indhi[nd];
4025 long i, *index, j, jsize, nd;
4036 for( nd=0; nd <
grid->ndim; nd++ )
4038 double pval =
grid->telg[0].par[nd];
4039 grid->val[nd][0] = pval;
4042 for( i=1; i <
grid->nmods; i++ )
4047 pval =
grid->telg[i].par[nd];
4057 for( j =
grid->nval[nd]-1; j >= i2; j-- )
4058 grid->val[nd][j+1] =
grid->val[nd][j];
4059 grid->val[nd][i2] = pval;
4064 jsize *=
grid->nval[nd];
4067 printf(
"%s[%ld]:",
grid->names[nd],
grid->nval[nd] );
4068 for( i=0; i <
grid->nval[nd]; i++ )
4069 printf(
" %g",
grid->val[nd][i] );
4074 index = (
long *)
MALLOC(
sizeof(
long)*
grid->ndim );
4075 val = (
double *)
MALLOC(
sizeof(
double)*
grid->ndim );
4078 grid->jlo = (
long *)
MALLOC(
sizeof(
long)*jsize );
4079 grid->jhi = (
long *)
MALLOC(
sizeof(
long)*jsize );
4111 for( index[nd]=0; index[nd] <
grid->nval[nd]; index[nd]++ )
4113 val[nd] =
grid->val[nd][index[nd]];
4118 if( lgList && nd ==
MIN2(
grid->ndim-1,1) )
4120 fprintf(
ioQQQ,
"\n" );
4121 if(
grid->ndim > 2 )
4123 fprintf(
ioQQQ,
"subgrid for" );
4124 for(
long n = nd+1; n <
grid->ndim; n++ )
4125 fprintf(
ioQQQ,
" %s=%g",
grid->names[n], val[n] );
4126 fprintf(
ioQQQ,
":\n\n" );
4128 if(
grid->ndim > 1 )
4130 fprintf(
ioQQQ,
"%6.6s\\%6.6s |",
grid->names[0],
grid->names[1] );
4131 for(
long n = 0; n <
grid->nval[1]; n++ )
4132 fprintf(
ioQQQ,
" %9.3g",
grid->val[1][n] );
4133 fprintf(
ioQQQ,
"\n" );
4134 fprintf(
ioQQQ,
"--------------|" );
4135 for(
long n = 0; n <
grid->nval[1]; n++ )
4136 fprintf(
ioQQQ,
"----------" );
4140 fprintf(
ioQQQ,
"%13.13s |\n",
grid->names[0] );
4141 fprintf(
ioQQQ,
"--------------|----------" );
4143 fprintf(
ioQQQ,
"\n" );
4144 for( index[0]=0; index[0] <
grid->nval[0]; index[0]++ )
4146 fprintf(
ioQQQ,
"%13.7g |",
grid->val[0][index[0]] );
4147 if(
grid->ndim > 1 )
4149 for( index[1]=0; index[1] <
grid->nval[1]; index[1]++ )
4154 fprintf(
ioQQQ,
" --" );
4160 fprintf(
ioQQQ,
"\n" );
4162 fprintf(
ioQQQ,
"\n" );
4176 for( i=0; i <
grid->ndim; i++ )
4178 ind += index[i]*mul;
4179 mul *=
grid->nval[i];
4185 bool lgIsTeffLoggGrid,
4193 double alogg_low = -DBL_MAX, alogg_high = DBL_MAX;
4211 *index_low = *index_high = -2;
4212 for( i=0; i < nmods; i++ )
4214 bool lgNext =
false;
4216 for( nd=0; nd < ndim; nd++ )
4218 if( nd != 1 && !
fp_equal(telg[i].par[nd],val[nd],10) )
4228 if( ndim == 1 ||
fp_equal(telg[i].par[1],val[1],10) )
4234 if( lgIsTeffLoggGrid )
4237 if( telg[i].par[1] < val[1] && telg[i].par[1] > alogg_low )
4240 alogg_low = telg[i].
par[1];
4243 if( telg[i].par[1] > val[1] && telg[i].par[1] < alogg_high )
4246 alogg_high = telg[i].
par[1];
4260 bool lgOutLo, lgOutHi;
4280 lgOutLo = ( x-xval[0] < -10.*DBL_EPSILON*fabs(xval[0]) );
4281 lgOutHi = ( x-xval[NVAL-1] > 10.*DBL_EPSILON*fabs(xval[NVAL-1]) );
4283 if( lgOutLo || lgOutHi )
4289 *ind1 = lgOutLo ? -1 : NVAL-1;
4290 *ind2 = lgOutLo ? 0 : NVAL;
4302 for( i=0; i < NVAL; i++ )
4313 for( i=0; i < NVAL-1; i++ )
4315 if( xval[i] < x && x < xval[i+1] )
4324 fprintf(
ioQQQ,
" insanity in FindIndex\n" );
4335 ioIN =
open_data( chFnam,
"r", scheme );
4353 vector<Energy> anu(
rfield.nupper);
4354 vector<realnum> flux(
rfield.nupper);
4358 if( strcmp(
grid->names[0],
"Teff" ) != 0 )
4365 for( i=0; i <
grid->nmods; i++ )
4367 fprintf(
ioQQQ,
"testing model %ld ", i+1 );
4368 for( nd=0; nd <
grid->npar; nd++ )
4369 fprintf(
ioQQQ,
" %s %g",
grid->names[nd],
grid->telg[i].par[nd] );
4374 fprintf(
ioQQQ,
" OK\n" );
4378 FILE *ioBUG = fopen(
"atmosphere_dump.txt", ( i == 0 ) ?
"w" :
"a" );
4380 fprintf( ioBUG,
"######## MODEL %ld", i+1 );
4381 for( nd=0; nd <
grid->npar; nd++ )
4382 fprintf( ioBUG,
" %s %g",
grid->names[nd],
grid->telg[i].par[nd] );
4383 fprintf( ioBUG,
"####################\n" );
4385 for( k=0; k <
rfield.nupper; ++k )
4386 fprintf( ioBUG,
"%e %e\n", anu[k].Ryd(), flux[k] );
4395 const vector<realnum>& flux,
4399 bool lgPassed =
true;
4409 for( k=1; k <
rfield.nupper; k++ )
4410 lumi += (anu[k].Ryd() - anu[k-1].Ryd())*(flux[k] + flux[k-1])/2.;
4415 if( fabs(Teff - chk) > toler*Teff ) {
4416 fprintf(
ioQQQ,
"\n*** WARNING, Teff discrepancy for this model, expected Teff %.2f, ", Teff);
4417 fprintf(
ioQQQ,
"integration yielded Teff %.2f, delta %.2f%%\n", chk, (chk/Teff-1.)*100. );
4454 for( j=0; j < nEdge; j++ )
4456 ind =
RebinFind(StarEner,nCont,AbsorbEdge[j]);
4459 ASSERT( ind >= 0 && ind+1 < nCont );
4461 EdgeLow[j] = StarEner[ind];
4462 EdgeHigh[j] = StarEner[ind+1];
4468 for( j=0; j < nCont; j++ )
4470 if( StarFlux[j] == 0.f )
4480 for( j=0; j < nCont-1; j++ )
4482 double ratio_x, ratio_y;
4485 ASSERT( StarEner[j+1] > StarEner[j] );
4489 ratio_x = (double)StarEner[j+1]/(
double)StarEner[j];
4490 ratio_y = (double)StarFlux[j+1]/(
double)StarFlux[j];
4491 StarPower[j] = (
realnum)(log(ratio_y)/log(ratio_x));
4494 for( j=0; j <
rfield.nupper; j++ )
4498 BinLow = ( j > 0 ) ?
4502 BinHigh = ( j+1 <
rfield.nupper ) ?
4506 BinNext = ( j+2 <
rfield.nupper ) ?
4514 for( k=0; k < nEdge; k++ )
4516 if( BinLow < EdgeLow[k] && BinNext > EdgeHigh[k] )
4518 BinMid = 0.99999f*EdgeLow[k];
4519 CloudyFlux[j] =
RebinSingleCell(BinLow,BinMid,StarEner,StarFlux,StarPower,nCont);
4525 BinMid = 1.00001f*EdgeHigh[k];
4526 CloudyFlux[j] =
RebinSingleCell(BinMid,BinNext,StarEner,StarFlux,StarPower,nCont);
4535 CloudyFlux[j] =
RebinSingleCell(BinLow,BinHigh,StarEner,StarFlux,StarPower,nCont);
4567 anu = sqrt(BinLow*BinHigh);
4569 widflx =
MIN2(BinHigh,StarEner[nCont-1])-BinLow;
4571 if( BinLow < StarEner[0] )
4575 retval = (
realnum)(StarFlux[0]*pow(anu/StarEner[0],2.));
4577 else if( BinLow > StarEner[nCont-1] )
4586 ipLo =
RebinFind(StarEner,nCont,BinLow);
4587 ipHi =
RebinFind(StarEner,nCont,BinHigh);
4589 ASSERT( ipLo >= 0 && ipLo < nCont-1 && ipHi >= ipLo );
4595 retval = (
realnum)(StarFlux[ipLo]*pow(anu/StarEner[ipLo],(
double)StarPower[ipLo]));
4609 for( i=ipLo; i <=
MIN2(ipHi,nCont-2); i++ )
4611 double pp1 = StarPower[i] + 1.;
4617 v1 = StarFlux[i]*pow(
x1/StarEner[i],(
double)StarPower[i]);
4621 else if( i == ipHi )
4638 if( fabs(pp1) < 0.001 )
4644 val = pow(
x2/
x1,pp1) - 1.;
4645 val = val*
x1*v1/pp1;
4650 retval = sum/widflx;
4678 if( val < array[0] )
4682 else if( val > array[nArr-1] )
4694 sgn =
sign3(val-array[i2]);
const int INPUT_LINE_LENGTH
bool fp_equal(sys_float x, sys_float y, int n=3)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
const ios_base::openmode mode_r
UNUSED const double FR1RYD
UNUSED const double STEFAN_BOLTZ
UNUSED const double SPEEDLIGHT
UNUSED const double RYDLAM
static const bool lgLINEAR
STATIC void InitGridCoStar(stellar_grid *)
STATIC void InterpolateModel(const stellar_grid *, const double[], double[], const long[], const long[], long[], long, vector< realnum > &, IntStage)
int Kurucz79Compile(process_counter &pc)
bool StarburstCompile(process_counter &pc)
static const bool lgTAKELOG
bool StarburstInitialize(const char chInName[], const char chOutName[], sb_mode mode)
STATIC void FindVCoStar(const stellar_grid *, double, realnum *, long[])
STATIC void InterpolateRectGrid(const stellar_grid *, const double[], double *, double *)
STATIC void FindHCoStar(const stellar_grid *, long, double, long, realnum *, long *, long *)
STATIC bool lgValidBinFile(const char *, process_counter &, access_scheme)
STATIC void InterpolateGridCoStar(const stellar_grid *, const double[], double *, double *)
long RauchInterpolateHydr(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC void SetLimits(const stellar_grid *, double, const long[], const long[], const long[], const realnum[], double *, double *)
static const long int VERSION_BIN
STATIC bool lgValidModel(const vector< Energy > &, const vector< realnum > &, double, double)
long Kurucz79Interpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
static const bool lgSILENT
int WMBASICCompile(process_counter &pc)
long WMBASICInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC void SetLimitsSub(const stellar_grid *, double, const long[], const long[], long[], long, double *, double *)
STATIC void RebinAtmosphere(long, const realnum[], const realnum[], realnum[], long, const realnum[])
int AtlasCompile(process_counter &pc)
STATIC void CoStarListModels(const stellar_grid *)
long AtlasInterpolate(double val[], long *nval, long *ndim, const char *chMetalicity, const char *chODFNew, bool lgList, double *Tlow, double *Thigh)
int RauchCompile(process_counter &pc)
STATIC void FindIndex(const double[], long, double, long *, long *, bool *)
static const int NMODS_HpHE
bool GridCompile(const char *InName)
long WernerInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC int RauchInitializeSub(const char[], const char[], const vector< mpp > &, long, long, long, const double[], int)
void getdataline(fstream &, string &)
STATIC void FillJ(const stellar_grid *, long[], double[], long, bool)
STATIC bool lgValidAsciiFile(const char *, access_scheme)
static const int NMODS_HYDR
int MihalasCompile(process_counter &pc)
STATIC void GetBins(const stellar_grid *, vector< Energy > &)
STATIC bool lgCompileAtmosphereCoStar(const char[], const char[], const realnum[], long, process_counter &)
static const int NMODS_PG1159
static const long int VERSION_ASCII
long RauchInterpolateHNi(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
STATIC bool lgFileReadable(const char *, process_counter &, access_scheme)
STATIC void ValidateGrid(const stellar_grid *, double)
STATIC void RauchReadMPP(vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &)
int CoStarCompile(process_counter &pc)
long CoStarInterpolate(double val[], long *nval, long *ndim, IntMode imode, bool lgHalo, bool lgList, double *val0_lo, double *val0_hi)
static const long int VERSION_RAUCH_MPP
STATIC void InitGrid(stellar_grid *, bool)
static const int NMODS_HCA
static const int NMODS_HNI
STATIC long JIndex(const stellar_grid *, const long[])
long GridInterpolate(double val[], long *nval, long *ndim, const char *FileName, bool lgList, double *Tlow, double *Thigh)
STATIC bool lgCompileAtmosphere(const char[], const char[], const realnum[], long, process_counter &)
long RauchInterpolateHpHe(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
int TlustyCompile(process_counter &pc)
long RauchInterpolateHelium(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
long TlustyInterpolate(double val[], long *nval, long *ndim, tl_grid tlg, const char *chMetalicity, bool lgList, double *Tlow, double *Thigh)
STATIC realnum RebinSingleCell(realnum, realnum, const realnum[], const realnum[], const realnum[], long)
long RauchInterpolateHCa(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
int WernerCompile(process_counter &pc)
STATIC void CheckVal(const stellar_grid *, double[], long *, long *)
static const int NMODS_HELIUM
STATIC void InterpolateModelCoStar(const stellar_grid *, const double[], double[], const long[], const long[], long[], long, long, vector< realnum > &)
STATIC void GetModel(const stellar_grid *, long, vector< realnum > &, bool, bool)
void AtmospheresAvail(void)
long RauchInterpolateCOWD(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
static const bool lgVERBOSE
STATIC void SearchModel(const mpp[], bool, long, const double[], long, long *, long *)
STATIC void InitIndexArrays(stellar_grid *, bool)
long MihalasInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC long RebinFind(const realnum[], long, realnum)
STATIC void FreeGrid(stellar_grid *)
long RauchInterpolatePG1159(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
static const unsigned int NMD5