16 static const int MNTS = 200;
36 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
37 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
149 long,
const double[],
int);
161 const long[],
long[],
long,
long,
realnum[]);
164 const realnum[],
double*,
double*);
201 fprintf(
ioQQQ,
"\n I will now list all stellar atmosphere grids that are ready to be used (if any).\n" );
202 fprintf(
ioQQQ,
" User-defined stellar atmosphere grids will not be included in this list.\n\n" );
211 fprintf(
ioQQQ,
" table star atlas Z+1.0 <Teff> [ <log(g)> ]\n" );
213 fprintf(
ioQQQ,
" table star atlas Z+0.5 <Teff> [ <log(g)> ]\n" );
215 fprintf(
ioQQQ,
" table star atlas Z+0.3 <Teff> [ <log(g)> ]\n" );
217 fprintf(
ioQQQ,
" table star atlas Z+0.2 <Teff> [ <log(g)> ]\n" );
219 fprintf(
ioQQQ,
" table star atlas Z+0.1 <Teff> [ <log(g)> ]\n" );
221 fprintf(
ioQQQ,
" table star atlas Z+0.0 <Teff> [ <log(g)> ]\n" );
223 fprintf(
ioQQQ,
" table star atlas Z-0.1 <Teff> [ <log(g)> ]\n" );
225 fprintf(
ioQQQ,
" table star atlas Z-0.2 <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.5 <Teff> [ <log(g)> ]\n" );
231 fprintf(
ioQQQ,
" table star atlas Z-1.0 <Teff> [ <log(g)> ]\n" );
233 fprintf(
ioQQQ,
" table star atlas Z-1.5 <Teff> [ <log(g)> ]\n" );
235 fprintf(
ioQQQ,
" table star atlas Z-2.0 <Teff> [ <log(g)> ]\n" );
237 fprintf(
ioQQQ,
" table star atlas Z-2.5 <Teff> [ <log(g)> ]\n" );
239 fprintf(
ioQQQ,
" table star atlas Z-3.0 <Teff> [ <log(g)> ]\n" );
241 fprintf(
ioQQQ,
" table star atlas Z-3.5 <Teff> [ <log(g)> ]\n" );
243 fprintf(
ioQQQ,
" table star atlas Z-4.0 <Teff> [ <log(g)> ]\n" );
245 fprintf(
ioQQQ,
" table star atlas Z-4.5 <Teff> [ <log(g)> ]\n" );
247 fprintf(
ioQQQ,
" table star atlas Z-5.0 <Teff> [ <log(g)> ]\n" );
250 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.5 <Teff> [ <log(g)> ]\n" );
252 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.2 <Teff> [ <log(g)> ]\n" );
254 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.0 <Teff> [ <log(g)> ]\n" );
256 fprintf(
ioQQQ,
" table star atlas odfnew Z-0.5 <Teff> [ <log(g)> ]\n" );
258 fprintf(
ioQQQ,
" table star atlas odfnew Z-1.0 <Teff> [ <log(g)> ]\n" );
260 fprintf(
ioQQQ,
" table star atlas odfnew Z-1.5 <Teff> [ <log(g)> ]\n" );
262 fprintf(
ioQQQ,
" table star atlas odfnew Z-2.0 <Teff> [ <log(g)> ]\n" );
264 fprintf(
ioQQQ,
" table star atlas odfnew Z-2.5 <Teff> [ <log(g)> ]\n" );
267 fprintf(
ioQQQ,
" table star atlas 3-dim <Teff> <log(g)> <log(Z)>\n" );
270 fprintf(
ioQQQ,
" table star atlas odfnew 3-dim <Teff> <log(g)> <log(Z)>\n" );
273 fprintf(
ioQQQ,
" table star costar solar (see Hazy for parameters)\n" );
275 fprintf(
ioQQQ,
" table star costar halo (see Hazy for parameters)\n" );
278 fprintf(
ioQQQ,
" table star kurucz79 <Teff>\n" );
281 fprintf(
ioQQQ,
" table star mihalas <Teff>\n" );
284 fprintf(
ioQQQ,
" table star rauch H-Ca solar <Teff> [ <log(g)> ]\n" );
286 fprintf(
ioQQQ,
" table star rauch H-Ca halo <Teff> [ <log(g)> ]\n" );
288 fprintf(
ioQQQ,
" table star rauch H-Ca 3-dim <Teff> <log(g)> <log(Z)>\n" );
291 fprintf(
ioQQQ,
" table star rauch H-Ni solar <Teff> [ <log(g)> ]\n" );
293 fprintf(
ioQQQ,
" table star rauch H-Ni halo <Teff> [ <log(g)> ]\n" );
295 fprintf(
ioQQQ,
" table star rauch H-Ni 3-dim <Teff> <log(g)> <log(Z)>\n" );
298 fprintf(
ioQQQ,
" table star rauch pg1159 <Teff> [ <log(g)> ]\n" );
301 fprintf(
ioQQQ,
" table star rauch hydrogen <Teff> [ <log(g)> ]\n" );
304 fprintf(
ioQQQ,
" table star rauch helium <Teff> [ <log(g)> ]\n" );
307 fprintf(
ioQQQ,
" table star rauch H+He <Teff> <log(g)> <frac(He)>\n" );
310 fprintf(
ioQQQ,
" table star \"starburst99.mod\" <age>\n" );
313 fprintf(
ioQQQ,
" table star tlusty Bstar Z+0.3 <Teff> [ <log(g)> ]\n" );
315 fprintf(
ioQQQ,
" table star tlusty Bstar Z+0.0 <Teff> [ <log(g)> ]\n" );
317 fprintf(
ioQQQ,
" table star tlusty Bstar Z-0.3 <Teff> [ <log(g)> ]\n" );
319 fprintf(
ioQQQ,
" table star tlusty Bstar Z-0.7 <Teff> [ <log(g)> ]\n" );
321 fprintf(
ioQQQ,
" table star tlusty Bstar Z-1.0 <Teff> [ <log(g)> ]\n" );
323 fprintf(
ioQQQ,
" table star tlusty Bstar Z-inf <Teff> [ <log(g)> ]\n" );
326 fprintf(
ioQQQ,
" table star tlusty Bstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
329 fprintf(
ioQQQ,
" table star tlusty Ostar Z+0.3 <Teff> [ <log(g)> ]\n" );
331 fprintf(
ioQQQ,
" table star tlusty Ostar Z+0.0 <Teff> [ <log(g)> ]\n" );
333 fprintf(
ioQQQ,
" table star tlusty Ostar Z-0.3 <Teff> [ <log(g)> ]\n" );
335 fprintf(
ioQQQ,
" table star tlusty Ostar Z-0.7 <Teff> [ <log(g)> ]\n" );
337 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.0 <Teff> [ <log(g)> ]\n" );
339 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.5 <Teff> [ <log(g)> ]\n" );
341 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.7 <Teff> [ <log(g)> ]\n" );
343 fprintf(
ioQQQ,
" table star tlusty Ostar Z-2.0 <Teff> [ <log(g)> ]\n" );
345 fprintf(
ioQQQ,
" table star tlusty Ostar Z-3.0 <Teff> [ <log(g)> ]\n" );
347 fprintf(
ioQQQ,
" table star tlusty Ostar Z-inf <Teff> [ <log(g)> ]\n" );
350 fprintf(
ioQQQ,
" table star tlusty Ostar 3-dim <Teff> <log(g)> <log(Z)>\n" );
353 fprintf(
ioQQQ,
" table star werner <Teff> [ <log(g)> ]\n" );
356 fprintf(
ioQQQ,
" table star wmbasic <Teff> <log(g)> <log(Z)>\n" );
379 fprintf(
ioQQQ,
" AtlasCompile on the job.\n" );
401 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp10k2.ascii",
"atlas_fp10k2.mod", Edges, 3L, pc );
403 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp05k2.ascii",
"atlas_fp05k2.mod", Edges, 3L, pc );
405 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp03k2.ascii",
"atlas_fp03k2.mod", Edges, 3L, pc );
407 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp02k2.ascii",
"atlas_fp02k2.mod", Edges, 3L, pc );
409 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp01k2.ascii",
"atlas_fp01k2.mod", Edges, 3L, pc );
411 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp00k2.ascii",
"atlas_fp00k2.mod", Edges, 3L, pc );
413 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm01k2.ascii",
"atlas_fm01k2.mod", Edges, 3L, pc );
415 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm02k2.ascii",
"atlas_fm02k2.mod", Edges, 3L, pc );
417 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm03k2.ascii",
"atlas_fm03k2.mod", Edges, 3L, pc );
419 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm05k2.ascii",
"atlas_fm05k2.mod", Edges, 3L, pc );
421 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm10k2.ascii",
"atlas_fm10k2.mod", Edges, 3L, pc );
423 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm15k2.ascii",
"atlas_fm15k2.mod", Edges, 3L, pc );
425 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm20k2.ascii",
"atlas_fm20k2.mod", Edges, 3L, pc );
427 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm25k2.ascii",
"atlas_fm25k2.mod", Edges, 3L, pc );
429 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm30k2.ascii",
"atlas_fm30k2.mod", Edges, 3L, pc );
431 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm35k2.ascii",
"atlas_fm35k2.mod", Edges, 3L, pc );
433 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm40k2.ascii",
"atlas_fm40k2.mod", Edges, 3L, pc );
435 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm45k2.ascii",
"atlas_fm45k2.mod", Edges, 3L, pc );
437 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm50k2.ascii",
"atlas_fm50k2.mod", Edges, 3L, pc );
442 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp05k2_odfnew.ascii",
"atlas_fp05k2_odfnew.mod",
446 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp02k2_odfnew.ascii",
"atlas_fp02k2_odfnew.mod",
450 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp00k2_odfnew.ascii",
"atlas_fp00k2_odfnew.mod",
454 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm05k2_odfnew.ascii",
"atlas_fm05k2_odfnew.mod",
458 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm10k2_odfnew.ascii",
"atlas_fm10k2_odfnew.mod",
462 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm15k2_odfnew.ascii",
"atlas_fm15k2_odfnew.mod",
466 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm20k2_odfnew.ascii",
"atlas_fm20k2_odfnew.mod",
470 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm25k2_odfnew.ascii",
"atlas_fm25k2_odfnew.mod",
478 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_3d_odfnew.ascii",
"atlas_3d_odfnew.mod", Edges, 3L, pc );
486 const char *chMetalicity,
487 const char *chODFNew,
497 grid.
name =
"atlas_";
503 grid.
name += chMetalicity;
506 grid.
name += chODFNew;
513 strcpy( chIdent,
"3-dim" );
517 strcpy( chIdent,
"Z " );
518 strcat( chIdent, chMetalicity );
520 strcat( chIdent, ( strlen(chODFNew) == 0 ?
" Kurucz" :
" ODFNew" ) );
521 grid.
ident = chIdent;
523 grid.
command =
"COMPILE STARS";
561 fprintf(
ioQQQ,
" CoStarCompile on the job.\n" );
604 grid.
name = ( lgHalo ?
"Sc1_costar_halo.mod" :
"Sc1_costar_solar.mod" );
608 grid.
ident =
" costar";
610 grid.
command =
"COMPILE STARS";
665 string OutName( InName );
670 fprintf(
ioQQQ,
" GridCompile on the job.\n" );
673 string::size_type ptr = OutName.find(
'.' );
674 ASSERT( ptr != string::npos );
675 OutName.replace( ptr, string::npos,
".mod" );
685 grid.
ident =
"bogus ident.";
686 grid.
command =
"bogus command.";
692 fprintf(
ioQQQ,
" GridCompile: checking effective temperatures...\n" );
705 const char *FileName,
716 string chTruncName( FileName );
717 string::size_type ptr = chTruncName.find(
'.' );
718 if( ptr != string::npos )
719 chTruncName.replace( ptr, string::npos,
"" );
721 grid.
name = FileName;
725 sprintf( chIdent,
"%12.12s", chTruncName.c_str() );
726 grid.
ident = chIdent;
728 string chString(
"COMPILE STARS \"" + chTruncName +
".ascii\"" );
729 grid.
command = chString.c_str();
750 fprintf(
ioQQQ,
" Kurucz79Compile on the job.\n" );
774 grid.
name =
"kurucz79.mod";
778 grid.
ident =
" Kurucz 1979";
780 grid.
command =
"COMPILE STARS";
801 fprintf(
ioQQQ,
" MihalasCompile on the job.\n" );
824 grid.
name =
"mihalas.mod";
828 grid.
ident =
" Mihalas";
830 grid.
command =
"COMPILE STARS";
859 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}},
860 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}},
861 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}},
862 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}},
863 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}},
864 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},
865 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},
866 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},
867 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},
868 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},
869 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},
870 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},
871 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},
872 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},
873 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},
874 {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}},
875 {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}},
876 {{400000.,8.0}},{{400000.,9.0}},
877 {{500000.,8.0}},{{500000.,9.0}},
886 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}},
887 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}},
888 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}},
889 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}},
890 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}},
891 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},
892 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},
893 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},
894 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},
895 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},
896 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},
897 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},
898 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},
899 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},
900 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}}
904 {{40000.,5.0}}, {{40000.,6.0}}, {{40000.,7.0}}, {{40000.,8.0}}, {{40000.,9.0}},
905 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}},
906 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}},
907 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}},
908 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}},
909 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}},
910 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}},
911 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}},
912 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}},
913 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}},
914 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}},
915 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}},
916 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}},
917 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}},
918 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}},
919 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}}
923 {{20000.,4.0}}, {{20000.,5.0}}, {{20000.,6.0}}, {{20000.,7.0}}, {{20000.,8.0}}, {{20000.,9.0}},
924 {{30000.,4.0}}, {{30000.,5.0}}, {{30000.,6.0}}, {{30000.,7.0}}, {{30000.,8.0}}, {{30000.,9.0}},
925 {{40000.,4.0}}, {{40000.,5.0}}, {{40000.,6.0}}, {{40000.,7.0}}, {{40000.,8.0}}, {{40000.,9.0}},
926 {{50000.,4.0}}, {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}},
927 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}},
928 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}},
929 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}},
930 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}},
931 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}},
932 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}},
933 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}},
934 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}},
935 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}},
936 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}},
937 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}},
938 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}},
939 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}},
940 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}},
941 {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}},
942 {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}},
943 {{400000.,8.0}},{{400000.,9.0}},
944 {{500000.,8.0}},{{500000.,9.0}},
953 {{50000.,5.0}}, {{50000.,6.0}}, {{50000.,7.0}}, {{50000.,8.0}}, {{50000.,9.0}},
954 {{60000.,5.0}}, {{60000.,6.0}}, {{60000.,7.0}}, {{60000.,8.0}}, {{60000.,9.0}},
955 {{70000.,5.0}}, {{70000.,6.0}}, {{70000.,7.0}}, {{70000.,8.0}}, {{70000.,9.0}},
956 {{80000.,5.0}}, {{80000.,6.0}}, {{80000.,7.0}}, {{80000.,8.0}}, {{80000.,9.0}},
957 {{90000.,5.0}}, {{90000.,6.0}}, {{90000.,7.0}}, {{90000.,8.0}}, {{90000.,9.0}},
958 {{100000.,5.0}},{{100000.,6.0}},{{100000.,7.0}},{{100000.,8.0}},{{100000.,9.0}},
959 {{110000.,6.0}},{{110000.,7.0}},{{110000.,8.0}},{{110000.,9.0}},
960 {{120000.,6.0}},{{120000.,7.0}},{{120000.,8.0}},{{120000.,9.0}},
961 {{130000.,6.0}},{{130000.,7.0}},{{130000.,8.0}},{{130000.,9.0}},
962 {{140000.,6.0}},{{140000.,7.0}},{{140000.,8.0}},{{140000.,9.0}},
963 {{150000.,6.0}},{{150000.,7.0}},{{150000.,8.0}},{{150000.,9.0}},
964 {{160000.,6.0}},{{160000.,7.0}},{{160000.,8.0}},{{160000.,9.0}},
965 {{170000.,6.0}},{{170000.,7.0}},{{170000.,8.0}},{{170000.,9.0}},
966 {{180000.,6.0}},{{180000.,7.0}},{{180000.,8.0}},{{180000.,9.0}},
967 {{190000.,6.0}},{{190000.,7.0}},{{190000.,8.0}},{{190000.,9.0}},
968 {{200000.,7.0}},{{200000.,8.0}},{{200000.,9.0}},
969 {{300000.,7.0}},{{300000.,8.0}},{{300000.,9.0}},
970 {{400000.,8.0}},{{400000.,9.0}},
971 {{500000.,8.0}},{{500000.,9.0}},
980 {{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}},
981 {{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}},
982 {{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}},
983 {{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}},
984 {{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}},
985 {{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}},
986 {{110000.,6.0}},{{110000.,6.5}},{{110000.,7.0}},{{110000.,7.5}},{{110000.,8.0}},{{110000.,8.5}},{{110000.,9.0}},
987 {{120000.,6.0}},{{120000.,6.5}},{{120000.,7.0}},{{120000.,7.5}},{{120000.,8.0}},{{120000.,8.5}},{{120000.,9.0}},
988 {{130000.,6.0}},{{130000.,6.5}},{{130000.,7.0}},{{130000.,7.5}},{{130000.,8.0}},{{130000.,8.5}},{{130000.,9.0}},
989 {{140000.,6.0}},{{140000.,6.5}},{{140000.,7.0}},{{140000.,7.5}},{{140000.,8.0}},{{140000.,8.5}},{{140000.,9.0}},
990 {{150000.,6.0}},{{150000.,6.5}},{{150000.,7.0}},{{150000.,7.5}},{{150000.,8.0}},{{150000.,8.5}},{{150000.,9.0}},
991 {{160000.,6.0}},{{160000.,6.5}},{{160000.,7.0}},{{160000.,7.5}},{{160000.,8.0}},{{160000.,8.5}},{{160000.,9.0}},
992 {{170000.,6.0}},{{170000.,6.5}},{{170000.,7.0}},{{170000.,7.5}},{{170000.,8.0}},{{170000.,8.5}},{{170000.,9.0}},
993 {{180000.,6.0}},{{180000.,6.5}},{{180000.,7.0}},{{180000.,7.5}},{{180000.,8.0}},{{180000.,8.5}},{{180000.,9.0}},
994 {{190000.,6.0}},{{190000.,6.5}},{{190000.,7.0}},{{190000.,7.5}},{{190000.,8.0}},{{190000.,8.5}},{{190000.,9.0}}
998 static const double par2[2] = { 0., -1. };
1001 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 };
1005 fprintf(
ioQQQ,
" RauchCompile on the job.\n" );
1013 fprintf(
ioQQQ,
" Creating rauch_h-ca_solar.ascii....\n" );
1020 fprintf(
ioQQQ,
" Creating rauch_h-ca_halo.ascii....\n" );
1029 fprintf(
ioQQQ,
" Creating rauch_h-ca_3d.ascii....\n" );
1037 if(
lgFileReadable(
"0050000_50_solar_iron.bin_0.1", dum, as ) &&
1040 fprintf(
ioQQQ,
" Creating rauch_h-ni_solar.ascii....\n" );
1041 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h-ni_solar.ascii",
"_solar_iron.bin_0.1",
1045 if(
lgFileReadable(
"0050000_50_halo__iron.bin_0.1", dum, as ) &&
1048 fprintf(
ioQQQ,
" Creating rauch_h-ni_halo.ascii....\n" );
1049 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h-ni_halo.ascii",
"_halo__iron.bin_0.1",
1053 if(
lgFileReadable(
"0050000_50_solar_iron.bin_0.1", dum, as ) &&
1057 fprintf(
ioQQQ,
" Creating rauch_h-ni_3d.ascii....\n" );
1065 if(
lgFileReadable(
"0040000_5.00_33_50_02_15.bin_0.1", dum, as ) &&
1068 fprintf(
ioQQQ,
" Creating rauch_pg1159.ascii....\n" );
1074 if(
lgFileReadable(
"0020000_4.00_H_00005-02000A.bin_0.1", dum, as ) &&
1077 fprintf(
ioQQQ,
" Creating rauch_hydr.ascii....\n" );
1078 lgFail = lgFail ||
RauchInitializeSub(
"rauch_hydr.ascii",
"_H_00005-02000A.bin_0.1",
1083 if(
lgFileReadable(
"0050000_5.00_He_00005-02000A.bin_0.1", dum, as ) &&
1086 fprintf(
ioQQQ,
" Creating rauch_helium.ascii....\n" );
1087 lgFail = lgFail ||
RauchInitializeSub(
"rauch_helium.ascii",
"_He_00005-02000A.bin_0.1",
1092 if(
lgFileReadable(
"0050000_5.00_H+He_1.000_0.000_00005-02000A.bin_0.1", dum, as ) &&
1095 fprintf(
ioQQQ,
" Creating rauch_h+he_3d.ascii....\n" );
1096 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_1.000_0.000_00005-02000A.bin_0.1",
1098 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.900_0.100_00005-02000A.bin_0.1",
1100 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.800_0.200_00005-02000A.bin_0.1",
1102 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.700_0.300_00005-02000A.bin_0.1",
1104 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.600_0.400_00005-02000A.bin_0.1",
1106 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.500_0.500_00005-02000A.bin_0.1",
1108 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.400_0.600_00005-02000A.bin_0.1",
1110 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.300_0.700_00005-02000A.bin_0.1",
1112 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.200_0.800_00005-02000A.bin_0.1",
1114 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.100_0.900_00005-02000A.bin_0.1",
1116 lgFail = lgFail ||
RauchInitializeSub(
"rauch_h+he_3d.ascii",
"_H+He_0.000_1.000_00005-02000A.bin_0.1",
1132 Edges[0] = 0.99946789f;
1133 Edges[1] = 1.8071406f;
1134 Edges[2] = 3.9996377f;
1137 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_solar.ascii",
"rauch_h-ca_solar.mod",Edges,3L, pc );
1139 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_halo.ascii",
"rauch_h-ca_halo.mod", Edges, 3L, pc );
1141 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_3d.ascii",
"rauch_h-ca_3d.mod", Edges, 3L, pc );
1144 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_solar.ascii",
"rauch_h-ni_solar.mod",Edges,3L, pc );
1146 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_halo.ascii",
"rauch_h-ni_halo.mod", Edges, 3L, pc );
1148 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_3d.ascii",
"rauch_h-ni_3d.mod", Edges, 3L, pc );
1151 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_pg1159.ascii",
"rauch_pg1159.mod", Edges, 3L, pc );
1154 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_hydr.ascii",
"rauch_hydr.mod", Edges, 3L, pc );
1157 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_helium.ascii",
"rauch_helium.mod", Edges, 3L, pc );
1160 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h+he_3d.ascii",
"rauch_h+he_3d.mod", Edges, 3L, pc );
1178 grid.
name =
"rauch_h-ca_3d.mod";
1180 grid.
name = ( lgHalo ?
"rauch_h-ca_halo.mod" :
"rauch_h-ca_solar.mod" );
1184 grid.
ident =
" H-Ca Rauch";
1186 grid.
command =
"COMPILE STARS";
1190 CheckVal( &grid, val, nval, ndim );
1212 grid.
name =
"rauch_h-ni_3d.mod";
1214 grid.
name = ( lgHalo ?
"rauch_h-ni_halo.mod" :
"rauch_h-ni_solar.mod" );
1218 grid.
ident =
" H-Ni Rauch";
1220 grid.
command =
"COMPILE STARS";
1224 CheckVal( &grid, val, nval, ndim );
1244 grid.
name =
"rauch_pg1159.mod";
1248 grid.
ident =
"PG1159 Rauch";
1250 grid.
command =
"COMPILE STARS";
1254 CheckVal( &grid, val, nval, ndim );
1274 grid.
name =
"rauch_hydr.mod";
1278 grid.
ident =
" Hydr Rauch";
1280 grid.
command =
"COMPILE STARS";
1284 CheckVal( &grid, val, nval, ndim );
1304 grid.
name =
"rauch_helium.mod";
1308 grid.
ident =
"Helium Rauch";
1310 grid.
command =
"COMPILE STARS";
1314 CheckVal( &grid, val, nval, ndim );
1334 grid.
name =
"rauch_h+he_3d.mod";
1338 grid.
ident =
" H+He Rauch";
1340 grid.
command =
"COMPILE STARS";
1344 CheckVal( &grid, val, nval, ndim );
1354 const char chOutName[])
1358 bool lgHeader =
true;
1359 long int i, j, nmods, ngp;
1361 size_t nsb_sz = (size_t)
NSB99;
1363 double *wavl, *fluxes[
MNTS], Age[
MNTS], lwavl;
1370 for( i=0; i <
MNTS; i++ )
1374 wavl = (
double *)
MALLOC(
sizeof(
double)*nsb_sz);
1386 double cage, cwavl, cfl1;
1390 if( sscanf( chLine,
" %le %le %le", &cage, &cwavl, &cfl1 ) != 3 )
1392 fprintf(
ioQQQ,
"syntax error in data of Starburst grid.\n" );
1403 fprintf(
ioQQQ,
"too many time steps in Starburst grid.\n" );
1404 fprintf(
ioQQQ,
"please increase MNTS and recompile.\n" );
1411 fluxes[nmods] = (
double *)
MALLOC(
sizeof(
double)*nsb_sz);
1415 if( ngp >= (
long)nsb_sz )
1421 fluxes[0] = (
double *)
REALLOC(fluxes[0],
sizeof(
double)*nsb_sz);
1422 wavl = (
double *)
REALLOC(wavl,
sizeof(
double)*nsb_sz);
1425 if( fabs(Age[nmods]-cage) > 10.*DBL_EPSILON*Age[nmods] )
1427 fprintf(
ioQQQ,
"age error in Starburst grid.\n" );
1435 if( fabs(wavl[ngp]-cwavl) > 10.*DBL_EPSILON*wavl[ngp] )
1437 fprintf(
ioQQQ,
"wavelength error in Starburst grid.\n" );
1444 fluxes[nmods][ngp] = pow( 10., cfl1 - 44.077911 );
1450 if( lgHeader && strncmp( &chLine[1],
"TIME [YR]", 9 ) == 0 )
1457 fprintf(
ioQQQ,
"syntax error in header of Starburst grid.\n" );
1469 fprintf( ioOut,
" %d\n", 1 );
1470 fprintf( ioOut,
" %d\n", 1 );
1471 fprintf( ioOut,
" Age\n" );
1472 fprintf( ioOut,
" %ld\n", nmods );
1473 fprintf( ioOut,
" %ld\n", ngp );
1475 fprintf( ioOut,
" lambda\n" );
1477 fprintf( ioOut,
" %.8e\n", 1. );
1479 fprintf( ioOut,
" F_lambda\n" );
1481 fprintf( ioOut,
" %.8e\n", 1. );
1483 for( i=0; i < nmods; i++ )
1485 fprintf( ioOut,
" %.3e", Age[i] );
1486 if( ((i+1)%4) == 0 )
1487 fprintf( ioOut,
"\n" );
1490 fprintf( ioOut,
"\n" );
1492 fprintf(
ioQQQ,
" Writing: " );
1495 for( j=0; j < ngp; j++ )
1497 fprintf( ioOut,
" %.4e", wavl[j] );
1499 if( ((j+1)%5) == 0 )
1500 fprintf( ioOut,
"\n" );
1504 fprintf( ioOut,
"\n" );
1507 fprintf(
ioQQQ,
"." );
1510 for( i=0; i < nmods; i++ )
1512 for( j=0; j < ngp; j++ )
1514 fprintf( ioOut,
" %.4e", fluxes[i][j] );
1516 if( ((j+1)%5) == 0 )
1517 fprintf( ioOut,
"\n" );
1521 fprintf( ioOut,
"\n" );
1524 fprintf(
ioQQQ,
"." );
1528 fprintf(
ioQQQ,
" done.\n" );
1533 for( i=0; i <
MNTS; i++ )
1539 for( i=0; i <
MNTS; i++ )
1550 bool lgFail =
false;
1554 fprintf(
ioQQQ,
" StarburstCompile on the job.\n" );
1562 lgFail = lgFail ||
lgCompileAtmosphere(
"starburst99.ascii",
"starburst99.mod", Edges, 0L, pc );
1572 bool lgFail =
false;
1576 fprintf(
ioQQQ,
" TlustyCompile on the job.\n" );
1581 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_p03.ascii",
"bstar2006_p03.mod", Edges, 0L, pc );
1583 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_p00.ascii",
"bstar2006_p00.mod", Edges, 0L, pc );
1585 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m03.ascii",
"bstar2006_m03.mod", Edges, 0L, pc );
1587 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m07.ascii",
"bstar2006_m07.mod", Edges, 0L, pc );
1589 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m10.ascii",
"bstar2006_m10.mod", Edges, 0L, pc );
1591 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m99.ascii",
"bstar2006_m99.mod", Edges, 0L, pc );
1594 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_3d.ascii",
"bstar2006_3d.mod", Edges, 0L, pc );
1597 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_p03.ascii",
"ostar2002_p03.mod", Edges, 0L, pc );
1599 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_p00.ascii",
"ostar2002_p00.mod", Edges, 0L, pc );
1601 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m03.ascii",
"ostar2002_m03.mod", Edges, 0L, pc );
1603 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m07.ascii",
"ostar2002_m07.mod", Edges, 0L, pc );
1605 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m10.ascii",
"ostar2002_m10.mod", Edges, 0L, pc );
1607 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m15.ascii",
"ostar2002_m15.mod", Edges, 0L, pc );
1609 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m17.ascii",
"ostar2002_m17.mod", Edges, 0L, pc );
1611 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m20.ascii",
"ostar2002_m20.mod", Edges, 0L, pc );
1613 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m30.ascii",
"ostar2002_m30.mod", Edges, 0L, pc );
1615 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m99.ascii",
"ostar2002_m99.mod", Edges, 0L, pc );
1618 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_3d.ascii",
"ostar2002_3d.mod", Edges, 0L, pc );
1627 const char *chMetalicity,
1638 grid.
name =
"bstar2006_";
1640 grid.
name =
"ostar2002_";
1646 grid.
name += chMetalicity;
1647 grid.
name +=
".mod";
1653 strcpy( chIdent,
"3-dim" );
1657 strcpy( chIdent,
"Z " );
1658 strcat( chIdent, chMetalicity );
1661 strcat( chIdent,
" Bstr06" );
1663 strcat( chIdent,
" Ostr02" );
1666 grid.
ident = chIdent;
1668 grid.
command =
"COMPILE STARS";
1672 CheckVal( &grid, val, nval, ndim );
1687 bool lgFail =
false;
1691 fprintf(
ioQQQ,
" WernerCompile on the job.\n" );
1705 Edges[0] = 0.99946789f;
1706 Edges[1] = 1.8071406f;
1707 Edges[2] = 3.9996377f;
1754 grid.
name =
"kwerner.mod";
1758 grid.
ident =
"Klaus Werner";
1760 grid.
command =
"COMPILE STARS";
1764 CheckVal( &grid, val, nval, ndim );
1798 bool lgFail =
false;
1802 fprintf(
ioQQQ,
" WMBASICCompile on the job.\n" );
1805 Edges[0] = 0.99946789f;
1806 Edges[1] = 1.8071406f;
1807 Edges[2] = 3.9996377f;
1828 grid.
name =
"wmbasic.mod";
1832 grid.
ident =
" WMBASIC";
1834 grid.
command =
"COMPILE STARS";
1838 CheckVal( &grid, val, nval, ndim );
1849 const char chFNameOut[],
1858 long int i, j, nskip, nModels, nWL;
1861 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL;
1888 fprintf(
ioQQQ,
" CompileAtmosphereCoStar got %s.\n", chFNameIn );
1893 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading nskip.\n" );
1896 sscanf( chLine,
"%li", &nskip );
1899 for( i=0; i < nskip; ++i )
1903 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails skipping header.\n" );
1911 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading nModels, nWL.\n" );
1914 sscanf( chLine,
"%li%li", &nModels, &nWL );
1916 if( nModels <= 0 || nWL <= 0 )
1918 fprintf(
ioQQQ,
" CompileAtmosphereCoStar scanned off impossible values for nModels=%li or nWL=%li\n",
1927 for( i=0; i < nModels; ++i )
1931 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading model parameters.\n" );
1935 telg[i].
chGrid = chLine[0];
1937 sscanf( chLine+1,
"%i", &telg[i].modid );
1939 sscanf( chLine+23,
"%lg", &telg[i].par[0] );
1941 sscanf( chLine+31,
"%lg", &telg[i].par[1] );
1943 sscanf( chLine+7,
"%lg", &telg[i].par[2] );
1945 sscanf( chLine+15,
"%lg", &telg[i].par[3] );
1948 ASSERT( telg[i].par[2] > 10. );
1949 ASSERT( telg[i].par[3] > 10. );
1952 telg[i].
par[2] = log10(telg[i].par[2]);
1970 val[1] = (int32)
MDIM;
1971 val[2] = (int32)
MNAM;
1974 val[5] = (int32)nModels;
1976 uval[0] =
sizeof(val) +
sizeof(uval) +
sizeof(names) + nModels*
sizeof(
mpp);
1979 strncpy( names[0],
"Teff\0\0",
MNAM+1 );
1980 strncpy( names[1],
"log(g)",
MNAM+1 );
1981 strncpy( names[2],
"log(M)",
MNAM+1 );
1982 strncpy( names[3],
"Age\0\0\0",
MNAM+1 );
1984 if( fwrite( val,
sizeof(val), 1, ioOUT ) != 1 ||
1985 fwrite( uval,
sizeof(uval), 1, ioOUT ) != 1 ||
1986 fwrite( names,
sizeof(names), 1, ioOUT ) != 1 ||
1988 fwrite( telg,
sizeof(
mpp), (
size_t)nModels, ioOUT ) != (
size_t)nModels ||
1990 fwrite(
rfield.
AnuOrg, (
size_t)uval[1], 1, ioOUT ) != 1 )
1992 fprintf(
ioQQQ,
" CompileAtmosphereCoStar failed writing header of output file.\n" );
2001 fprintf(
ioQQQ,
" Compiling: " );
2004 for( i=0; i < nModels; ++i )
2009 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading the skip to next spectrum.\n" );
2012 sscanf( chLine,
"%li", &nskip );
2014 for( j=0; j < nskip; ++j )
2018 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails doing the skip.\n" );
2026 for( j=nWL-1; j >= 0; --j )
2030 fprintf(
ioQQQ,
" CompileAtmosphereCoStar fails reading the spectral data.\n" );
2033 double help1, help2;
2034 sscanf( chLine,
"%lg %lg", &help1, &help2 );
2038 StarFlux[j] = (
realnum)(
PI*pow(10.,help2));
2044 ASSERT( StarEner[j] < StarEner[j+1] );
2050 RebinAtmosphere(nWL-1, StarEner+1, StarFlux+1, CloudyFlux, nedges, Edges );
2053 if( fwrite( CloudyFlux, (
size_t)uval[1], 1, ioOUT ) != 1 )
2055 fprintf(
ioQQQ,
" CompileAtmosphereCoStar failed writing star flux.\n" );
2059 fprintf(
ioQQQ,
"." );
2063 fprintf(
ioQQQ,
" done.\n" );
2073 fprintf(
ioQQQ,
"\n CompileAtmosphereCoStar completed ok.\n\n" );
2095 long i, j, k, nmodid, off, ptr;
2096 long *indloTr, *indhiTr, useTr[2];
2097 long indlo[2], indhi[2], index[2];
2099 double lval[2], aval[4];
2103 switch( grid->
imode )
2112 lval[0] = log10(val[0]);
2118 lval[0] = log10(val[1]);
2123 fprintf(
ioQQQ,
" InterpolateGridCoStar called with insane value for imode: %d.\n", grid->
imode );
2127 nmodid = (long)(lval[1]+0.5);
2145 for( j=0; j < grid->
nTracks; j++ )
2149 if( grid->
trackLen[j] >= nmodid ) {
2150 index[0] = nmodid - 1;
2161 ValTr[j] = -FLT_MAX;
2166 FindHCoStar( grid, j, lval[1], off, ValTr, indloTr, indhiTr );
2171 for( j=0; j < grid->
nTracks; j++ )
2173 if( indloTr[j] >= 0 )
2174 printf(
"track %c: models %c%d, %c%d, val %g\n",
2189 fprintf(
ioQQQ,
" The parameters for the requested CoStar model are out of range.\n" );
2195 ASSERT( indloTr[useTr[0]] >= 0 && indloTr[useTr[0]] < (
int)grid->
nmods );
2196 ASSERT( indhiTr[useTr[0]] >= 0 && indhiTr[useTr[0]] < (
int)grid->
nmods );
2197 ASSERT( indloTr[useTr[1]] >= 0 && indloTr[useTr[1]] < (
int)grid->
nmods );
2198 ASSERT( indhiTr[useTr[1]] >= 0 && indhiTr[useTr[1]] < (
int)grid->
nmods );
2201 printf(
"interpolate between tracks %c and %c\n", (
char)(
'A'+useTr[0]), (
char)(
'A'+useTr[1]) );
2204 indlo[0] = indloTr[useTr[0]];
2205 indhi[0] = indhiTr[useTr[0]];
2206 indlo[1] = indloTr[useTr[1]];
2207 indhi[1] = indhiTr[useTr[1]];
2220 FILE *ioBUG = fopen(
"interpolated.txt",
"w" );
2231 SetLimits( grid, lval[0], NULL, NULL, useTr, ValTr, val0_lo, val0_hi );
2236 fprintf(
ioQQQ,
" * c<< FINAL: T_eff = %7.1f, ", aval[0] );
2237 fprintf(
ioQQQ,
"log(g) = %4.2f, M(ZAMS) = %5.1f, age = ", aval[1], pow(10.,aval[2]) );
2239 fprintf(
ioQQQ,
" >>> *\n" );
2257 long index[2], j, mod1, mod2;
2261 indloTr[track] = -2;
2262 indhiTr[track] = -2;
2263 ValTr[track] = -FLT_MAX;
2267 for( j=0; j < grid->
trackLen[track]; j++ )
2273 if( fabs(par2-grid->
telg[mod1].
par[off+1]) <= 10.*FLT_EPSILON*fabs(grid->
telg[mod1].
par[off+1]) )
2275 indloTr[track] = mod1;
2276 indhiTr[track] = mod1;
2282 for( j=0; j < grid->
trackLen[track]-1; j++ )
2290 if( (par2 - grid->
telg[mod1].
par[off+1])*(par2 - grid->
telg[mod2].
par[off+1]) < 0. )
2294 indloTr[track] = mod1;
2295 indhiTr[track] = mod2;
2296 frac = (par2 - grid->
telg[mod2].
par[off+1])/
2299 (1.-frac)*grid->
telg[mod2].
par[off] );
2323 for( j=0; j < grid->
nTracks; j++ )
2326 if( ValTr[j] != -FLT_MAX && fabs(par1-(
double)ValTr[j]) <= 10.*FLT_EPSILON*fabs(ValTr[j]) )
2339 for( j=0; j < grid->
nTracks-1; j++ )
2341 if( ValTr[j] != -FLT_MAX )
2347 for( i = j+1; i < grid->
nTracks; i++ )
2349 if( ValTr[i] != -FLT_MAX )
2357 if( j2 > 0 && ((
realnum)par1-ValTr[j])*((
realnum)par1-ValTr[j2]) < 0.f )
2374 long index[2], maxlen, n;
2379 for( n=0; n < grid->
nTracks; n++ )
2382 fprintf(
ioQQQ,
"\n" );
2383 fprintf(
ioQQQ,
" Track\\Index |" );
2384 for( n = 0; n < maxlen; n++ )
2385 fprintf(
ioQQQ,
" %5ld ", n+1 );
2386 fprintf(
ioQQQ,
"\n" );
2387 fprintf(
ioQQQ,
"--------------|" );
2388 for( n = 0; n < maxlen; n++ )
2389 fprintf(
ioQQQ,
"----------------" );
2390 fprintf(
ioQQQ,
"\n" );
2392 for( index[1]=0; index[1] < grid->
nTracks; ++index[1] )
2395 double Teff, alogg, Mass;
2397 fprintf(
ioQQQ,
" %c", (
char)(
'A'+index[1]) );
2400 Mass = pow(10.,grid->
telg[ptr].
par[2]);
2401 fprintf(
ioQQQ,
" (%3.0f Msol) |", Mass );
2403 for( index[0]=0; index[0] < grid->
trackLen[index[1]]; ++index[0] )
2406 Teff = grid->
telg[ptr].
par[0];
2407 alogg = grid->
telg[ptr].
par[1];
2408 fprintf(
ioQQQ,
" (%6.1f,%4.2f)", Teff, alogg );
2410 fprintf(
ioQQQ,
"\n" );
2417 const char chSuff[],
2422 const double par2[],
2432 double *wavl, *fluxes;
2455 fprintf( ioOut,
" %d\n", ( ngrids == 1 ? 2 : 3 ) );
2456 fprintf( ioOut,
" %d\n", ( ngrids == 1 ? 2 : 3 ) );
2457 fprintf( ioOut,
" Teff\n" );
2458 fprintf( ioOut,
" log(g)\n" );
2460 fprintf( ioOut,
" log(Z)\n" );
2461 else if( ngrids == 11 )
2462 fprintf( ioOut,
" f(He)\n" );
2464 fprintf( ioOut,
" %ld\n", nmods*ngrids );
2465 fprintf( ioOut,
" %d\n",
NRAUCH );
2467 fprintf( ioOut,
" lambda\n" );
2469 fprintf( ioOut,
" %.8e\n", 1. );
2471 fprintf( ioOut,
" F_lambda\n" );
2473 fprintf( ioOut,
" %.8e\n",
PI*1.e-8 );
2475 for( j=0; j < ngrids; j++ )
2478 for( i=0; i < nmods; i++ )
2481 fprintf( ioOut,
" %.0f %.1f", telg[i].par[0], telg[i].par[1] );
2483 fprintf( ioOut,
" %.0f %.1f %.1f", telg[i].par[0], telg[i].par[1], par2[j] );
2484 if( ((i+1)%4) == 0 )
2485 fprintf( ioOut,
"\n" );
2488 fprintf( ioOut,
"\n" );
2491 fprintf(
ioQQQ,
" Writing: " );
2494 for( i=0; i < nmods; i++ )
2498 sprintf( chLine,
"%7.7ld_%2ld", (
long)(telg[i].par[0]+0.5), (
long)(10.*telg[i].par[1]+0.5) );
2499 else if( format == 2 )
2500 sprintf( chLine,
"%7.7ld_%.2f", (
long)(telg[i].par[0]+0.5), telg[i].par[1] );
2503 fprintf(
ioQQQ,
" insanity in RauchInitializeSub\n" );
2507 string chFileName( chLine );
2508 chFileName += chSuff;
2523 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2529 while( chLine[0] ==
'*' )
2533 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2540 for( j=0; j <
NRAUCH; j++ )
2549 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2556 if( sscanf( chLine,
"%lf %le", &wl, &ttemp ) != 2 )
2558 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2568 if( fabs(wavl[j]-wl) > 10.*DBL_EPSILON*wl )
2570 fprintf(
ioQQQ,
" RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2582 if( i == 0 && n == 1 )
2585 for( j=0; j <
NRAUCH; j++ )
2587 fprintf( ioOut,
" %.4e", wavl[j] );
2589 if( ((j+1)%5) == 0 )
2590 fprintf( ioOut,
"\n" );
2594 fprintf( ioOut,
"\n" );
2597 for( j=0; j <
NRAUCH; j++ )
2599 fprintf( ioOut,
" %.4e", fluxes[j] );
2601 if( ((j+1)%5) == 0 )
2602 fprintf( ioOut,
"\n" );
2606 fprintf( ioOut,
"\n" );
2609 fprintf(
ioQQQ,
"." );
2614 fprintf(
ioQQQ,
" done.\n" );
2632 const char chFNameOut[],
2640 char chDataType[11];
2643 bool lgFreqX, lgFreqY, lgFlip;
2646 long int i, imod, version, nd, ndim, npar, nmods, ngrid;
2649 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL, *scratch = NULL;
2651 double convert_wavl, convert_flux;
2665 fprintf(
ioQQQ,
" lgCompileAtmosphere got %s.\n", chFNameIn );
2668 if( fscanf( ioIN,
"%ld", &version ) != 1 )
2670 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading VERSION.\n" );
2676 fprintf(
ioQQQ,
" lgCompileAtmosphere: there is a version number mismatch in"
2677 " the ascii atmosphere file: %s.\n", chFNameIn );
2678 fprintf(
ioQQQ,
" lgCompileAtmosphere: Please recreate this file or download the"
2679 " latest version following the instructions on the Cloudy website.\n" );
2684 if( fscanf( ioIN,
"%ld", &ndim ) != 1 || ndim <= 0 || ndim >
MDIM )
2686 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid dimension of grid.\n" );
2691 if( fscanf( ioIN,
"%ld", &npar ) != 1 || npar <= 0 || npar < ndim || npar > MDIM )
2693 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid no. of model parameters.\n" );
2698 memset( names,
'\0', MDIM*(
MNAM+1) );
2700 for( nd=0; nd < npar; nd++ )
2702 if( fscanf( ioIN,
"%6s", names[nd] ) != 1 )
2704 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading parameter label.\n" );
2708 if( strcmp( names[0],
"Teff" ) != 0 && strcmp( names[0],
"Age" ) != 0 )
2710 fprintf(
ioQQQ,
" lgCompileAtmosphere: first parameter must be \"Teff\" or \"Age\".\n" );
2713 if( ndim > 1 && strcmp( names[0],
"Age" ) == 0 )
2715 fprintf(
ioQQQ,
" First parameter is \"Age\", but ndim is not 1.\n" );
2718 if( ndim >= 2 && strcmp( names[1],
"log(g)" ) != 0 )
2720 fprintf(
ioQQQ,
" lgCompileAtmosphere: second parameter must be \"log(g)\".\n" );
2725 if( fscanf( ioIN,
"%ld", &nmods ) != 1 || nmods <= 0 )
2727 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid number of models.\n" );
2731 if( fscanf( ioIN,
"%ld", &ngrid ) != 1 || ngrid <= 1 )
2733 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid number of grid points.\n" );
2738 if( fscanf( ioIN,
"%10s", chDataType ) != 1 )
2740 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading wavl DataType string.\n" );
2744 if( strcmp( chDataType,
"lambda" ) == 0 )
2746 else if( strcmp( chDataType,
"nu" ) == 0 )
2749 fprintf(
ioQQQ,
" lgCompileAtmosphere found illegal wavl DataType: %s.\n", chDataType );
2753 if( fscanf( ioIN,
"%le", &convert_wavl ) != 1 || convert_wavl <= 0. )
2755 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid wavl conversion factor.\n" );
2760 if( fscanf( ioIN,
"%10s", chDataType ) != 1 )
2762 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading flux DataType string.\n" );
2766 if( strcmp( chDataType,
"F_lambda" ) == 0 || strcmp( chDataType,
"H_lambda" ) == 0 )
2768 else if( strcmp( chDataType,
"F_nu" ) == 0 || strcmp( chDataType,
"H_nu" ) == 0 )
2771 fprintf(
ioQQQ,
" lgCompileAtmosphere found illegal flux DataType: %s.\n", chDataType );
2775 if( fscanf( ioIN,
"%le", &convert_flux ) != 1 || convert_flux <= 0. )
2777 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid flux conversion factor.\n" );
2783 for( i=0; i < nmods; i++ )
2785 for( nd=0; nd < npar; nd++ )
2787 if( fscanf( ioIN,
"%le", &telg[i].par[nd] ) != 1 )
2789 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid model parameter.\n" );
2793 if( telg[i].par[0] <= 0. )
2795 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading valid %s.\n", names[0] );
2813 val[1] = (int32)MDIM;
2814 val[2] = (int32)
MNAM;
2815 val[3] = (int32)ndim;
2816 val[4] = (int32)npar;
2817 val[5] = (int32)nmods;
2819 uval[0] =
sizeof(val) +
sizeof(uval) +
sizeof(names) + nmods*
sizeof(
mpp);
2822 if( fwrite( val,
sizeof(val), 1, ioOUT ) != 1 ||
2823 fwrite( uval,
sizeof(uval), 1, ioOUT ) != 1 ||
2824 fwrite( names,
sizeof(names), 1, ioOUT ) != 1 ||
2826 fwrite( telg,
sizeof(
mpp), (size_t)nmods, ioOUT ) != (size_t)nmods ||
2828 fwrite(
rfield.
AnuOrg, (
size_t)uval[1], 1, ioOUT ) != 1 )
2830 fprintf(
ioQQQ,
" lgCompileAtmosphere failed writing header of output file.\n" );
2841 for( i=0; i < ngrid; i++ )
2844 if( fscanf( ioIN,
"%lg", &help ) != 1 )
2846 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading wavelength.\n" );
2851 scratch[i] = (
realnum)(help*convert_wavl);
2853 ASSERT( scratch[i] > 0.f );
2856 lgFlip = ( !lgFreqX && scratch[0] < scratch[1] ) || ( lgFreqX && scratch[0] > scratch[1] );
2859 for( i=0; i < ngrid; i++ )
2868 StarEner[ngrid-i-1] = scratch[i];
2870 StarEner[i] = scratch[i];
2873 ASSERT( StarEner[0] > 0.f );
2875 for( i=1; i < ngrid; i++ )
2876 ASSERT( StarEner[i] > StarEner[i-1] );
2878 fprintf(
ioQQQ,
" Compiling: " );
2880 for( imod=0; imod < nmods; imod++ )
2885 for( i=0; i < ngrid; i++ )
2888 if( fscanf( ioIN,
"%lg", &help ) != 1 )
2890 fprintf(
ioQQQ,
" lgCompileAtmosphere failed reading star flux.\n" );
2895 scratch[i] = (
realnum)(help*convert_flux);
2898 ASSERT( scratch[i] >= 0.f );
2901 for( i=0; i < ngrid; i++ )
2904 StarFlux[ngrid-i-1] = scratch[i];
2906 StarFlux[i] = scratch[i];
2909 for( i=0; i < ngrid; i++ )
2913 StarFlux[i] *= CONVERT_FNU/
POW2(StarEner[i]);
2914 ASSERT( StarFlux[i] >= 0.f );
2918 RebinAtmosphere( ngrid, StarEner, StarFlux, CloudyFlux, nedges, Edges );
2921 if( fwrite( CloudyFlux, (
size_t)uval[1], 1, ioOUT ) != 1 )
2923 fprintf(
ioQQQ,
" lgCompileAtmosphere failed writing star flux.\n" );
2927 fprintf(
ioQQQ,
"." );
2931 fprintf(
ioQQQ,
" done.\n" );
2943 fprintf(
ioQQQ,
" lgCompileAtmosphere completed ok.\n\n" );
2961 int32 version, mdim, mnam;
2973 fprintf(
ioQQQ,
" Error: stellar atmosphere file not found.\n" );
2974 fprintf(
ioQQQ ,
"\n\n If the path is set then it is possible that the stellar atmosphere data files do not exist.\n");
2975 fprintf(
ioQQQ ,
" Have the stellar data files been downloaded and compiled with the COMPILE STARS command?\n");
2976 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");
2981 if( fread( &version,
sizeof(version), 1, grid->
ioIN ) != 1 ||
2982 fread( &mdim,
sizeof(mdim), 1, grid->
ioIN ) != 1 ||
2983 fread( &mnam,
sizeof(mnam), 1, grid->
ioIN ) != 1 ||
2984 fread( &grid->
ndim,
sizeof(grid->
ndim), 1, grid->
ioIN ) != 1 ||
2985 fread( &grid->
npar,
sizeof(grid->
npar), 1, grid->
ioIN ) != 1 ||
2986 fread( &grid->
nmods,
sizeof(grid->
nmods), 1, grid->
ioIN ) != 1 ||
2987 fread( &grid->
ngrid,
sizeof(grid->
ngrid), 1, grid->
ioIN ) != 1 ||
2991 fprintf(
ioQQQ,
" InitGrid failed reading header.\n" );
2998 fprintf(
ioQQQ,
" InitGrid: there is a version mismatch between"
2999 " the compiled atmospheres file I expected and the one I found.\n" );
3000 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3001 " atmospheres file with the command: %s.\n", grid->
command );
3007 fprintf(
ioQQQ,
" InitGrid: the compiled atmospheres file is produced"
3008 " with an incompatible version of Cloudy.\n" );
3009 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3010 " atmospheres file with the command: %s.\n", grid->
command );
3023 if( fread( &grid->
names,
sizeof(grid->
names), 1, grid->
ioIN ) != 1 )
3025 fprintf(
ioQQQ,
" InitGrid failed reading names array.\n" );
3030 grid->
val = (
double **)
MALLOC(
sizeof(
double*)*grid->
ndim );
3031 for( nd=0; nd < grid->
ndim; nd++ )
3039 fprintf(
ioQQQ,
" InitGrid failed reading model parameter block.\n" );
3047 int res = fseek( grid->
ioIN, 0, SEEK_END );
3050 long End = ftell( grid->
ioIN );
3052 if( End != Expected )
3054 fprintf(
ioQQQ,
" InitGrid: Problem performing sanity check for size of binary file.\n" );
3055 fprintf(
ioQQQ,
" InitGrid: I expected to find %ld bytes, but actually found %ld bytes.\n",
3057 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3058 " atmospheres file with the command: %s.\n", grid->
command );
3078 int32 version, mdim, mnam;
3083 grid.
name = binName;
3088 if( fread( &version,
sizeof(version), 1, grid.
ioIN ) != 1 ||
3089 fread( &mdim,
sizeof(mdim), 1, grid.
ioIN ) != 1 ||
3090 fread( &mnam,
sizeof(mnam), 1, grid.
ioIN ) != 1 ||
3091 fread( &grid.
ndim,
sizeof(grid.
ndim), 1, grid.
ioIN ) != 1 ||
3092 fread( &grid.
npar,
sizeof(grid.
npar), 1, grid.
ioIN ) != 1 ||
3098 fclose( grid.
ioIN );
3105 fclose( grid.
ioIN );
3113 int res = fseek( grid.
ioIN, 0, SEEK_END );
3116 long End = ftell( grid.
ioIN );
3118 if( End != Expected )
3120 fclose( grid.
ioIN );
3126 fclose( grid.
ioIN );
3140 if( (ioIN =
open_data( ascName,
"r", scheme )) == NULL )
3144 if( fscanf( ioIN,
"%ld", &version ) != 1 || version !=
VERSION_ASCII )
3168 grid->
jlo = grid->
jhi = NULL;
3171 memset( grid->
jval, 0xff, (
size_t)(grid->
nval[0]*grid->
nval[1]*
sizeof(
long)) );
3179 track = (char)(
'A'+index[1]);
3183 for( i=0; i < grid->
nmods; i++ )
3199 grid->
trackLen[index[1]] = index[0];
3215 *ndim = (long)grid->
ndim;
3216 if( *ndim == 2 && *nval == 1 )
3219 val[*nval] = grid->
val[1][grid->
nval[1]-1];
3222 if( *ndim != (
long)grid->
ndim )
3224 fprintf(
ioQQQ,
" A %ld-dim grid was requested, but a %ld-dim grid was found.\n",
3225 *ndim, (
long)grid->
ndim );
3230 fprintf(
ioQQQ,
" A %ld-dim grid was requested, but only %ld parameters were entered.\n",
3253 indlo = (
long *)
MALLOC((
size_t)(grid->
ndim*
sizeof(long)) );
3254 indhi = (
long *)
MALLOC((
size_t)(grid->
ndim*
sizeof(long)) );
3255 index = (
long *)
MALLOC((
size_t)(grid->
ndim*
sizeof(long)) );
3256 aval = (
double *)
MALLOC((
size_t)(grid->
npar*
sizeof(double)) );
3270 for( nd=0; nd < grid->
ndim; nd++ )
3272 FindIndex( grid->
val[nd], grid->
nval[nd], val[nd], &indlo[nd], &indhi[nd], &lgInvalid );
3276 " Requested parameter %s = %.2f is not within the range %.2f to %.2f\n",
3277 grid->
names[nd], val[nd], grid->
val[nd][0], grid->
val[nd][grid->
nval[nd]-1] );
3287 if( grid->
npar == 1 )
3289 " * c<< FINAL: %6s = %13.2f"
3291 grid->
names[0], aval[0] );
3292 else if( grid->
npar == 2 )
3294 " * c<< FINAL: %6s = %10.2f"
3295 " %6s = %8.5f >>> *\n",
3296 grid->
names[0], aval[0], grid->
names[1], aval[1] );
3297 else if( grid->
npar == 3 )
3299 " * c<< FINAL: %6s = %7.0f"
3300 " %6s = %5.2f %6s = %5.2f >>> *\n",
3301 grid->
names[0], aval[0], grid->
names[1], aval[1],
3302 grid->
names[2], aval[2] );
3303 else if( grid->
npar >= 4 )
3306 " * c<< FINAL: %4s = %7.0f"
3307 " %6s = %4.2f %6s = %5.2f %6s = ",
3308 grid->
names[0], aval[0], grid->
names[1], aval[1],
3311 fprintf(
ioQQQ,
" >>> *\n" );
3324 FILE *ioBUG = fopen(
"interpolated.txt",
"w" );
3330 if( strcmp( grid->
names[0],
"Teff" ) == 0 )
3337 SetLimits( grid, val[0], indlo, indhi, NULL, NULL, Tlow, Thigh );
3354 fclose( grid->
ioIN );
3356 for( i = 0; i < grid->
ndim; i++ )
3384 lgDryRun = ( flux1 == NULL );
3388 long n =
JIndex(grid,index);
3390 ind = ( grid->
jlo[n] >= 0 ) ? grid->
jlo[n] : grid->
jhi[n];
3392 ind = ( grid->
jhi[n] >= 0 ) ? grid->
jhi[n] : grid->
jlo[n];
3393 else if( grid->
ndim == 1 )
3401 fprintf(
ioQQQ,
" The requested interpolation could not be completed, sorry.\n" );
3402 fprintf(
ioQQQ,
" No suitable match was found for a model with" );
3403 for( i=0; i < grid->
ndim; i++ )
3404 fprintf(
ioQQQ,
" %s=%.6g ", grid->
names[i], grid->
val[i][index[i]] );
3405 fprintf(
ioQQQ,
"\n" );
3409 for( i=0; i < grid->
npar; i++ )
3410 aval[i] = grid->
telg[ind].
par[i];
3416 if( fabs(grid->
val[i][index[i]]-aval[i]) > 10.*DBL_EPSILON*fabs(aval[i]) )
3418 fprintf(
ioQQQ,
" No exact match was found for a model with" );
3419 for( j=0; j < grid->
ndim; j++ )
3420 fprintf(
ioQQQ,
" %s=%.6g ", grid->
names[j], grid->
val[j][index[j]] );
3421 fprintf(
ioQQQ,
"- using the following model instead:\n" );
3435 const realnum SECURE = 10.f*FLT_EPSILON;
3439 aval2 = (
double*)
MALLOC((
size_t)(grid->
npar*
sizeof(double)) );
3450 index[nd] = indlo[nd];
3451 InterpolateModel( grid, val, aval, indlo, indhi, index, nd, flux1, stage );
3456 index[nd] = indhi[nd];
3457 InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, NULL, stage );
3460 if( fabs(aval2[nd]-aval[nd]) > 10.*DBL_EPSILON*fabs(aval[nd]) )
3462 double fr1, fr2, fc1 = 0., fc2 = 0.;
3465 InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, flux2, stage );
3467 fr1 = (aval2[nd]-val[nd])/(aval2[nd]-aval[nd]);
3475 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3480 fprintf(
ioQQQ,
"interpolation nd=%ld fr1=%g\n", nd, fr1 );
3484 if( nd == 0 && strcmp( grid->
names[nd],
"Teff" ) == 0 )
3496 fc1 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->
val[nd][indlo[nd]])*4. : 0.;
3497 fc2 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->
val[nd][indhi[nd]])*4. : 0.;
3501 flux1[i] = (
realnum)(fr1*(flux1[i]+fc1) + fr2*(flux2[i]+fc2));
3504 for( i=0; i < grid->
npar; i++ )
3505 aval[i] = fr1*aval[i] + fr2*aval2[i];
3530 ind = ( index[1] == 0 ) ? indlo[index[0]] : indhi[index[0]];
3534 for( i=0; i < grid->
npar; i++ )
3535 aval[i] = grid->
telg[ind].
par[i];
3541 const realnum SECURE = 10.f*FLT_EPSILON;
3552 lgSkip = ( nd == 1 ) ? ( indhi[index[0]] == indlo[index[0]] ) :
3553 ( indlo[0] == indlo[1] && indhi[0] == indhi[1] );
3558 double fr1, fr2, *aval2;
3561 aval2 = (
double*)
MALLOC((
size_t)(grid->
npar*
sizeof(double)) );
3566 fr1 = (aval2[nd+off]-val[nd])/(aval2[nd+off]-aval[nd+off]);
3570 fprintf(
ioQQQ,
"interpolation nd=%ld fr1=%g\n", nd, fr1 );
3573 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3576 flux1[i] = (
realnum)(fr1*flux1[i] + fr2*flux2[i]);
3578 for( i=0; i < grid->
npar; i++ )
3579 aval[i] = fr1*aval[i] + fr2*aval2[i];
3604 ASSERT( ind >= 0 && ind <= grid->nmods );
3610 fprintf(
ioQQQ,
" Error seeking atmosphere %ld\n", ind );
3616 fprintf(
ioQQQ,
" Error trying to read atmosphere %ld\n", ind );
3624 if( grid->
npar == 1 )
3626 " * c<< %s model%5ld read. "
3627 " %6s = %13.2f >>> *\n",
3629 else if( grid->
npar == 2 )
3631 " * c<< %s model%5ld read. "
3632 " %6s = %10.2f %6s = %8.5f >>> *\n",
3635 else if( grid->
npar == 3 )
3637 " * c<< %s model%5ld read. "
3638 " %6s=%7.0f %6s=%5.2f %6s=%5.2f >>> *\n",
3642 else if( grid->
npar >= 4 )
3646 " %4s=%5.0f %6s=%4.2f %6s=%5.2f %6s=",
3651 fprintf(
ioQQQ,
" >> *\n" );
3659 flux[i] = (
realnum)log10(
MAX2( 1e-37, (
double)flux[i] ) );
3678 long index[
MDIM], j;
3679 const double SECURE = (1. + 20.*(double)FLT_EPSILON);
3684 switch( grid->
imode )
3694 for( j=0; j < grid->
nTracks; j++ )
3696 if( ValTr[j] != -FLT_MAX )
3700 pow(10.,(
double)ValTr[j]) : ValTr[j];
3701 *loLim =
MIN2(*loLim,temp);
3702 *hiLim =
MAX2(*hiLim,temp);
3708 index[1] = useTr[0];
3710 index[1] = useTr[1];
3714 printf(
"set limit 0: (models %d, %d) %f %f\n",
3715 ptr0+1, ptr1+1, grid->
telg[ptr0].
par[3], grid->
telg[ptr1].
par[3] );
3717 index[0] = grid->
trackLen[useTr[0]]-1;
3718 index[1] = useTr[0];
3720 index[0] = grid->
trackLen[useTr[1]]-1;
3721 index[1] = useTr[1];
3725 printf(
"set limit 1: (models %d, %d) %f %f\n",
3726 ptr0+1, ptr1+1, grid->
telg[ptr0].
par[3], grid->
telg[ptr1].
par[3] );
3730 fprintf(
ioQQQ,
" SetLimits called with insane value for imode: %d.\n", grid->
imode );
3734 ASSERT( fabs(*loLim) < DBL_MAX && fabs(*hiLim) < DBL_MAX );
3737 if( *hiLim <= *loLim )
3739 fprintf(
ioQQQ,
" no room to optimize: lower limit %.4f, upper limit %.4f.\n",
3749 printf(
"set limits: %g %g\n",*loLim,*hiLim);
3777 double loLoc = +DBL_MAX;
3778 double hiLoc = -DBL_MAX;
3781 for( index[0]=0; index[0] < grid->
nval[0]; ++index[0] )
3790 if( grid->
jhi[n] < 0 )
3794 if( grid->
val[0][index[0]] < val )
3796 if( grid->
val[0][index[0]] > val )
3803 if( grid->
val[0][index[0]] <= val )
3807 if( loLoc == DBL_MAX )
3808 loLoc = grid->
val[0][index[0]];
3810 if( grid->
val[0][index[0]] >= val )
3814 hiLoc = grid->
val[0][index[0]];
3819 ASSERT( fabs(loLoc) < DBL_MAX && fabs(hiLoc) < DBL_MAX && loLoc < hiLoc );
3821 *loLim =
MAX2(*loLim,loLoc);
3822 *hiLim =
MIN2(*hiLim,hiLoc);
3826 index[nd] = indlo[nd];
3827 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
3829 if( indhi[nd] != indlo[nd] )
3831 index[nd] = indhi[nd];
3832 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
3841 long i, *index, j, jsize, nd;
3852 for( nd=0; nd < grid->
ndim; nd++ )
3854 double pval = grid->
telg[0].
par[nd];
3855 grid->
val[nd][0] = pval;
3858 for( i=1; i < grid->
nmods; i++ )
3873 for( j = grid->
nval[nd]-1; j >= i2; j-- )
3874 grid->
val[nd][j+1] = grid->
val[nd][j];
3875 grid->
val[nd][i2] = pval;
3880 jsize *= grid->
nval[nd];
3883 printf(
"%s[%ld]:", grid->
names[nd], grid->
nval[nd] );
3884 for( i=0; i < grid->
nval[nd]; i++ )
3885 printf(
" %g", grid->
val[nd][i] );
3890 index = (
long *)
MALLOC(
sizeof(
long)*grid->
ndim );
3891 val = (
double *)
MALLOC(
sizeof(
double)*grid->
ndim );
3894 grid->
jlo = (
long *)
MALLOC(
sizeof(
long)*jsize );
3895 grid->
jhi = (
long *)
MALLOC(
sizeof(
long)*jsize );
3899 FillJ( grid, index, val, grid->
ndim, lgList );
3921 long n =
JIndex(grid,index);
3926 for( index[nd]=0; index[nd] < grid->
nval[nd]; index[nd]++ )
3928 val[nd] = grid->
val[nd][index[nd]];
3929 FillJ( grid, index, val, nd, lgList );
3933 if( lgList && nd ==
MIN2(grid->
ndim-1,1) )
3936 bool lgAge = ( strcmp( grid->
names[0],
"Age" ) == 0 );
3938 fprintf(
ioQQQ,
"\n" );
3939 if( grid->
ndim > 2 )
3941 fprintf(
ioQQQ,
"subgrid for" );
3942 for( n = nd+1; n < grid->
ndim; n++ )
3943 fprintf(
ioQQQ,
" %s=%g", grid->
names[n], val[n] );
3944 fprintf(
ioQQQ,
":\n\n" );
3946 if( grid->
ndim > 1 )
3948 fprintf(
ioQQQ,
"Teff\\lg g|" );
3949 for( n = 0; n < grid->
nval[1]; n++ )
3950 fprintf(
ioQQQ,
" %5.2f", grid->
val[1][n] );
3951 fprintf(
ioQQQ,
"\n" );
3952 fprintf(
ioQQQ,
"---------|" );
3953 for( n = 0; n < grid->
nval[1]; n++ )
3954 fprintf(
ioQQQ,
"------" );
3959 fprintf(
ioQQQ,
" Age |\n" );
3961 fprintf(
ioQQQ,
" Teff |\n" );
3962 fprintf(
ioQQQ,
"---------|------" );
3964 fprintf(
ioQQQ,
"\n" );
3965 for( index[0]=0; index[0] < grid->
nval[0]; index[0]++ )
3968 fprintf(
ioQQQ,
"%8.2e |", grid->
val[0][index[0]] );
3970 fprintf(
ioQQQ,
"%8.0f |", grid->
val[0][index[0]] );
3971 if( grid->
ndim > 1 )
3973 for( index[1]=0; index[1] < grid->
nval[1]; index[1]++ )
3978 fprintf(
ioQQQ,
" --" );
3984 fprintf(
ioQQQ,
"\n" );
3999 for( i=0; i < grid->
ndim; i++ )
4001 ind += index[i]*mul;
4002 mul *= grid->
nval[i];
4015 double alogg_low = -DBL_MAX, alogg_high = DBL_MAX;
4033 *index_low = *index_high = -2;
4034 for( i=0; i < nmods; i++ )
4036 bool lgNext =
false;
4038 for( nd=0; nd < ndim; nd++ )
4040 if( nd != 1 && fabs(telg[i].par[nd]-val[nd]) > 10.*DBL_EPSILON*fabs(val[nd]) )
4051 if( ndim == 1 || fabs(telg[i].par[1]-val[1]) <= 10.*DBL_EPSILON*fabs(val[1]) )
4058 if( telg[i].par[1] < val[1] && telg[i].par[1] > alogg_low )
4061 alogg_low = telg[i].
par[1];
4064 if( telg[i].par[1] > val[1] && telg[i].par[1] < alogg_high )
4067 alogg_high = telg[i].
par[1];
4080 bool lgOutLo, lgOutHi;
4100 lgOutLo = ( x-xval[0] < -10.*DBL_EPSILON*fabs(xval[0]) );
4101 lgOutHi = ( x-xval[NVAL-1] > 10.*DBL_EPSILON*fabs(xval[NVAL-1]) );
4103 if( lgOutLo || lgOutHi )
4109 *ind1 = lgOutLo ? -1 : NVAL-1;
4110 *ind2 = lgOutLo ? 0 : NVAL;
4122 for( i=0; i < NVAL; i++ )
4125 if( fabs(xval[i]-x) <= 10.*DBL_EPSILON*fabs(xval[i]) )
4134 for( i=0; i < NVAL-1; i++ )
4136 if( xval[i] < x && x < xval[i+1] )
4145 fprintf(
ioQQQ,
" insanity in FindIndex\n" );
4156 ioIN =
open_data( chFnam,
"r", scheme );
4178 if( strcmp( grid->
names[0],
"Teff" ) != 0 )
4188 for( i=0; i < grid->
nmods; i++ )
4190 fprintf(
ioQQQ,
"testing model %ld ", i+1 );
4191 for( nd=0; nd < grid->
npar; nd++ )
4197 fprintf(
ioQQQ,
" OK\n" );
4201 FILE *ioBUG = fopen(
"atmosphere_dump.txt", ( i == 0 ) ?
"w" :
"a" );
4203 fprintf( ioBUG,
"######## MODEL %ld", i+1 );
4204 for( nd=0; nd < grid->
npar; nd++ )
4205 fprintf( ioBUG,
" %s %g", grid->
names[nd], grid->
telg[i].
par[nd] );
4206 fprintf( ioBUG,
" ####################\n" );
4209 fprintf( ioBUG,
"%e %e\n", anu[k], flux[k] );
4225 bool lgPassed =
true;
4236 lumi += (anu[k] - anu[k-1])*(flux[k] + flux[k-1])/2.;
4241 if( fabs(Teff - chk) > toler*Teff ) {
4242 fprintf(
ioQQQ,
"\n*** WARNING, Teff discrepancy for this model, expected Teff %.2f, ", Teff);
4243 fprintf(
ioQQQ,
"integration yielded Teff %.2f, delta %.2f%%\n", chk, (chk/Teff-1.)*100. );
4280 for( j=0; j < nEdge; j++ )
4282 ind =
RebinFind(StarEner,nCont,AbsorbEdge[j]);
4285 ASSERT( ind >= 0 && ind+1 < nCont );
4287 EdgeLow[j] = StarEner[ind];
4288 EdgeHigh[j] = StarEner[ind+1];
4294 for( j=0; j < nCont; j++ )
4296 if( StarFlux[j] == 0.f )
4306 for( j=0; j < nCont-1; j++ )
4308 double ratio_x, ratio_y;
4311 ASSERT( StarEner[j+1] > StarEner[j] );
4315 ratio_x = (double)StarEner[j+1]/(
double)StarEner[j];
4316 ratio_y = (double)StarFlux[j+1]/(
double)StarFlux[j];
4317 StarPower[j] = (
realnum)(log(ratio_y)/log(ratio_x));
4324 BinLow = ( j > 0 ) ?
4340 for( k=0; k < nEdge; k++ )
4342 if( BinLow < EdgeLow[k] && BinNext > EdgeHigh[k] )
4344 BinMid = 0.99999f*EdgeLow[k];
4345 CloudyFlux[j] =
RebinSingleCell(BinLow,BinMid,StarEner,StarFlux,StarPower,nCont);
4351 BinMid = 1.00001f*EdgeHigh[k];
4352 CloudyFlux[j] =
RebinSingleCell(BinMid,BinNext,StarEner,StarFlux,StarPower,nCont);
4361 CloudyFlux[j] =
RebinSingleCell(BinLow,BinHigh,StarEner,StarFlux,StarPower,nCont);
4393 anu = sqrt(BinLow*BinHigh);
4395 widflx =
MIN2(BinHigh,StarEner[nCont-1])-BinLow;
4397 if( BinLow < StarEner[0] )
4401 retval = (
realnum)(StarFlux[0]*pow(anu/StarEner[0],2.));
4403 else if( BinLow > StarEner[nCont-1] )
4412 ipLo =
RebinFind(StarEner,nCont,BinLow);
4413 ipHi =
RebinFind(StarEner,nCont,BinHigh);
4415 ASSERT( ipLo >= 0 && ipLo < nCont-1 && ipHi >= ipLo );
4421 retval = (
realnum)(StarFlux[ipLo]*pow(anu/StarEner[ipLo],(
double)StarPower[ipLo]));
4435 for( i=ipLo; i <=
MIN2(ipHi,nCont-2); i++ )
4437 double pp1 = StarPower[i] + 1.;
4443 v1 = StarFlux[i]*pow(x1/StarEner[i],(
double)StarPower[i]);
4447 else if( i == ipHi )
4464 if( fabs(pp1) < 0.001 )
4466 val = x1*v1*log(x2/x1);
4470 val = pow(x2/x1,pp1) - 1.;
4471 val = val*x1*v1/pp1;
4476 retval = sum/widflx;
4504 if( val < array[0] )
4508 else if( val > array[nArr-1] )
4520 sgn =
sign3(val-array[i2]);