14 static const int MNTS = 200;
58 mpp() { memset(
this, 0,
sizeof(
mpp)); }
181 long,
const double[],
int);
182 STATIC void RauchReadMPP(vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&);
184 STATIC void WriteASCIIHead(FILE*,
long,
long,
const vector<string>&,
long,
long,
const string&,
double,
185 const string&,
double,
const vector<mpp>&,
const char*,
int);
200 vector<long>&,
long,vector<realnum>&,
bool,
const realnum[]=NULL,
long=0L);
202 const vector<long>&,vector<long>&,
long,vector<realnum>&,
int);
204 const long[],vector<long>&,
long,
long,vector<realnum>&);
205 template<
class T>
void SortUnique(vector<T>&,vector<T>&);
209 const vector<realnum>&,
double*,
double*);
211 long,
double*,
double*);
215 STATIC void SearchModel(
const vector<mpp>&,
bool,
long,
const vector<double>&,
long,
long*,
long*);
225 long,
const T[],
const realnum[]);
265 fprintf(
ioQQQ,
"\n I will now list all stellar atmosphere grids that are ready to be used (if any).\n" );
266 fprintf(
ioQQQ,
" User-defined stellar atmosphere grids will not be included in this list.\n\n" );
275 fprintf(
ioQQQ,
" table star atlas Z+1.0 <Teff> [ <log(g)> ]\n" );
277 fprintf(
ioQQQ,
" table star atlas Z+0.5 <Teff> [ <log(g)> ]\n" );
279 fprintf(
ioQQQ,
" table star atlas Z+0.3 <Teff> [ <log(g)> ]\n" );
281 fprintf(
ioQQQ,
" table star atlas Z+0.2 <Teff> [ <log(g)> ]\n" );
283 fprintf(
ioQQQ,
" table star atlas Z+0.1 <Teff> [ <log(g)> ]\n" );
285 fprintf(
ioQQQ,
" table star atlas Z+0.0 <Teff> [ <log(g)> ]\n" );
287 fprintf(
ioQQQ,
" table star atlas Z-0.1 <Teff> [ <log(g)> ]\n" );
289 fprintf(
ioQQQ,
" table star atlas Z-0.2 <Teff> [ <log(g)> ]\n" );
291 fprintf(
ioQQQ,
" table star atlas Z-0.3 <Teff> [ <log(g)> ]\n" );
293 fprintf(
ioQQQ,
" table star atlas Z-0.5 <Teff> [ <log(g)> ]\n" );
295 fprintf(
ioQQQ,
" table star atlas Z-1.0 <Teff> [ <log(g)> ]\n" );
297 fprintf(
ioQQQ,
" table star atlas Z-1.5 <Teff> [ <log(g)> ]\n" );
299 fprintf(
ioQQQ,
" table star atlas Z-2.0 <Teff> [ <log(g)> ]\n" );
301 fprintf(
ioQQQ,
" table star atlas Z-2.5 <Teff> [ <log(g)> ]\n" );
303 fprintf(
ioQQQ,
" table star atlas Z-3.0 <Teff> [ <log(g)> ]\n" );
305 fprintf(
ioQQQ,
" table star atlas Z-3.5 <Teff> [ <log(g)> ]\n" );
307 fprintf(
ioQQQ,
" table star atlas Z-4.0 <Teff> [ <log(g)> ]\n" );
309 fprintf(
ioQQQ,
" table star atlas Z-4.5 <Teff> [ <log(g)> ]\n" );
311 fprintf(
ioQQQ,
" table star atlas Z-5.0 <Teff> [ <log(g)> ]\n" );
314 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.5 <Teff> [ <log(g)> ]\n" );
316 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.2 <Teff> [ <log(g)> ]\n" );
318 fprintf(
ioQQQ,
" table star atlas odfnew Z+0.0 <Teff> [ <log(g)> ]\n" );
320 fprintf(
ioQQQ,
" table star atlas odfnew Z-0.5 <Teff> [ <log(g)> ]\n" );
322 fprintf(
ioQQQ,
" table star atlas odfnew Z-1.0 <Teff> [ <log(g)> ]\n" );
324 fprintf(
ioQQQ,
" table star atlas odfnew Z-1.5 <Teff> [ <log(g)> ]\n" );
326 fprintf(
ioQQQ,
" table star atlas odfnew Z-2.0 <Teff> [ <log(g)> ]\n" );
328 fprintf(
ioQQQ,
" table star atlas odfnew Z-2.5 <Teff> [ <log(g)> ]\n" );
331 fprintf(
ioQQQ,
" table star atlas 3-dim <Teff> <log(g)> <log(Z)>\n" );
334 fprintf(
ioQQQ,
" table star atlas odfnew 3-dim <Teff> <log(g)> <log(Z)>\n" );
337 fprintf(
ioQQQ,
" table star costar solar (see Hazy for parameters)\n" );
339 fprintf(
ioQQQ,
" table star costar halo (see Hazy for parameters)\n" );
348 fprintf(
ioQQQ,
" table star rauch H-Ca solar <Teff> [ <log(g)> ]\n" );
350 fprintf(
ioQQQ,
" table star rauch H-Ca halo <Teff> [ <log(g)> ]\n" );
352 fprintf(
ioQQQ,
" table star rauch H-Ca 3-dim <Teff> <log(g)> <log(Z)>\n" );
355 fprintf(
ioQQQ,
" table star rauch H-Ni solar <Teff> [ <log(g)> ]\n" );
357 fprintf(
ioQQQ,
" table star rauch H-Ni halo <Teff> [ <log(g)> ]\n" );
359 fprintf(
ioQQQ,
" table star rauch H-Ni 3-dim <Teff> <log(g)> <log(Z)>\n" );
362 fprintf(
ioQQQ,
" table star rauch pg1159 <Teff> [ <log(g)> ]\n" );
367 fprintf(
ioQQQ,
" table star rauch hydrogen <Teff> [ <log(g)> ]\n" );
370 fprintf(
ioQQQ,
" table star rauch helium <Teff> [ <log(g)> ]\n" );
373 fprintf(
ioQQQ,
" table star rauch H+He <Teff> <log(g)> <frac(He)>\n" );
376 fprintf(
ioQQQ,
" table star \"starburst99.mod\" <age>\n" );
378 fprintf(
ioQQQ,
" table star \"starburst99_2d.mod\" <age> <Z>\n" );
381 fprintf(
ioQQQ,
" table star tlusty OBstar Z+0.3 <Teff> [ <log(g)> ]\n" );
383 fprintf(
ioQQQ,
" table star tlusty OBstar Z+0.0 <Teff> [ <log(g)> ]\n" );
385 fprintf(
ioQQQ,
" table star tlusty OBstar Z-0.3 <Teff> [ <log(g)> ]\n" );
387 fprintf(
ioQQQ,
" table star tlusty OBstar Z-0.7 <Teff> [ <log(g)> ]\n" );
389 fprintf(
ioQQQ,
" table star tlusty OBstar Z-1.0 <Teff> [ <log(g)> ]\n" );
391 fprintf(
ioQQQ,
" table star tlusty OBstar Z-inf <Teff> [ <log(g)> ]\n" );
394 fprintf(
ioQQQ,
" table star tlusty OBstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
397 fprintf(
ioQQQ,
" table star tlusty Bstar Z+0.3 <Teff> [ <log(g)> ]\n" );
399 fprintf(
ioQQQ,
" table star tlusty Bstar Z+0.0 <Teff> [ <log(g)> ]\n" );
401 fprintf(
ioQQQ,
" table star tlusty Bstar Z-0.3 <Teff> [ <log(g)> ]\n" );
403 fprintf(
ioQQQ,
" table star tlusty Bstar Z-0.7 <Teff> [ <log(g)> ]\n" );
405 fprintf(
ioQQQ,
" table star tlusty Bstar Z-1.0 <Teff> [ <log(g)> ]\n" );
407 fprintf(
ioQQQ,
" table star tlusty Bstar Z-inf <Teff> [ <log(g)> ]\n" );
410 fprintf(
ioQQQ,
" table star tlusty Bstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
413 fprintf(
ioQQQ,
" table star tlusty Ostar Z+0.3 <Teff> [ <log(g)> ]\n" );
415 fprintf(
ioQQQ,
" table star tlusty Ostar Z+0.0 <Teff> [ <log(g)> ]\n" );
417 fprintf(
ioQQQ,
" table star tlusty Ostar Z-0.3 <Teff> [ <log(g)> ]\n" );
419 fprintf(
ioQQQ,
" table star tlusty Ostar Z-0.7 <Teff> [ <log(g)> ]\n" );
421 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.0 <Teff> [ <log(g)> ]\n" );
423 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.5 <Teff> [ <log(g)> ]\n" );
425 fprintf(
ioQQQ,
" table star tlusty Ostar Z-1.7 <Teff> [ <log(g)> ]\n" );
427 fprintf(
ioQQQ,
" table star tlusty Ostar Z-2.0 <Teff> [ <log(g)> ]\n" );
429 fprintf(
ioQQQ,
" table star tlusty Ostar Z-3.0 <Teff> [ <log(g)> ]\n" );
431 fprintf(
ioQQQ,
" table star tlusty Ostar Z-inf <Teff> [ <log(g)> ]\n" );
434 fprintf(
ioQQQ,
" table star tlusty Ostar 3-dim <Teff> <log(g)> <log(Z)>\n" );
437 fprintf(
ioQQQ,
" table star werner <Teff> [ <log(g)> ]\n" );
440 fprintf(
ioQQQ,
" table star wmbasic <Teff> <log(g)> <log(Z)>\n" );
445 fprintf(
ioQQQ,
" table HM05 quasar <z> [ <factor> ]\n" );
471 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp10k2.ascii",
"atlas_fp10k2.mod", NULL, 0L, pc );
473 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp05k2.ascii",
"atlas_fp05k2.mod", NULL, 0L, pc );
475 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp03k2.ascii",
"atlas_fp03k2.mod", NULL, 0L, pc );
477 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp02k2.ascii",
"atlas_fp02k2.mod", NULL, 0L, pc );
479 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp01k2.ascii",
"atlas_fp01k2.mod", NULL, 0L, pc );
481 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp00k2.ascii",
"atlas_fp00k2.mod", NULL, 0L, pc );
483 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm01k2.ascii",
"atlas_fm01k2.mod", NULL, 0L, pc );
485 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm02k2.ascii",
"atlas_fm02k2.mod", NULL, 0L, pc );
487 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm03k2.ascii",
"atlas_fm03k2.mod", NULL, 0L, pc );
489 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm05k2.ascii",
"atlas_fm05k2.mod", NULL, 0L, pc );
491 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm10k2.ascii",
"atlas_fm10k2.mod", NULL, 0L, pc );
493 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm15k2.ascii",
"atlas_fm15k2.mod", NULL, 0L, pc );
495 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm20k2.ascii",
"atlas_fm20k2.mod", NULL, 0L, pc );
497 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm25k2.ascii",
"atlas_fm25k2.mod", NULL, 0L, pc );
499 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm30k2.ascii",
"atlas_fm30k2.mod", NULL, 0L, pc );
501 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm35k2.ascii",
"atlas_fm35k2.mod", NULL, 0L, pc );
503 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm40k2.ascii",
"atlas_fm40k2.mod", NULL, 0L, pc );
505 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm45k2.ascii",
"atlas_fm45k2.mod", NULL, 0L, pc );
507 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm50k2.ascii",
"atlas_fm50k2.mod", NULL, 0L, pc );
512 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp05k2_odfnew.ascii",
"atlas_fp05k2_odfnew.mod",
516 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp02k2_odfnew.ascii",
"atlas_fp02k2_odfnew.mod",
520 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fp00k2_odfnew.ascii",
"atlas_fp00k2_odfnew.mod",
524 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm05k2_odfnew.ascii",
"atlas_fm05k2_odfnew.mod",
528 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm10k2_odfnew.ascii",
"atlas_fm10k2_odfnew.mod",
532 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm15k2_odfnew.ascii",
"atlas_fm15k2_odfnew.mod",
536 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm20k2_odfnew.ascii",
"atlas_fm20k2_odfnew.mod",
540 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_fm25k2_odfnew.ascii",
"atlas_fm25k2_odfnew.mod",
548 lgFail = lgFail ||
lgCompileAtmosphere(
"atlas_3d_odfnew.ascii",
"atlas_3d_odfnew.mod", NULL, 0L, pc );
556 const char *chMetalicity,
557 const char *chODFNew,
565 grid.
name =
"atlas_";
571 grid.
name += chMetalicity;
574 grid.
name += chODFNew;
582 strcpy( chIdent,
"3-dim" );
586 strcpy( chIdent,
"Z " );
587 strcat( chIdent, chMetalicity );
589 strcat( chIdent, ( strlen(chODFNew) == 0 ?
" Kurucz" :
" ODFNew" ) );
590 grid.
ident = chIdent;
592 grid.
command =
"COMPILE STARS";
633 fprintf(
ioQQQ,
" Creating Sc1_costar_solar.ascii....\n" );
634 lgFail = lgFail ||
CoStarInitialize(
"Sc1_costar_z020_lb.fluxes",
"Sc1_costar_solar.ascii" );
638 fprintf(
ioQQQ,
" Creating Sc1_costar_halo.ascii....\n" );
639 lgFail = lgFail ||
CoStarInitialize(
"Sc1_costar_z004_lb.fluxes",
"Sc1_costar_halo.ascii" );
665 grid.
name = ( lgHalo ?
"Sc1_costar_halo.mod" :
"Sc1_costar_solar.mod" );
669 grid.
ident =
" costar";
671 grid.
command =
"COMPILE STARS";
728 string OutName( InName );
729 string::size_type ptr = OutName.find(
'.' );
730 ASSERT( ptr != string::npos );
731 OutName.replace( ptr, string::npos,
".mod" );
743 grid.
ident =
"bogus ident.";
744 grid.
command =
"bogus command.";
750 if( strcmp( grid.
names[0],
"Teff" ) == 0 )
752 fprintf(
ioQQQ,
" GridCompile: checking effective temperatures...\n" );
763 const char *FileName,
771 string chTruncName( FileName );
772 string::size_type ptr = chTruncName.find(
'.' );
773 if( ptr != string::npos )
774 chTruncName.replace( ptr, string::npos,
"" );
777 grid.
name = FileName;
782 grid.
ident = chTruncName.substr(0,12);
783 if( grid.
ident.length() < 12 )
784 grid.
ident.append(12-grid.
ident.length(),
' ');
787 grid.
command =
"COMPILE STARS \"" + chTruncName +
".ascii\"";
809 if( version == 2005 )
814 grid.
ident =
" HM05 ";
816 else if( version == 2012 )
819 grid.
ident =
" HM12 ";
825 grid.
name +=
"quasar.ascii";
826 grid.
ident +=
"QUASAR";
830 grid.
name +=
"galaxy.ascii";
831 grid.
ident +=
"GALAXY";
841 vector<realnum> Edges;
842 Edges.push_back(
realnum(RYDLAM/1216.0) );
843 if( version == 2012 )
845 Edges.push_back(
realnum(RYDLAM/1026.0) );
846 Edges.push_back(
realnum(RYDLAM/972.5) );
847 Edges.push_back(
realnum(RYDLAM/949.7) );
848 Edges.push_back(
realnum(RYDLAM/937.8) );
851 Edges.push_back(
realnum(RYDLAM/304.0) );
852 if( version == 2012 )
854 Edges.push_back(
realnum(RYDLAM/256.4) );
855 Edges.push_back(
realnum(RYDLAM/243.1) );
856 Edges.push_back(
realnum(RYDLAM/237.4) );
857 Edges.push_back(
realnum(RYDLAM/234.4) );
859 Edges.push_back(
realnum(RYDLAM/228.0) );
863 long nval = 1, ndim = 1;
864 CheckVal( &grid, &val, &nval, &ndim );
900 grid.
name =
"kurucz79.mod";
904 grid.
ident =
" Kurucz 1979";
906 grid.
command =
"COMPILE STARS";
929 Edges[0] = (
realnum)(RYDLAM/911.204);
930 Edges[1] = (
realnum)(RYDLAM/503.968);
931 Edges[2] = (
realnum)(RYDLAM/227.806);
952 grid.
name =
"mihalas.mod";
956 grid.
ident =
" Mihalas";
958 grid.
command =
"COMPILE STARS";
975 static const double par2[2] = { 0., -1. };
978 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 };
997 RauchReadMPP( telg1, telg2, telg3, telg4, telg5, telg6 );
1001 bool lgFail =
false;
1006 fprintf(
ioQQQ,
" Creating rauch_h-ca_solar.ascii....\n" );
1007 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ca_solar.ascii",
"_solar_bin_0.1",
1013 fprintf(
ioQQQ,
" Creating rauch_h-ca_halo.ascii....\n" );
1014 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ca_halo.ascii",
"_halo__bin_0.1",
1022 fprintf(
ioQQQ,
" Creating rauch_h-ca_3d.ascii....\n" );
1023 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ca_3d.ascii",
"_solar_bin_0.1",
1025 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ca_3d.ascii",
"_halo__bin_0.1",
1030 if(
lgFileReadable(
"0050000_50_solar_iron.bin_0.1", dum, as ) &&
1033 fprintf(
ioQQQ,
" Creating rauch_h-ni_solar.ascii....\n" );
1034 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ni_solar.ascii",
"_solar_iron.bin_0.1",
1038 if(
lgFileReadable(
"0050000_50_halo__iron.bin_0.1", dum, as ) &&
1041 fprintf(
ioQQQ,
" Creating rauch_h-ni_halo.ascii....\n" );
1042 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ni_halo.ascii",
"_halo__iron.bin_0.1",
1046 if(
lgFileReadable(
"0050000_50_solar_iron.bin_0.1", dum, as ) &&
1050 fprintf(
ioQQQ,
" Creating rauch_h-ni_3d.ascii....\n" );
1051 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ni_3d.ascii",
"_solar_iron.bin_0.1",
1053 lgFail = lgFail ||
RauchInitialize(
"rauch_h-ni_3d.ascii",
"_halo__iron.bin_0.1",
1058 if(
lgFileReadable(
"0040000_5.00_33_50_02_15.bin_0.1", dum, as ) &&
1061 fprintf(
ioQQQ,
" Creating rauch_pg1159.ascii....\n" );
1062 lgFail = lgFail ||
RauchInitialize(
"rauch_pg1159.ascii",
"_33_50_02_15.bin_0.1",
1067 if(
lgFileReadable(
"0020000_4.00_H_00005-02000A.bin_0.1", dum, as ) &&
1071 lgFail = lgFail ||
RauchInitialize(
"rauch_hydr.ascii",
"_H_00005-02000A.bin_0.1",
1076 if(
lgFileReadable(
"0050000_5.00_He_00005-02000A.bin_0.1", dum, as ) &&
1079 fprintf(
ioQQQ,
" Creating rauch_helium.ascii....\n" );
1080 lgFail = lgFail ||
RauchInitialize(
"rauch_helium.ascii",
"_He_00005-02000A.bin_0.1",
1085 if(
lgFileReadable(
"0050000_5.00_H+He_1.000_0.000_00005-02000A.bin_0.1", dum, as ) &&
1088 fprintf(
ioQQQ,
" Creating rauch_h+he_3d.ascii....\n" );
1089 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_1.000_0.000_00005-02000A.bin_0.1",
1091 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.900_0.100_00005-02000A.bin_0.1",
1093 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.800_0.200_00005-02000A.bin_0.1",
1095 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.700_0.300_00005-02000A.bin_0.1",
1097 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.600_0.400_00005-02000A.bin_0.1",
1099 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.500_0.500_00005-02000A.bin_0.1",
1101 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.400_0.600_00005-02000A.bin_0.1",
1103 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.300_0.700_00005-02000A.bin_0.1",
1105 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.200_0.800_00005-02000A.bin_0.1",
1107 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.100_0.900_00005-02000A.bin_0.1",
1109 lgFail = lgFail ||
RauchInitialize(
"rauch_h+he_3d.ascii",
"_H+He_0.000_1.000_00005-02000A.bin_0.1",
1114 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_solar.ascii",
"rauch_h-ca_solar.mod", NULL,0L, pc );
1116 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_halo.ascii",
"rauch_h-ca_halo.mod", NULL, 0L, pc );
1118 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ca_3d.ascii",
"rauch_h-ca_3d.mod", NULL, 0L, pc );
1121 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_solar.ascii",
"rauch_h-ni_solar.mod", NULL,0L, pc );
1123 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_halo.ascii",
"rauch_h-ni_halo.mod", NULL, 0L, pc );
1125 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h-ni_3d.ascii",
"rauch_h-ni_3d.mod", NULL, 0L, pc );
1128 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_pg1159.ascii",
"rauch_pg1159.mod", NULL, 0L, pc );
1130 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_cowd.ascii",
"rauch_cowd.mod", NULL, 0L, pc );
1133 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_hydr.ascii",
"rauch_hydr.mod", NULL, 0L, pc );
1136 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_helium.ascii",
"rauch_helium.mod", NULL, 0L, pc );
1139 lgFail = lgFail ||
lgCompileAtmosphere(
"rauch_h+he_3d.ascii",
"rauch_h+he_3d.mod", NULL, 0L, pc );
1156 grid.
name =
"rauch_h-ca_3d.mod";
1158 grid.
name = ( lgHalo ?
"rauch_h-ca_halo.mod" :
"rauch_h-ca_solar.mod" );
1162 grid.
ident =
" H-Ca Rauch";
1164 grid.
command =
"COMPILE STARS";
1168 CheckVal( &grid, val, nval, ndim );
1188 grid.
name =
"rauch_h-ni_3d.mod";
1190 grid.
name = ( lgHalo ?
"rauch_h-ni_halo.mod" :
"rauch_h-ni_solar.mod" );
1194 grid.
ident =
" H-Ni Rauch";
1196 grid.
command =
"COMPILE STARS";
1200 CheckVal( &grid, val, nval, ndim );
1218 grid.
name =
"rauch_pg1159.mod";
1222 grid.
ident =
"PG1159 Rauch";
1224 grid.
command =
"COMPILE STARS";
1228 CheckVal( &grid, val, nval, ndim );
1246 grid.
name =
"rauch_cowd.mod";
1250 grid.
ident =
"C/O WD Rauch";
1252 grid.
command =
"COMPILE STARS";
1256 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 );
1302 grid.
name =
"rauch_helium.mod";
1306 grid.
ident =
"Helium Rauch";
1308 grid.
command =
"COMPILE STARS";
1312 CheckVal( &grid, val, nval, ndim );
1330 grid.
name =
"rauch_h+he_3d.mod";
1334 grid.
ident =
" H+He Rauch";
1336 grid.
command =
"COMPILE STARS";
1340 CheckVal( &grid, val, nval, ndim );
1349 const char chOutName[],
1355 vector<mpp> telg(
MNTS);
1356 vector<double> wavl, fluxes[
MNTS];
1357 wavl.reserve(
NSB99);
1365 bool lgHeader =
true;
1373 double cage, cwavl, cfl, cfl1, cfl2, cfl3;
1374 if( sscanf( chLine,
" %le %le %le %le %le", &cage, &cwavl, &cfl1, &cfl2, &cfl3 ) != 5 )
1376 fprintf(
ioQQQ,
"syntax error in data of Starburst grid.\n" );
1396 fprintf(
ioQQQ,
"too many time steps in Starburst grid.\n" );
1397 fprintf(
ioQQQ,
"please increase MNTS and recompile.\n" );
1405 fluxes[nmods].reserve(
NSB99);
1407 fluxes[nmods].reserve(fluxes[nmods-1].size());
1408 telg[nmods].par[0] = cage;
1411 if( !
fp_equal(telg[nmods].par[0],cage,10) )
1418 wavl.push_back(cwavl);
1421 if( !
fp_equal(wavl[ngp],cwavl,10) )
1423 fprintf(
ioQQQ,
"wavelength error in Starburst grid.\n" );
1430 fluxes[nmods].push_back(
exp10(cfl - 44.077911) );
1436 if( lgHeader && strncmp( &chLine[1],
"TIME [YR]", 9 ) == 0 )
1443 fprintf(
ioQQQ,
"syntax error in header of Starburst grid.\n" );
1453 FILE *ioOut =
open_data( chOutName,
"w" );
1455 vector<string> names;
1456 names.push_back(
"Age" );
1457 WriteASCIIHead(ioOut,
VERSION_ASCII, 1, names, nmods, ngp,
"lambda", 1.,
"F_lambda", 1., telg,
" %.3e", 4);
1459 for(
long i=0; i < nmods; i++ )
1476 bool lgFail =
false;
1483 lgFail = lgFail ||
lgCompileAtmosphere(
"starburst99.ascii",
"starburst99.mod", NULL, 0L, pc );
1486 lgFail = lgFail ||
lgCompileAtmosphere(
"starburst99_2d.ascii",
"starburst99_2d.mod", NULL, 0L, pc );
1499 bool lgFail =
false;
1501 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_p03.ascii",
"obstar_merged_p03.mod",NULL,0L,pc);
1503 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_p00.ascii",
"obstar_merged_p00.mod",NULL,0L,pc);
1505 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m03.ascii",
"obstar_merged_m03.mod",NULL,0L,pc);
1507 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m07.ascii",
"obstar_merged_m07.mod",NULL,0L,pc);
1509 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m10.ascii",
"obstar_merged_m10.mod",NULL,0L,pc);
1511 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_m99.ascii",
"obstar_merged_m99.mod",NULL,0L,pc);
1514 lgFail = lgFail ||
lgCompileAtmosphere(
"obstar_merged_3d.ascii",
"obstar_merged_3d.mod", NULL, 0L, pc);
1517 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_p03.ascii",
"bstar2006_p03.mod", NULL, 0L, pc );
1519 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_p00.ascii",
"bstar2006_p00.mod", NULL, 0L, pc );
1521 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m03.ascii",
"bstar2006_m03.mod", NULL, 0L, pc );
1523 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m07.ascii",
"bstar2006_m07.mod", NULL, 0L, pc );
1525 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m10.ascii",
"bstar2006_m10.mod", NULL, 0L, pc );
1527 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_m99.ascii",
"bstar2006_m99.mod", NULL, 0L, pc );
1530 lgFail = lgFail ||
lgCompileAtmosphere(
"bstar2006_3d.ascii",
"bstar2006_3d.mod", NULL, 0L, pc );
1533 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_p03.ascii",
"ostar2002_p03.mod", NULL, 0L, pc );
1535 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_p00.ascii",
"ostar2002_p00.mod", NULL, 0L, pc );
1537 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m03.ascii",
"ostar2002_m03.mod", NULL, 0L, pc );
1539 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m07.ascii",
"ostar2002_m07.mod", NULL, 0L, pc );
1541 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m10.ascii",
"ostar2002_m10.mod", NULL, 0L, pc );
1543 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m15.ascii",
"ostar2002_m15.mod", NULL, 0L, pc );
1545 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m17.ascii",
"ostar2002_m17.mod", NULL, 0L, pc );
1547 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m20.ascii",
"ostar2002_m20.mod", NULL, 0L, pc );
1549 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m30.ascii",
"ostar2002_m30.mod", NULL, 0L, pc );
1551 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_m99.ascii",
"ostar2002_m99.mod", NULL, 0L, pc );
1554 lgFail = lgFail ||
lgCompileAtmosphere(
"ostar2002_3d.ascii",
"ostar2002_3d.mod", NULL, 0L, pc );
1563 const char *chMetalicity,
1572 grid.
name =
"obstar_merged_";
1574 grid.
name =
"bstar2006_";
1576 grid.
name =
"ostar2002_";
1582 grid.
name += chMetalicity;
1583 grid.
name +=
".mod";
1590 strcpy( chIdent,
"3-dim" );
1594 strcpy( chIdent,
"Z " );
1595 strcat( chIdent, chMetalicity );
1598 strcat( chIdent,
" OBstar" );
1600 strcat( chIdent,
" Bstr06" );
1602 strcat( chIdent,
" Ostr02" );
1605 grid.
ident = chIdent;
1607 grid.
command =
"COMPILE STARS";
1611 CheckVal( &grid, val, nval, ndim );
1628 realnum Edges[3] = { 0.99946789f, 1.8071406f, 3.9996377f };
1652 bool lgFail =
false;
1675 grid.
name =
"kwerner.mod";
1679 grid.
ident =
"Klaus Werner";
1681 grid.
command =
"COMPILE STARS";
1685 CheckVal( &grid, val, nval, ndim );
1721 realnum Edges[3] = { 0.9994665f, 1.807140f, 3.999632f };
1725 bool lgFail =
false;
1742 grid.
name =
"wmbasic.mod";
1746 grid.
ident =
" WMBASIC";
1748 grid.
command =
"COMPILE STARS";
1752 CheckVal( &grid, val, nval, ndim );
1761 const char chFNameOut[])
1782 fprintf(
ioQQQ,
" CoStarInitialize fails reading nskip.\n" );
1785 sscanf( chLine,
"%li", &nskip );
1788 for(
long i=0; i < nskip; ++i )
1792 fprintf(
ioQQQ,
" CoStarInitialize fails skipping header.\n" );
1801 fprintf(
ioQQQ,
" CoStarInitialize fails reading nModels, nWL.\n" );
1804 sscanf( chLine,
"%li%li", &nModels, &nWL );
1806 if( nModels <= 0 || nWL <= 0 )
1808 fprintf(
ioQQQ,
" CoStarInitialize scanned off impossible values for nModels=%li or nWL=%li\n",
1814 vector<mpp> telg(nModels);
1817 for(
long i=0; i < nModels; ++i )
1821 fprintf(
ioQQQ,
" CoStarInitialize fails reading model parameters.\n" );
1825 telg[i].chGrid = chLine[0];
1827 sscanf( chLine+1,
"%i", &telg[i].modid );
1829 sscanf( chLine+23,
"%lg", &telg[i].par[0] );
1831 sscanf( chLine+31,
"%lg", &telg[i].par[1] );
1833 sscanf( chLine+7,
"%lg", &telg[i].par[2] );
1835 sscanf( chLine+15,
"%lg", &telg[i].par[3] );
1838 ASSERT( telg[i].par[2] > 10. );
1839 ASSERT( telg[i].par[3] > 10. );
1842 telg[i].par[2] = log10(telg[i].par[2]);
1857 vector<string> names;
1858 names.push_back(
"Teff" );
1859 names.push_back(
"log(g)" );
1860 names.push_back(
"log(M)" );
1861 names.push_back(
"Age" );
1862 WriteASCIIHead(ioOUT,
VERSION_COSTAR, 2, names, nModels, nWL-1,
"lambda", 1.,
"F_nu", PI, telg,
" %.6e", 1);
1865 vector<double> StarWavl(nWL);
1866 vector<double> StarFlux(nWL);
1869 for(
long i=0; i < nModels; ++i )
1874 fprintf(
ioQQQ,
" CoStarInitialize fails reading the skip to next spectrum.\n" );
1877 sscanf( chLine,
"%li", &nskip );
1879 for(
long j=0; j < nskip; ++j )
1883 fprintf(
ioQQQ,
" CoStarInitialize fails doing the skip.\n" );
1889 for(
long j=0; j < nWL; ++j )
1893 fprintf(
ioQQQ,
" CoStarInitialize fails reading the spectral data.\n" );
1896 double help1, help2;
1897 sscanf( chLine,
"%lg %lg", &help1, &help2 );
1901 StarWavl[j] = help1;
1905 StarFlux[j] =
exp10(help2);
1909 ASSERT( StarWavl[j] > StarWavl[j-1] );
1937 switch( grid->
imode )
1946 lval[0] = log10(val[0]);
1952 lval[0] = log10(val[1]);
1957 fprintf(
ioQQQ,
" InterpolateGridCoStar called with insane value for imode: %d.\n", grid->
imode );
1961 long nmodid = (long)(lval[1]+0.5);
1987 vector<realnum> ValTr(grid->
nTracks);
1988 vector<long> indloTr(grid->
nTracks);
1989 vector<long> indhiTr(grid->
nTracks);
1990 vector<long> index(2);
1993 for(
long j=0; j < grid->
nTracks; j++ )
1997 if( grid->
trackLen[j] >= nmodid ) {
1998 index[0] = nmodid - 1;
2009 ValTr[j] = -FLT_MAX;
2014 FindHCoStar( grid, j, lval[1], off, ValTr, indloTr, indhiTr );
2020 for(
long j=0; j < grid->
nTracks; j++ )
2022 if( indloTr[j] >= 0 )
2023 printf(
"track %c: models %c%d, %c%d, val %g\n",
2024 (
char)(
'A'+j), grid->
telg[indloTr[j]].chGrid, grid->
telg[indloTr[j]].modid,
2025 grid->
telg[indhiTr[j]].chGrid, grid->
telg[indhiTr[j]].modid, ValTr[j]);
2029 long useTr[2], indlo[2], indhi[2];
2040 fprintf(
ioQQQ,
" The parameters for the requested CoStar model are out of range.\n" );
2046 ASSERT( indloTr[useTr[0]] >= 0 && indloTr[useTr[0]] < (
int)grid->
nmods );
2047 ASSERT( indhiTr[useTr[0]] >= 0 && indhiTr[useTr[0]] < (
int)grid->
nmods );
2048 ASSERT( indloTr[useTr[1]] >= 0 && indloTr[useTr[1]] < (
int)grid->
nmods );
2049 ASSERT( indhiTr[useTr[1]] >= 0 && indhiTr[useTr[1]] < (
int)grid->
nmods );
2052 printf(
"interpolate between tracks %c and %c\n", (
char)(
'A'+useTr[0]), (
char)(
'A'+useTr[1]) );
2054 indlo[0] = indloTr[useTr[0]];
2055 indhi[0] = indhiTr[useTr[0]];
2056 indlo[1] = indloTr[useTr[1]];
2057 indhi[1] = indhiTr[useTr[1]];
2071 fprintf(
ioQQQ,
"Failed to read atmosphere models.\n" );
2075 for(
size_t i=0; i < grid->
index_list2.size(); ++i )
2078 vector<double> aval(4);
2090 FILE *ioBUG =
open_data(
"interpolated.txt",
"w" );
2102 SetLimits( grid, lval[0], dum, dum, useTr, ValTr, val0_lo, val0_hi );
2107 fprintf(
ioQQQ,
" * #<< FINAL: T_eff = %7.1f, ", aval[0] );
2108 fprintf(
ioQQQ,
"log(g) = %4.2f, M(ZAMS) = %5.1f, age = ", aval[1],
exp10(aval[2]) );
2119 vector<realnum>& ValTr,
2120 vector<long>& indloTr,
2121 vector<long>& indhiTr)
2125 indloTr[track] = -2;
2126 indhiTr[track] = -2;
2127 ValTr[track] = -FLT_MAX;
2130 vector<long> index(2);
2133 for(
long j=0; j < grid->
trackLen[track]; j++ )
2139 if( fabs(par2-grid->
telg[mod1].par[off+1]) <= 10.*FLT_EPSILON*fabs(grid->
telg[mod1].par[off+1]) )
2141 indloTr[track] = mod1;
2142 indhiTr[track] = mod1;
2143 ValTr[track] = (
realnum)grid->
telg[mod1].par[off];
2148 for(
long j=0; j < grid->
trackLen[track]-1; j++ )
2156 if( (par2 - grid->
telg[mod1].par[off+1])*(par2 - grid->
telg[mod2].par[off+1]) < 0. )
2160 indloTr[track] = mod1;
2161 indhiTr[track] = mod2;
2162 frac = (par2 - grid->
telg[mod2].par[off+1])/
2163 (grid->
telg[mod1].par[off+1] - grid->
telg[mod2].par[off+1]);
2164 ValTr[track] = (
realnum)(frac*grid->
telg[mod1].par[off] +
2165 (1.-frac)*grid->
telg[mod2].par[off] );
2174 vector<realnum>& ValTr,
2186 for(
long j=0; j < grid->
nTracks; j++ )
2189 if( ValTr[j] != -FLT_MAX && fabs(par1-(
double)ValTr[j]) <= 10.*FLT_EPSILON*fabs(ValTr[j]) )
2200 for(
long j=0; j < grid->
nTracks-1; j++ )
2202 if( ValTr[j] != -FLT_MAX )
2206 for(
long i = j+1; i < grid->
nTracks; i++ )
2208 if( ValTr[i] != -FLT_MAX )
2216 if( j2 > 0 && ((
realnum)par1-ValTr[j])*((
realnum)par1-ValTr[j2]) < 0.f )
2235 for(
long n=0; n < grid->
nTracks; n++ )
2240 for(
long n = 0; n < maxlen; n++ )
2244 for(
long n = 0; n < maxlen; n++ )
2248 vector<long> index(2);
2249 for( index[1]=0; index[1] < grid->
nTracks; ++index[1] )
2252 double Teff, alogg, Mass;
2260 for( index[0]=0; index[0] < grid->
trackLen[index[1]]; ++index[0] )
2263 Teff = grid->
telg[ptr].par[0];
2264 alogg = grid->
telg[ptr].par[1];
2273 const char chSuff[],
2274 const vector<mpp>& telg,
2278 const double par2[],
2284 vector<double> wavl(
NRAUCH);
2285 vector<double> fluxes(
NRAUCH);
2302 long ndim = ( ngrids == 1 ) ? 2 : 3;
2303 vector<string> names;
2304 names.push_back(
"Teff" );
2305 names.push_back(
"log(g)" );
2307 names.push_back(
"log(Z)" );
2308 else if( ngrids == 11 )
2309 names.push_back(
"f(He)" );
2310 else if( ngrids != 1 )
2313 vector<mpp> mytelg(ngrids*telg.size());
2316 for(
long j=0; j < ngrids; j++ )
2317 for(
size_t i=0; i < telg.size(); i++ )
2319 mytelg[k] = telg[i];
2321 mytelg[k].par[2] = par2[j];
2327 "F_lambda", PI*1.e-8, mytelg,
" %.1f", 4);
2330 for(
long i=0; i < nmods; i++ )
2335 sprintf( chLine,
"%7.7ld_%2ld", (
long)(telg[i].par[0]+0.5), (
long)(10.*telg[i].par[1]+0.5) );
2336 else if( format == 2 )
2337 sprintf( chLine,
"%7.7ld_%.2f", (
long)(telg[i].par[0]+0.5), telg[i].par[1] );
2340 string chFileName( chLine );
2341 chFileName += chSuff;
2357 fprintf(
ioQQQ,
" RauchInitialize error in atmosphere file %ld %ld\n", i, j );
2362 while( chLine[0] ==
'*' )
2366 fprintf(
ioQQQ,
" RauchInitialize error in atmosphere file %ld %ld\n", i, j );
2372 for( j=0; j <
NRAUCH; j++ )
2381 fprintf(
ioQQQ,
" RauchInitialize error in atmosphere file %ld %ld\n", i, j );
2387 if( sscanf( chLine,
"%lf %le", &wl, &ttemp ) != 2 )
2389 fprintf(
ioQQQ,
" RauchInitialize error in atmosphere file %ld %ld\n", i, j );
2400 fprintf(
ioQQQ,
" RauchInitialize error in atmosphere file %ld %ld\n", i, j );
2411 if( i == 0 && n == 1 )
2430 const char fnam[] =
"rauch_models.dat";
2437 istringstream iss( line );
2441 fprintf(
ioQQQ,
" RauchReadMPP: the version of %s is not the current version.\n", fnam );
2442 fprintf(
ioQQQ,
" Please obtain the current version from the Cloudy web site.\n" );
2443 fprintf(
ioQQQ,
" I expected to find version %ld and got %ld instead.\n",
2449 unsigned long ndata;
2450 istringstream iss2( line );
2452 ASSERT( ndata == telg1.size() );
2455 getline( ioDATA, line );
2457 for(
unsigned long i=0; i < ndata; ++i )
2458 ioDATA >> telg1[i].par[0] >> telg1[i].par[1];
2459 getline( ioDATA, line );
2462 istringstream iss3( line );
2464 ASSERT( ndata == telg2.size() );
2465 getline( ioDATA, line );
2467 for(
unsigned long i=0; i < ndata; ++i )
2468 ioDATA >> telg2[i].par[0] >> telg2[i].par[1];
2469 getline( ioDATA, line );
2472 istringstream iss4( line );
2474 ASSERT( ndata == telg3.size() );
2475 getline( ioDATA, line );
2477 for(
unsigned long i=0; i < ndata; ++i )
2478 ioDATA >> telg3[i].par[0] >> telg3[i].par[1];
2479 getline( ioDATA, line );
2482 istringstream iss5( line );
2484 ASSERT( ndata == telg4.size() );
2485 getline( ioDATA, line );
2487 for(
unsigned long i=0; i < ndata; ++i )
2488 ioDATA >> telg4[i].par[0] >> telg4[i].par[1];
2489 getline( ioDATA, line );
2492 istringstream iss6( line );
2494 ASSERT( ndata == telg5.size() );
2495 getline( ioDATA, line );
2497 for(
unsigned long i=0; i < ndata; ++i )
2498 ioDATA >> telg5[i].par[0] >> telg5[i].par[1];
2499 getline( ioDATA, line );
2502 istringstream iss7( line );
2504 ASSERT( ndata == telg6.size() );
2505 getline( ioDATA, line );
2507 for(
unsigned long i=0; i < ndata; ++i )
2508 ioDATA >> telg6[i].par[0] >> telg6[i].par[1];
2509 getline( ioDATA, line );
2512 istringstream iss8( line );
2522 getline( ioDATA, line );
2524 while( line[0] ==
'#' );
2530 const vector<string>& names,
2533 const string& wtype,
2535 const string& ftype,
2537 const vector<mpp>& telg,
2543 long npar = (long)names.size();
2546 ASSERT( wtype ==
"lambda" || wtype ==
"nu" );
2547 ASSERT( ftype ==
"F_lambda" || ftype ==
"F_nu" ||
2548 ftype ==
"H_lambda" || ftype ==
"H_nu" );
2551 fprintf( ioOut,
" %ld\n", version );
2552 fprintf( ioOut,
" %ld\n", ndim );
2553 fprintf( ioOut,
" %ld\n", npar );
2554 for(
long i=0; i < npar; ++i )
2555 fprintf( ioOut,
" %s\n", names[i].c_str() );
2556 fprintf( ioOut,
" %ld\n", nmods );
2557 fprintf( ioOut,
" %ld\n", ngrid );
2558 fprintf( ioOut,
" %s\n", wtype.c_str() );
2559 fprintf( ioOut,
" %.8e\n", wfac );
2560 fprintf( ioOut,
" %s\n", ftype.c_str() );
2561 fprintf( ioOut,
" %.8e\n", ffac );
2564 for( i=0; i < nmods; i++ )
2567 for(
long j=0; j < npar; ++j )
2568 fprintf( ioOut, format, telg[i].par[j] );
2570 fprintf( ioOut,
" %c%d", telg[i].chGrid, telg[i].modid );
2571 if( ((i+1)%nmult) == 0 )
2574 if( (i%nmult) != 0 )
2579 const vector<double>& data,
2587 for( i=0; i < ngrid; i++ )
2589 fprintf( ioOut, format, data[i] );
2590 if( ((i+1)%nmult) == 0 )
2593 if( (i%nmult) != 0 )
2607 fprintf(
ioQQQ,
" ReadAtmosphere fails reading line.\n" );
2611 while( chLine[0] ==
'#' );
2615 if( sscanf( chLine,
"%ld", &version ) != 1 )
2617 fprintf(
ioQQQ,
" ReadAtmosphere failed reading VERSION.\n" );
2623 fprintf(
ioQQQ,
" ReadAtmosphere: there is a version number mismatch in"
2624 " the ascii atmosphere file: %s.\n", grid->
name.c_str() );
2625 fprintf(
ioQQQ,
" ReadAtmosphere: Please recreate this file or download the"
2626 " latest version following the instructions on the Cloudy website.\n" );
2633 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid dimension of grid.\n" );
2638 if( fscanf( grid->
ioIN,
"%d", &grid->
npar ) != 1 || grid->
npar <= 0 ||
2641 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid no. of model parameters.\n" );
2645 for(
long nd=0; nd < grid->
npar; nd++ )
2647 if( fscanf( grid->
ioIN,
"%6s", grid->
names[nd] ) != 1 )
2649 fprintf(
ioQQQ,
" ReadAtmosphere failed reading parameter label.\n" );
2655 if( fscanf( grid->
ioIN,
"%d", &grid->
nmods ) != 1 || grid->
nmods <= 0 )
2657 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid number of models.\n" );
2661 if( fscanf( grid->
ioIN,
"%d", &grid->
ngrid ) != 1 || grid->
ngrid <= 1 )
2663 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid number of grid points.\n" );
2667 char chDataType[11];
2669 if( fscanf( grid->
ioIN,
"%10s", chDataType ) != 1 )
2671 fprintf(
ioQQQ,
" ReadAtmosphere failed reading wavl DataType string.\n" );
2675 if( strcmp( chDataType,
"lambda" ) == 0 )
2677 else if( strcmp( chDataType,
"nu" ) == 0 )
2680 fprintf(
ioQQQ,
" ReadAtmosphere found illegal wavl DataType: %s.\n", chDataType );
2686 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid wavl conversion factor.\n" );
2691 if( fscanf( grid->
ioIN,
"%10s", chDataType ) != 1 )
2693 fprintf(
ioQQQ,
" ReadAtmosphere failed reading flux DataType string.\n" );
2697 if( strcmp( chDataType,
"F_lambda" ) == 0 || strcmp( chDataType,
"H_lambda" ) == 0 )
2699 else if( strcmp( chDataType,
"F_nu" ) == 0 || strcmp( chDataType,
"H_nu" ) == 0 )
2702 fprintf(
ioQQQ,
" ReadAtmosphere found illegal flux DataType: %s.\n", chDataType );
2708 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid flux conversion factor.\n" );
2714 for(
long i=0; i < grid->
nmods; i++ )
2716 for(
long nd=0; nd < grid->
npar; nd++ )
2718 if( fscanf( grid->
ioIN,
"%le", &grid->
telg[i].par[nd] ) != 1 )
2720 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid model parameter.\n" );
2726 if( fscanf( grid->
ioIN,
" %c%d", &grid->
telg[i].chGrid, &grid->
telg[i].modid ) != 2 )
2728 fprintf(
ioQQQ,
" ReadAtmosphere failed reading valid track ID.\n" );
2739 const vector<long>& index_list)
2744 vector<realnum> StarEner(grid->
ngrid);
2745 vector<realnum> scratch(grid->
ngrid);
2746 vector<realnum> StarFlux(grid->
ngrid);
2749 for(
long i=0; i < grid->
ngrid; i++ )
2752 if( fscanf( grid->
ioIN,
"%lg", &help ) != 1 )
2754 fprintf(
ioQQQ,
" ReadAtmosphere failed reading wavelength.\n" );
2761 if( scratch[i] <= 0.f )
2763 fprintf(
ioQQQ,
" PROBLEM: a non-positive %s was found, value: %e\n",
2764 grid->
lgFreqX ?
"frequency" :
"wavelength", help );
2769 bool lgFlip = ( !grid->
lgFreqX && scratch[0] < scratch[1] ) || ( grid->
lgFreqX && scratch[0] > scratch[1] );
2772 for(
long i=0; i < grid->
ngrid; i++ )
2776 scratch[i] /= (
realnum)FR1RYD;
2778 scratch[i] = (
realnum)(RYDLAM/scratch[i]);
2781 StarEner[grid->
ngrid-i-1] = scratch[i];
2783 StarEner[i] = scratch[i];
2786 ASSERT( StarEner[0] > 0.f );
2788 for(
long i=1; i < grid->
ngrid; i++ )
2790 if( StarEner[i] <= StarEner[i-1] )
2792 fprintf(
ioQQQ,
" PROBLEM: the %s grid is not strictly monotonically increasing/decreasing\n",
2793 grid->
lgFreqX ?
"frequency" :
"wavelength" );
2799 size_t ni_end = ( index_list.size() == 0 ) ? grid->
nmods : index_list.size();
2801 for(
long imod=0; imod < grid->
nmods; imod++ )
2806 for(
long i=0; i < grid->
ngrid; i++ )
2809 if( fscanf( grid->
ioIN,
"%lg", &help ) != 1 )
2811 fprintf(
ioQQQ,
" ReadAtmosphere failed reading star flux.\n" );
2819 if( scratch[i] < 0.f )
2821 fprintf(
ioQQQ,
"\n PROBLEM: a negative flux was found, model number %ld, value: %e\n",
2827 for(
long i=0; i < grid->
ngrid; i++ )
2830 StarFlux[grid->
ngrid-i-1] = scratch[i];
2832 StarFlux[i] = scratch[i];
2835 for(
long i=0; i < grid->
ngrid; i++ )
2839 StarFlux[i] *= CONVERT_FNU/
POW2(StarEner[i]);
2840 ASSERT( StarFlux[i] >= 0.f );
2849 if( index_list.size() == 0 || imod == index_list[ni] )
2871 const char chFNameOut[],
2880 grid.
name = chFNameIn;
2889 fprintf(
ioQQQ,
" lgCompileAtmosphere got %s.\n", chFNameIn );
2910 val[1] = (int32)
MDIM;
2911 val[2] = (int32)
MNAM;
2912 val[3] = (int32)grid.
ndim;
2913 val[4] = (int32)grid.
npar;
2914 val[5] = (int32)grid.
nmods;
2916 uval[0] =
sizeof(val) +
sizeof(uval) +
sizeof(dval) +
sizeof(md5sum) +
2929 if( fwrite( val,
sizeof(val), 1, ioOUT ) != 1 ||
2930 fwrite( uval,
sizeof(uval), 1, ioOUT ) != 1 ||
2932 fwrite( dval,
sizeof(dval), 1, ioOUT ) != 1 ||
2934 fwrite( md5sum,
sizeof(md5sum), 1, ioOUT ) != 1 ||
2935 fwrite( grid.
names,
sizeof(grid.
names), 1, ioOUT ) != 1 ||
2939 fwrite(
get_ptr(SaveAnu), (
size_t)uval[1], 1, ioOUT ) != 1 )
2941 fprintf(
ioQQQ,
" lgCompileAtmosphere failed writing header of output file.\n" );
2950 for(
long imod=0; imod < grid.
nmods; imod++ )
2953 if( fwrite( &grid.
CloudyFlux[imod][0], (
size_t)uval[1], 1, ioOUT ) != 1 )
2955 fprintf(
ioQQQ,
" lgCompileAtmosphere failed writing star flux.\n" );
2962 fprintf(
ioQQQ,
" lgCompileAtmosphere completed ok.\n\n" );
2977 const char* mode = lgASCII ?
"r" :
"rb";
2984 const char* compilemsg = lgASCII ?
"" :
" and compiled with the COMPILE STARS command";
2985 fprintf(
ioQQQ,
" Error: stellar atmosphere file not found.\n" );
2986 fprintf(
ioQQQ,
"\n\n If the path is set then it is possible that the stellar"
2987 " atmosphere data files do not exist.\n");
2988 fprintf(
ioQQQ,
" Have the stellar data files been downloaded%s?\n", compilemsg );
2989 fprintf(
ioQQQ,
" If you are simply running the test suite and do not need the"
2990 " stellar continua then you should simply ignore this failure\n");
3003 strcmp( grid->
names[0],
"Teff" ) == 0 &&
3004 strcmp( grid->
names[1],
"log(g)" ) == 0 );
3019 fprintf(
ioQQQ,
" InitGrid failed reading header.\n" );
3029 int32 version, mdim, mnam;
3030 double mesh_elo, mesh_ehi;
3033 if( fread( &version,
sizeof(version), 1, grid->
ioIN ) != 1 ||
3034 fread( &mdim,
sizeof(mdim), 1, grid->
ioIN ) != 1 ||
3035 fread( &mnam,
sizeof(mnam), 1, grid->
ioIN ) != 1 ||
3036 fread( &grid->
ndim,
sizeof(grid->
ndim), 1, grid->
ioIN ) != 1 ||
3037 fread( &grid->
npar,
sizeof(grid->
npar), 1, grid->
ioIN ) != 1 ||
3038 fread( &grid->
nmods,
sizeof(grid->
nmods), 1, grid->
ioIN ) != 1 ||
3039 fread( &grid->
ngrid,
sizeof(grid->
ngrid), 1, grid->
ioIN ) != 1 ||
3042 fread( &mesh_elo,
sizeof(mesh_elo), 1, grid->
ioIN ) != 1 ||
3043 fread( &mesh_ehi,
sizeof(mesh_ehi), 1, grid->
ioIN ) != 1 ||
3045 fread( md5sum,
sizeof(md5sum), 1, grid->
ioIN ) != 1 )
3047 fprintf(
ioQQQ,
" InitGrid failed reading header.\n" );
3054 fprintf(
ioQQQ,
" InitGrid: there is a version mismatch between"
3055 " the compiled atmospheres file I expected and the one I found.\n" );
3056 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3057 " atmospheres file with the command: %s.\n", grid->
command.c_str() );
3063 fprintf(
ioQQQ,
" InitGrid: the compiled atmospheres file is produced"
3064 " with an incompatible version of Cloudy.\n" );
3065 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3066 " atmospheres file with the command: %s.\n", grid->
command.c_str() );
3075 fprintf(
ioQQQ,
" InitGrid: the compiled atmospheres file is produced"
3076 " with an incompatible frequency grid.\n" );
3077 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3078 " atmospheres file with the command: %s.\n", grid->
command.c_str() );
3089 if( fread( &grid->
names,
sizeof(grid->
names), 1, grid->
ioIN ) != 1 )
3091 fprintf(
ioQQQ,
" InitGrid failed reading names array.\n" );
3099 fprintf(
ioQQQ,
" InitGrid failed reading model parameter block.\n" );
3107 int res = fseek( grid->
ioIN, 0, SEEK_END );
3110 long End = ftell( grid->
ioIN );
3112 if( End != Expected )
3114 fprintf(
ioQQQ,
" InitGrid: Problem performing sanity check for size of binary file.\n" );
3115 fprintf(
ioQQQ,
" InitGrid: I expected to find %ld bytes, but actually found %ld bytes.\n",
3117 fprintf(
ioQQQ,
" InitGrid: Please recompile the stellar"
3118 " atmospheres file with the command: %s.\n", grid->
command.c_str() );
3138 grid.
name = binName;
3143 int32 version, mdim, mnam;
3144 double mesh_elo, mesh_ehi, mesh_res_factor;
3146 if( fread( &version,
sizeof(version), 1, grid.
ioIN ) != 1 ||
3147 fread( &mdim,
sizeof(mdim), 1, grid.
ioIN ) != 1 ||
3148 fread( &mnam,
sizeof(mnam), 1, grid.
ioIN ) != 1 ||
3149 fread( &grid.
ndim,
sizeof(grid.
ndim), 1, grid.
ioIN ) != 1 ||
3150 fread( &grid.
npar,
sizeof(grid.
npar), 1, grid.
ioIN ) != 1 ||
3155 fread( &mesh_elo,
sizeof(mesh_elo), 1, grid.
ioIN ) != 1 ||
3156 fread( &mesh_ehi,
sizeof(mesh_ehi), 1, grid.
ioIN ) != 1 ||
3157 fread( &mesh_res_factor,
sizeof(mesh_res_factor), 1, grid.
ioIN ) != 1 ||
3158 fread( md5sum,
sizeof(md5sum), 1, grid.
ioIN ) != 1 )
3177 int res = fseek( grid.
ioIN, 0, SEEK_END );
3180 long End = ftell( grid.
ioIN );
3182 if( End != Expected )
3201 if( !ioIN.is_open() )
3206 while( getline( ioIN, chLine ) )
3208 if( chLine[0] !=
'#' )
3216 if( sscanf( chLine.c_str(),
"%ld", &version ) != 1 ||
3240 vector<long> index(2);
3246 char track = (char)(
'A'+index[1]);
3250 for(
long i=0; i < grid->
nmods; i++ )
3252 if( grid->
telg[i].chGrid == track && grid->
telg[i].modid == index[0]+1 )
3266 grid->
trackLen[index[1]] = index[0];
3281 *ndim = (long)grid->
ndim;
3285 val[*nval] = grid->
val[1][grid->
nval[1]-1];
3288 if( *ndim != (
long)grid->
ndim )
3290 fprintf(
ioQQQ,
" A %ld-dim grid was requested, but a %ld-dim grid was found.\n",
3291 *ndim, (
long)grid->
ndim );
3296 fprintf(
ioQQQ,
" A %ld-dim grid was requested, but only %ld parameters were entered.\n",
3313 vector<long> indlo(grid->
ndim);
3314 vector<long> indhi(grid->
ndim);
3315 vector<long> index(grid->
ndim);
3316 vector<double> aval(grid->
npar);
3344 for(
long nd=0; nd < grid->
ndim; nd++ )
3347 FindIndex( grid->
val, nd, grid->
nval[nd], val[nd], &indlo[nd], &indhi[nd], &lgInvalid );
3351 " Requested parameter %s = %.2f is not within the range %.2f to %.2f\n",
3352 grid->
names[nd], val[nd], grid->
val[nd][0], grid->
val[nd][grid->
nval[nd]-1] );
3358 lgTakeLog, Edges, nEdges );
3363 if( grid->
npar == 1 )
3365 " * #<< FINAL: %6s = %13.2f"
3367 grid->
names[0], aval[0] );
3368 else if( grid->
npar == 2 )
3370 " * #<< FINAL: %6s = %10.2f"
3371 " %6s = %8.5f >>> *\n",
3372 grid->
names[0], aval[0], grid->
names[1], aval[1] );
3373 else if( grid->
npar == 3 )
3375 " * #<< FINAL: %6s = %7.0f"
3376 " %6s = %5.2f %6s = %5.2f >>> *\n",
3377 grid->
names[0], aval[0], grid->
names[1], aval[1],
3378 grid->
names[2], aval[2] );
3379 else if( grid->
npar >= 4 )
3382 " * #<< FINAL: %4s = %7.0f"
3383 " %6s = %4.2f %6s = %5.2f %6s = ",
3384 grid->
names[0], aval[0], grid->
names[1], aval[1],
3401 FILE *ioBUG =
open_data(
"interpolated.txt",
"w" );
3407 if( strcmp( grid->
names[0],
"Teff" ) == 0 )
3414 vector<realnum> dum;
3415 SetLimits( grid, val[0], indlo, indhi, NULL, dum, Tlow, Thigh );
3420 vector<double>& aval,
3421 const vector<long>& indlo,
3422 const vector<long>& indhi,
3423 vector<long>& index,
3425 vector<realnum>& flux1,
3436 for( map<string,int>::const_iterator p=grid->
caution.begin(); p != grid->
caution.end(); ++p )
3446 fprintf(
ioQQQ,
"Failed to read atmosphere models.\n" );
3450 for(
size_t i=0; i < grid->
index_list2.size(); ++i )
3459 vector<double>& aval,
3460 const vector<long>& indlo,
3461 const vector<long>& indhi,
3462 vector<long>& index,
3464 vector<realnum>& flux1,
3473 long ind, n =
JIndex(grid,index);
3475 ind = ( grid->
jlo[n] >= 0 ) ? grid->
jlo[n] : grid->
jhi[n];
3477 ind = ( grid->
jhi[n] >= 0 ) ? grid->
jhi[n] : grid->
jlo[n];
3478 else if( grid->
ndim == 1 )
3486 fprintf(
ioQQQ,
" The requested interpolation could not be completed, sorry.\n" );
3487 fprintf(
ioQQQ,
" No suitable match was found for a model with" );
3488 for(
long i=0; i < grid->
ndim; i++ )
3494 for(
long i=0; i < grid->
npar; i++ )
3495 aval[i] = grid->
telg[ind].par[i];
3501 if( !
fp_equal(grid->
val[i][index[i]],aval[i],10) )
3504 oss <<
" No exact match was found for a model with";
3505 for(
long j=0; j < grid->
ndim; j++ )
3506 oss <<
" " << grid->
names[j] <<
"=" << grid->
val[j][index[j]] <<
" ";
3507 oss <<
"- using model " << ind+1 <<
" instead.";
3519 while( i < grid->index_list2.size() && grid->
index_list2[i] != ind )
3531 # if !defined NDEBUG
3532 const realnum SECURE = 10.f*FLT_EPSILON;
3541 index[nd] = indlo[nd];
3542 int next_stage = ( nd == 1 ) ? stage|
IS_FIRST : stage;
3543 InterpolateModel( grid, val, aval, indlo, indhi, index, nd, flux1, next_stage );
3546 vector<double> aval2(grid->
npar);
3548 index[nd] = indhi[nd];
3549 next_stage = ( nd == 1 ) ? stage|
IS_SECOND : stage;
3550 InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, flux2, next_stage );
3554 double fr1 = (aval2[nd]-val[nd])/(aval2[nd]-aval[nd]);
3560 double fr2 = 1. - fr1;
3562 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3565 fprintf(
ioQQQ,
"interpolation nd=%ld fr1=%g\n", nd, fr1 );
3568 double fc1 = 0.,
fc2 = 0.;
3569 if( nd == 0 && strcmp( grid->
names[nd],
"Teff" ) == 0 )
3581 fc1 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->
val[nd][indlo[nd]])*4. : 0.;
3582 fc2 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->
val[nd][indhi[nd]])*4. : 0.;
3586 flux1[i] = (
realnum)(fr1*(flux1[i]+fc1) + fr2*(flux2[i]+
fc2));
3588 for(
long i=0; i < grid->
npar; i++ )
3589 aval[i] = fr1*aval[i] + fr2*aval2[i];
3596 vector<double>& aval,
3599 vector<long>& index,
3602 vector<realnum>& flux1)
3608 long ind = ( index[1] == 0 ) ? indlo[index[0]] : indhi[index[0]];
3611 while( i < grid->index_list2.size() && grid->
index_list2[i] != ind )
3618 for(
long j=0; j < grid->npar; j++ )
3619 aval[j] = grid->
telg[ind].par[j];
3623 # if !defined NDEBUG
3624 const realnum SECURE = 10.f*FLT_EPSILON;
3635 bool lgSkip = ( nd == 1 ) ? ( indhi[index[0]] == indlo[index[0]] ) :
3636 ( indlo[0] == indlo[1] && indhi[0] == indhi[1] );
3641 vector<double> aval2(grid->
npar);
3646 double fr1 = (aval2[nd+off]-val[nd])/(aval2[nd+off]-aval[nd+off]);
3647 double fr2 = 1. - fr1;
3650 fprintf(
ioQQQ,
"interpolation nd=%ld fr1=%g\n", nd, fr1 );
3652 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3655 flux1[i] = (
realnum)(fr1*flux1[i] + fr2*flux2[i]);
3657 for(
long i=0; i < grid->
npar; i++ )
3658 aval[i] = fr1*aval[i] + fr2*aval2[i];
3670 sort( in.begin(), in.end() );
3673 out.push_back( lastval );
3674 for(
size_t i=1; i < in.size(); ++i )
3675 if( in[i] != lastval )
3678 out.push_back( lastval );
3683 vector<Energy>& ener)
3694 if( fseek( grid->
ioIN, (
long)(grid->
nOffset), SEEK_SET ) != 0 )
3696 fprintf(
ioQQQ,
" Error finding atmosphere frequency bins\n");
3703 fprintf(
ioQQQ,
" Error reading atmosphere frequency bins\n" );
3708 ener[i].set(data[i]);
3725 ASSERT( ind >= 0 && ind <= grid->nmods );
3733 fprintf(
ioQQQ,
" Error seeking atmosphere %ld\n", ind );
3739 fprintf(
ioQQQ,
" Error trying to read atmosphere %ld\n", ind );
3748 if( grid->
npar == 1 )
3750 " * #<< %s model%5ld read. "
3751 " %6s = %13.2f >>> *\n",
3752 grid->
ident.c_str(), ind, grid->
names[0], grid->
telg[ind-1].par[0] );
3753 else if( grid->
npar == 2 )
3755 " * #<< %s model%5ld read. "
3756 " %6s = %10.2f %6s = %8.5f >>> *\n",
3757 grid->
ident.c_str(), ind, grid->
names[0], grid->
telg[ind-1].par[0],
3758 grid->
names[1], grid->
telg[ind-1].par[1] );
3759 else if( grid->
npar == 3 )
3761 " * #<< %s model%5ld read. "
3762 " %6s=%7.0f %6s=%5.2f %6s=%5.2f >>> *\n",
3763 grid->
ident.c_str(), ind, grid->
names[0], grid->
telg[ind-1].par[0],
3764 grid->
names[1], grid->
telg[ind-1].par[1],
3765 grid->
names[2], grid->
telg[ind-1].par[2] );
3766 else if( grid->
npar >= 4 )
3770 " %4s=%5.0f %6s=%4.2f %6s=%5.2f %6s=",
3771 grid->
ident.c_str(), ind, grid->
names[0], grid->
telg[ind-1].par[0],
3772 grid->
names[1], grid->
telg[ind-1].par[1],
3787 volatile double help = flux[i];
3799 const vector<long>& indlo,
3800 const vector<long>& indhi,
3802 const vector<realnum>& ValTr,
3810 const double SECURE = (1. + 20.*(double)FLT_EPSILON);
3813 vector<long> index(
MDIM);
3817 switch( grid->
imode )
3827 for(
long j=0; j < grid->
nTracks; j++ )
3829 if( ValTr[j] != -FLT_MAX )
3833 exp10((
double)ValTr[j]) : ValTr[j];
3834 *loLim =
MIN2(*loLim,temp);
3835 *hiLim =
MAX2(*hiLim,temp);
3841 index[1] = useTr[0];
3843 index[1] = useTr[1];
3845 *loLim =
MAX2(grid->
telg[ptr0].par[3],grid->
telg[ptr1].par[3]);
3848 printf(
"set limit 0: (models %d, %d) %f %f\n",
3849 ptr0+1, ptr1+1, grid->
telg[ptr0].par[3], grid->
telg[ptr1].par[3] );
3851 index[0] = grid->
trackLen[useTr[0]]-1;
3852 index[1] = useTr[0];
3854 index[0] = grid->
trackLen[useTr[1]]-1;
3855 index[1] = useTr[1];
3857 *hiLim =
MIN2(grid->
telg[ptr0].par[3],grid->
telg[ptr1].par[3]);
3860 printf(
"set limit 1: (models %d, %d) %f %f\n",
3861 ptr0+1, ptr1+1, grid->
telg[ptr0].par[3], grid->
telg[ptr1].par[3] );
3865 fprintf(
ioQQQ,
" SetLimits called with insane value for imode: %d.\n", grid->
imode );
3869 ASSERT( fabs(*loLim) < DBL_MAX && fabs(*hiLim) < DBL_MAX );
3872 if( *hiLim <= *loLim )
3874 fprintf(
ioQQQ,
" no room to optimize: lower limit %.4f, upper limit %.4f.\n",
3884 printf(
"set limits: %g %g\n",*loLim,*hiLim);
3895 const vector<long>& indlo,
3896 const vector<long>& indhi,
3897 vector<long>& index,
3908 double loLoc = +DBL_MAX;
3909 double hiLoc = -DBL_MAX;
3911 for( index[0]=0; index[0] < grid->
nval[0]; ++index[0] )
3919 long n =
JIndex(grid,index);
3920 if( grid->
jlo[n] < 0 && grid->
jhi[n] < 0 )
3924 if( grid->
val[0][index[0]] < val )
3927 if( grid->
val[0][index[0]] > val )
3934 if( grid->
val[0][index[0]] <= val )
3938 if( loLoc == DBL_MAX )
3939 loLoc = grid->
val[0][index[0]];
3941 if( grid->
val[0][index[0]] >= val )
3945 hiLoc = grid->
val[0][index[0]];
3950 ASSERT( fabs(loLoc) < DBL_MAX && fabs(hiLoc) < DBL_MAX && loLoc <= hiLoc );
3952 *loLim =
MAX2(*loLim,loLoc);
3953 *hiLim =
MIN2(*hiLim,hiLoc);
3957 index[nd] = indlo[nd];
3958 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
3960 if( indhi[nd] != indlo[nd] )
3962 index[nd] = indhi[nd];
3963 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
3979 for(
long nd=0; nd < grid->
ndim; nd++ )
3981 double pval = grid->
telg[0].par[nd];
3982 grid->
val[nd][0] = pval;
3985 for(
long i=1; i < grid->
nmods; i++ )
3990 pval = grid->
telg[i].par[nd];
4000 for(
long j = grid->
nval[nd]-1; j >= i2; j-- )
4001 grid->
val[nd][j+1] = grid->
val[nd][j];
4002 grid->
val[nd][i2] = pval;
4007 jsize *= grid->
nval[nd];
4011 printf(
"%s[%ld]:", grid->
names[nd], grid->
nval[nd] );
4012 for(
long i=0; i < grid->
nval[nd]; i++ )
4013 printf(
" %g", grid->
val[nd][i] );
4018 vector<long> index(grid->
ndim);
4019 vector<double> val(grid->
ndim);
4021 grid->
jlo.resize(jsize);
4022 grid->
jhi.resize(jsize);
4026 FillJ( grid, index, val, grid->
ndim, lgList );
4033 vector<long>& index,
4034 vector<double>& val,
4044 long n =
JIndex(grid,index);
4046 &grid->
jlo[n], &grid->
jhi[n] );
4050 for( index[nd]=0; index[nd] < grid->
nval[nd]; index[nd]++ )
4052 val[nd] = grid->
val[nd][index[nd]];
4053 FillJ( grid, index, val, nd, lgList );
4057 if( lgList && nd ==
MIN2(grid->
ndim-1,1) )
4060 if( grid->
ndim > 2 )
4063 for(
long n = nd+1; n < grid->
ndim; n++ )
4067 if( grid->
ndim > 1 )
4070 for(
long n = 0; n < grid->
nval[1]; n++ )
4074 for(
long n = 0; n < grid->
nval[1]; n++ )
4083 for( index[0]=0; index[0] < grid->
nval[0]; index[0]++ )
4086 if( grid->
ndim > 1 )
4088 for( index[1]=0; index[1] < grid->
nval[1]; index[1]++ )
4106 const vector<long>& index)
4112 for(
long i=0; i < grid->
ndim; i++ )
4114 ind += index[i]*mul;
4115 mul *= grid->
nval[i];
4121 bool lgIsTeffLoggGrid,
4123 const vector<double>& val,
4144 *index_low = *index_high = -2;
4145 double alogg_low = -DBL_MAX, alogg_high = DBL_MAX;
4146 for(
long i=0; i < nmods; i++ )
4148 bool lgNext =
false;
4150 for(
long nd=0; nd < ndim; nd++ )
4152 if( nd != 1 && !
fp_equal(telg[i].par[nd],val[nd],10) )
4162 if( ndim == 1 ||
fp_equal(telg[i].par[1],val[1],10) )
4168 if( lgIsTeffLoggGrid )
4171 if( telg[i].par[1] < val[1] && telg[i].par[1] > alogg_low )
4174 alogg_low = telg[i].par[1];
4177 if( telg[i].par[1] > val[1] && telg[i].par[1] < alogg_high )
4180 alogg_high = telg[i].par[1];
4211 bool lgOutLo = ( x-xval[nd][0] < -10.*DBL_EPSILON*fabs(xval[nd][0]) );
4212 bool lgOutHi = ( x-xval[nd][NVAL-1] > 10.*DBL_EPSILON*fabs(xval[nd][NVAL-1]) );
4214 if( lgOutLo || lgOutHi )
4220 *ind1 = lgOutLo ? -1 : NVAL-1;
4221 *ind2 = lgOutLo ? 0 : NVAL;
4233 for(
long i=0; i < NVAL; i++ )
4244 for(
long i=0; i < NVAL-1; i++ )
4246 if( xval[nd][i] < x && x < xval[nd][i+1] )
4262 FILE *ioIN =
open_data( chFnam,
"r", scheme );
4277 const vector<Energy>& anu)
4288 fprintf(
ioQQQ,
" ValidateMesh: the compiled atmospheres file is produced"
4289 " with an incompatible version of Cloudy.\n" );
4290 fprintf(
ioQQQ,
" ValidateMesh: Please recompile the stellar"
4291 " atmospheres file with the command: %s.\n", grid->
command.c_str() );
4303 if( strcmp( grid->
names[0],
"Teff" ) != 0 )
4311 for(
long i=0; i < grid->
nmods; i++ )
4314 for(
long nd=0; nd < grid->
npar; nd++ )
4325 const vector<realnum>& flux,
4336 lumi += (anu[k].Ryd() - anu[k-1].Ryd())*(flux[k] + flux[k-1])/2.;
4339 double chk =
powpq(lumi*FR1RYD/STEFAN_BOLTZ,1,4);
4341 bool lgPassed =
true;
4342 if( fabs(Teff - chk) > toler*Teff ) {
4343 fprintf(
ioQQQ,
"\n*** WARNING, Teff discrepancy for this model, expected Teff %.2f, ", Teff);
4344 fprintf(
ioQQQ,
"integration yielded Teff %.2f, delta %.2f%%\n", chk, (chk/Teff-1.)*100. );
4352 const vector<realnum>& StarFlux,
4361 for(
long j=0; j < nCont; j++ )
4363 if( StarFlux[j] == 0.f )
4371 vector<long> EdgeInd;
4372 vector<realnum> EdgeLow, EdgeHigh;
4373 for(
long j=0; j < nEdges; j++ )
4377 if( ind >= 1 && ind+2 < nCont )
4379 EdgeInd.push_back( ind );
4380 EdgeLow.push_back( StarEner[ind] );
4381 EdgeHigh.push_back( StarEner[ind+1] );
4385 EdgeInd.push_back( -1 );
4386 EdgeLow.push_back( 0. );
4387 EdgeHigh.push_back( 0. );
4391 vector<realnum> StarPower(nCont-1);
4393 for(
long j=0; j < nCont-1; j++ )
4396 ASSERT( StarEner[j+1] > StarEner[j] );
4400 double ratio_x = (double)StarEner[j+1]/(
double)StarEner[j];
4401 double ratio_y = (double)StarFlux[j+1]/(
double)StarFlux[j];
4402 StarPower[j] = (
realnum)(log(ratio_y)/log(ratio_x));
4410 bool lgDone =
false;
4411 for(
long k=0; k < nEdges; k++ )
4417 ipLo = EdgeInd[k]-1;
4419 ipLo = EdgeInd[k]+1;
4421 if( fabs(StarPower[ipLo]) < fabs(StarPower[EdgeInd[k]]) )
4423 CloudyFlux[j] = StarFlux[ipLo]*
4424 pow(
rfield.
anu(j)/StarEner[ipLo],(double)StarPower[ipLo]);
4440 const vector<realnum>& StarPower,
4449 double widflx =
MIN2(BinHigh,StarEner[nCont-1])-BinLow;
4452 if( BinLow < StarEner[0] )
4456 retval = (
realnum)(StarFlux[0]*
pow2(anu/StarEner[0]));
4458 else if( BinLow > StarEner[nCont-1] )
4467 long ipLo =
RebinFind(StarEner,nCont,BinLow);
4468 long ipHi =
RebinFind(StarEner,nCont,BinHigh);
4470 ASSERT( ipLo >= 0 && ipLo < nCont-1 && ipHi >= ipLo );
4476 retval = (
realnum)(StarFlux[ipLo]*
pow(anu/StarEner[ipLo],(
double)StarPower[ipLo]));
4490 for(
long i=ipLo; i <=
MIN2(ipHi,nCont-2); i++ )
4492 double v1, val,
x1,
x2;
4493 double pp1 = StarPower[i] + 1.;
4499 v1 = StarFlux[i]*
pow(x1/StarEner[i],(
double)StarPower[i]);
4502 else if( i == ipHi )
4517 if( fabs(pp1) < 0.001 )
4519 double z = log(x2/x1);
4521 val = x1*v1*z*(((zp/24.+1./6.)*zp+1./2.)*zp+1.);
4525 val =
pow(x2/x1,pp1) - 1.;
4531 retval = sum/widflx;
4549 if( val < array[0] )
4551 else if( val > array[nArr-1] )
4562 const vector<mpp>& telg,
4569 FILE *ioBUG =
open_data( fnam, ( imod == 0 ) ?
"w" :
"a" );
4571 fprintf( ioBUG,
"######## MODEL %ld", imod+1 );
4572 for(
long nd=0; nd < npar; nd++ )
4573 fprintf( ioBUG,
" %s %g", names[nd], telg[imod].par[nd] );
4574 fprintf( ioBUG,
" ####################\n" );
4576 for(
long k=0; k < nflux; ++k )
4577 fprintf( ioBUG,
"%e %e\n", anu[k], flux[k] );
STATIC realnum RebinSingleCell(long, const realnum[], const realnum[], const vector< realnum > &, long)
static const long int VERSION_BIN
STATIC void RauchReadMPP(vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &)
static const bool lgTAKELOG
int WernerCompile(process_counter &pc)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
static const int IS_SECOND
bool GridCompile(const char *InName)
STATIC void CoStarListModels(const stellar_grid *)
static const int IS_COLLECT
long RauchInterpolateHydr(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
void DumpAtmosphere(const char *fnam, long, long, char[MDIM][MNAM+1], const vector< mpp > &, long, const T[], const realnum[])
STATIC void WriteASCIIHead(FILE *, long, long, const vector< string > &, long, long, const string &, double, const string &, double, const vector< mpp > &, const char *, int)
STATIC void InterpolateRectGrid(stellar_grid *, const double[], double *, double *, bool=lgTAKELOG, const realnum[]=NULL, long=0L)
map< string, int > caution
NORETURN void TotalInsanity(void)
static const realnum Edges_CoStar[]
static const int NMODS_PG1159
string mesh_md5sum() const
multi_arr< double, 2 > val
long CoStarInterpolate(double val[], long *nval, long *ndim, IntMode imode, bool lgHalo, bool lgList, double *val0_lo, double *val0_hi)
static const int NMODS_HYDR
long TlustyInterpolate(double val[], long *nval, long *ndim, tl_grid tlg, const char *chMetalicity, bool lgList, double *Tlow, double *Thigh)
STATIC void InitGrid(stellar_grid *, bool, bool=lgREAD_BIN)
static const int NMODS_HpHE
void invalidate_array(T *p, size_t size)
STATIC bool lgValidBinFile(const char *, process_counter &, access_scheme)
multi_arr< realnum, 2 > CloudyFlux
long HaardtMadauInterpolate(double val, int version, bool lgQuasar, double *zlow, double *zhigh)
int TlustyCompile(process_counter &pc)
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
STATIC bool lgCompileAtmosphere(const char[], const char[], const realnum[], long, process_counter &)
STATIC void InitGridBin(stellar_grid *)
long RauchInterpolateHCa(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
STATIC long JIndex(const stellar_grid *, const vector< long > &)
STATIC void InitGridCoStar(stellar_grid *)
vector< Energy > tNu[LIMSPC]
STATIC void CheckVal(const stellar_grid *, double[], long *, long *)
STATIC bool lgReadAtmosphereHead(stellar_grid *)
long WMBASICInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC void SetLimitsSub(const stellar_grid *, double, const vector< long > &, const vector< long > &, vector< long > &, long, double *, double *)
vector< realnum > tslop[LIMSPC]
static const long int VERSION_RAUCH_MPP
static const long int VERSION_COSTAR
STATIC void InterpolateGridCoStar(stellar_grid *, const double[], double *, double *)
STATIC void InterpolateModelCoStar(const stellar_grid *, const double[], vector< double > &, const long[], const long[], vector< long > &, long, long, vector< realnum > &)
double anu(size_t i) const
long hunt_bisect(const T x[], long n, T xval)
STATIC void SearchModel(const vector< mpp > &, bool, long, const vector< double > &, long, long *, long *)
static const long int VERSION_ASCII
int RauchCompile(process_counter &pc)
long int nflux_with_check
bool lgCoStarInterpolationCaution
STATIC void FindVCoStar(const stellar_grid *, double, vector< realnum > &, long[])
const double * anuptr() const
bool fp_equal(sys_float x, sys_float y, int n=3)
STATIC void GetModel(const stellar_grid *, long, realnum *, bool, bool)
long RauchInterpolateHpHe(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
const ios_base::openmode mode_r
vector< long > index_list
void swap(count_ptr< T > &a, count_ptr< T > &b)
int AtlasCompile(process_counter &pc)
int Kurucz79Compile(process_counter &pc)
STATIC void RebinAtmosphere(const vector< realnum > &, const vector< realnum > &, long, const realnum[], long, realnum[])
static const bool lgREAD_BIN
bool StarburstCompile(process_counter &pc)
long MihalasInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC bool lgReadAtmosphereTail(stellar_grid *, const realnum[], long, const vector< long > &)
int CoStarCompile(process_counter &pc)
void SortUnique(vector< T > &, vector< T > &)
STATIC bool lgValidASCIIFile(const char *, access_scheme)
const int INPUT_LINE_LENGTH
STATIC void FindHCoStar(const stellar_grid *, long, double, long, vector< realnum > &, vector< long > &, vector< long > &)
STATIC void GetBins(const stellar_grid *, vector< Energy > &)
STATIC void SetLimits(const stellar_grid *, double, const vector< long > &, const vector< long > &, const long[], const vector< realnum > &, double *, double *)
double anumin(size_t i) const
STATIC void InitIndexArrays(stellar_grid *, bool)
void InitGridASCII(stellar_grid *)
STATIC void ValidateMesh(const stellar_grid *, const vector< Energy > &)
void getdataline(fstream &, string &)
int MihalasCompile(process_counter &pc)
static const bool DEBUGPRT
static const int NMODS_HELIUM
long RauchInterpolateCOWD(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC void WriteASCIIData(FILE *, const vector< double > &, long, const char *, int)
STATIC void FindIndex(const multi_arr< double, 2 > &, long, long, double, long *, long *, bool *)
STATIC void FillJ(stellar_grid *, vector< long > &, vector< double > &, long, bool)
long WernerInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
#define DEBUG_ENTRY(funcname)
double powpq(double x, int p, int q)
STATIC void InterpolateModel(stellar_grid *, const double[], vector< double > &, const vector< long > &, const vector< long > &, vector< long > &, long, vector< realnum > &, bool, const realnum[]=NULL, long=0L)
STATIC bool RauchInitialize(const char[], const char[], const vector< mpp > &, long, long, long, const double[], int)
static const unsigned int NMD5
static const bool lgVERBOSE
long RauchInterpolateHelium(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
long GridInterpolate(double val[], long *nval, long *ndim, const char *FileName, bool lgList, double *Tlow, double *Thigh)
STATIC void ValidateGrid(const stellar_grid *, double)
static const bool lgLINEAR
static const int NMODS_HCA
int fprintf(const Output &stream, const char *format,...)
int WMBASICCompile(process_counter &pc)
double pow(double x, int i)
static const int IS_FIRST
STATIC bool lgValidModel(const vector< Energy > &, const vector< realnum > &, double, double)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
long RebinFind(const realnum[], long, realnum)
static const int NMODS_HNI
double anumax(size_t i) const
static const int IS_EXECUTE
vector< long > index_list2
long Kurucz79Interpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
static const bool lgSILENT
static const bool lgREAD_ASCII
long RauchInterpolatePG1159(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
STATIC bool CoStarInitialize(const char[], const char[])
STATIC bool lgFileReadable(const char *, process_counter &, access_scheme)
long RauchInterpolateHNi(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
double getResolutionScaleFactor() const
long AtlasInterpolate(double val[], long *nval, long *ndim, const char *chMetalicity, const char *chODFNew, bool lgList, double *Tlow, double *Thigh)
bool StarburstInitialize(const char chInName[], const char chOutName[], sb_mode mode)