cloudy trunk
Loading...
Searching...
No Matches
stars.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3#include "cddefines.h"
4#include "physconst.h"
5#include "optimize.h"
6#include "continuum.h"
7#include "called.h"
8#include "rfield.h"
9#include "thirdparty.h"
10#include "stars.h"
11/*lint -e785 too few initializers */
12/*lint -e801 use of go to depreciated */
13
15static const int NSB99 = 1250;
17static const int MNTS = 200;
18
20static const int NRAUCH = 19951;
22static const int NMODS_HCA = 66;
24static const int NMODS_HNI = 51;
26static const int NMODS_PG1159 = 71;
28static const int NMODS_HYDR = 100;
30static const int NMODS_HELIUM = 81;
32static const int NMODS_HpHE = 117;
33
34/* set to 1 to turn on debug print statements in these routines */
35#define DEBUGPRT 0
36
37#define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
38#define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
39
40static const bool lgSILENT = false;
41static const bool lgVERBOSE = true;
42
43static const bool lgLINEAR = false;
44static const bool lgTAKELOG = true;
45
49
51typedef struct
52{
53 double par[MDIM];
54 int modid;
55 char chGrid;
56} mpp;
57
65
66/* this is the structure of the binary atmosphere file (VERSION 20100902[01]):
67 *
68 * ============================
69 * * int32 VERSION *
70 * * int32 MDIM *
71 * * int32 MNAM *
72 * * int32 ndim *
73 * * int32 npar *
74 * * int32 nmods *
75 * * int32 ngrid *
76 * * uint32 nOffset *
77 * * uint32 nBlocksize *
78 * * double mesh_elo *
79 * * double mesh_ehi *
80 * * double mesh_res_factor *
81 * * char md5sum[NMD5] *
82 * * char names[MDIM][MNAM+1] *
83 * * mpp telg[nmods] *
84 * * realnum anu[ngrid] *
85 * * realnum mod1[ngrid] *
86 * * ... *
87 * * realnum modn[ngrid] *
88 * ============================
89 *
90 * nOffset == 7*sizeof(int32) + 2*sizeof(uint32) + 3*sizeof(double) +
91 * (NMD5 + MDIM*(MNAM+1))*sizeof(char) + nmods*sizeof(mpp)
92 * nBlocksize == ngrid*size(realnum) */
93
95typedef struct
96{
98 string name;
102 access_scheme scheme;
104 FILE *ioIN;
107 const char *ident;
109 const char *command;
111 IntMode imode;
113 int32 ndim;
115 int32 npar;
117 int32 nmods;
119 int32 ngrid;
121 uint32 nOffset;
126 mpp *telg; /* telg[nmods] */
128 double **val; /* val[ndim][nval[n]] */
130 long *nval; /* nval[ndim] */
138 long *jlo; /* jlo(nval[0],...,nval[ndim-1]) */
139 long *jhi; /* jhi(nval[0],...,nval[ndim-1]) */
141 char names[MDIM][MNAM+1];
143 long *trackLen; /* trackLen[nTracks] */
147 long *jval;
149
150/* internal routines */
151STATIC bool lgCompileAtmosphereCoStar(const char[],const char[],const realnum[],long,process_counter&);
152STATIC void InterpolateGridCoStar(const stellar_grid*,const double[],double*,double*);
153STATIC void FindHCoStar(const stellar_grid*,long,double,long,realnum*,long*,long*);
154STATIC void FindVCoStar(const stellar_grid*,double,realnum*,long[]);
156STATIC int RauchInitializeSub(const char[],const char[],const vector<mpp>&,long,long,
157 long,const double[],int);
158STATIC void RauchReadMPP(vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&);
159inline void getdataline(fstream&,string&);
160STATIC bool lgCompileAtmosphere(const char[],const char[],const realnum[],long,process_counter&);
161STATIC void InitGrid(stellar_grid*,bool);
163STATIC bool lgValidAsciiFile(const char*,access_scheme);
165STATIC void CheckVal(const stellar_grid*,double[],long*,long*);
166STATIC void InterpolateRectGrid(const stellar_grid*,const double[],double*,double*);
168STATIC void InterpolateModel(const stellar_grid*,const double[],double[],const long[],
169 const long[],long[],long,vector<realnum>&,IntStage);
170STATIC void InterpolateModelCoStar(const stellar_grid*,const double[],double[],const long[],
171 const long[],long[],long,long,vector<realnum>&);
172STATIC void GetBins(const stellar_grid*,vector<Energy>&);
173STATIC void GetModel(const stellar_grid*,long,vector<realnum>&,bool,bool);
174STATIC void SetLimits(const stellar_grid*,double,const long[],const long[],const long[],
175 const realnum[],double*,double*);
176STATIC void SetLimitsSub(const stellar_grid*,double,const long[],const long[],long[],long,
177 double*,double*);
179STATIC void FillJ(const stellar_grid*,long[],double[],long,bool);
180STATIC long JIndex(const stellar_grid*,const long[]);
181STATIC void SearchModel(const mpp[],bool,long,const double[],long,long*,long*);
182STATIC void FindIndex(const double[],long,double,long*,long*,bool*);
184STATIC void ValidateGrid(const stellar_grid*,double);
185STATIC bool lgValidModel(const vector<Energy>&,const vector<realnum>&,double,double);
186STATIC void RebinAtmosphere(long,const realnum[],const realnum[],realnum[],long,const realnum[]);
187STATIC realnum RebinSingleCell(realnum,realnum,const realnum[],const realnum[],const realnum[],long);
188STATIC long RebinFind(const realnum[],long,realnum);
189
190
191/* the version number for the ascii/binary atmosphere files */
192static const long int VERSION_ASCII = 20060612L;
193/* binary files are incompatible when floats are converted to doubles */
194#ifdef FLT_IS_DBL
195static const long int VERSION_BIN = 201009020L;
196#else
197static const long int VERSION_BIN = 201009021L;
198#endif
199static const long int VERSION_RAUCH_MPP = 20090324;
200
203{
204 DEBUG_ENTRY( "AtmospheresAvail()" );
205
206 /* This routine makes a list of all the stellar atmosphere grids that are valid,
207 * giving the parameters for use in the input script as well. It is simply a long
208 * list of if-statements, so if any grid is added to Cloudy, it should be added in
209 * this routine as well.
210 *
211 * NB NB NB -- test this routine regularly to see if the list is still complete! */
212
213 fprintf( ioQQQ, "\n I will now list all stellar atmosphere grids that are ready to be used (if any).\n" );
214 fprintf( ioQQQ, " User-defined stellar atmosphere grids will not be included in this list.\n\n" );
215
216 process_counter dum;
217
218 /* we always look in the data directory regardless of where we are,
219 * it would be very confusing to the user if we did otherwise... */
221
222 if( lgValidBinFile( "atlas_fp10k2.mod", dum, as ) )
223 fprintf( ioQQQ, " table star atlas Z+1.0 <Teff> [ <log(g)> ]\n" );
224 if( lgValidBinFile( "atlas_fp05k2.mod", dum, as ) )
225 fprintf( ioQQQ, " table star atlas Z+0.5 <Teff> [ <log(g)> ]\n" );
226 if( lgValidBinFile( "atlas_fp03k2.mod", dum, as ) )
227 fprintf( ioQQQ, " table star atlas Z+0.3 <Teff> [ <log(g)> ]\n" );
228 if( lgValidBinFile( "atlas_fp02k2.mod", dum, as ) )
229 fprintf( ioQQQ, " table star atlas Z+0.2 <Teff> [ <log(g)> ]\n" );
230 if( lgValidBinFile( "atlas_fp01k2.mod", dum, as ) )
231 fprintf( ioQQQ, " table star atlas Z+0.1 <Teff> [ <log(g)> ]\n" );
232 if( lgValidBinFile( "atlas_fp00k2.mod", dum, as ) )
233 fprintf( ioQQQ, " table star atlas Z+0.0 <Teff> [ <log(g)> ]\n" );
234 if( lgValidBinFile( "atlas_fm01k2.mod", dum, as ) )
235 fprintf( ioQQQ, " table star atlas Z-0.1 <Teff> [ <log(g)> ]\n" );
236 if( lgValidBinFile( "atlas_fm02k2.mod", dum, as ) )
237 fprintf( ioQQQ, " table star atlas Z-0.2 <Teff> [ <log(g)> ]\n" );
238 if( lgValidBinFile( "atlas_fm03k2.mod", dum, as ) )
239 fprintf( ioQQQ, " table star atlas Z-0.3 <Teff> [ <log(g)> ]\n" );
240 if( lgValidBinFile( "atlas_fm05k2.mod", dum, as ) )
241 fprintf( ioQQQ, " table star atlas Z-0.5 <Teff> [ <log(g)> ]\n" );
242 if( lgValidBinFile( "atlas_fm10k2.mod", dum, as ) )
243 fprintf( ioQQQ, " table star atlas Z-1.0 <Teff> [ <log(g)> ]\n" );
244 if( lgValidBinFile( "atlas_fm15k2.mod", dum, as ) )
245 fprintf( ioQQQ, " table star atlas Z-1.5 <Teff> [ <log(g)> ]\n" );
246 if( lgValidBinFile( "atlas_fm20k2.mod", dum, as ) )
247 fprintf( ioQQQ, " table star atlas Z-2.0 <Teff> [ <log(g)> ]\n" );
248 if( lgValidBinFile( "atlas_fm25k2.mod", dum, as ) )
249 fprintf( ioQQQ, " table star atlas Z-2.5 <Teff> [ <log(g)> ]\n" );
250 if( lgValidBinFile( "atlas_fm30k2.mod", dum, as ) )
251 fprintf( ioQQQ, " table star atlas Z-3.0 <Teff> [ <log(g)> ]\n" );
252 if( lgValidBinFile( "atlas_fm35k2.mod", dum, as ) )
253 fprintf( ioQQQ, " table star atlas Z-3.5 <Teff> [ <log(g)> ]\n" );
254 if( lgValidBinFile( "atlas_fm40k2.mod", dum, as ) )
255 fprintf( ioQQQ, " table star atlas Z-4.0 <Teff> [ <log(g)> ]\n" );
256 if( lgValidBinFile( "atlas_fm45k2.mod", dum, as ) )
257 fprintf( ioQQQ, " table star atlas Z-4.5 <Teff> [ <log(g)> ]\n" );
258 if( lgValidBinFile( "atlas_fm50k2.mod", dum, as ) )
259 fprintf( ioQQQ, " table star atlas Z-5.0 <Teff> [ <log(g)> ]\n" );
260
261 if( lgValidBinFile( "atlas_fp05k2_odfnew.mod", dum, as ) )
262 fprintf( ioQQQ, " table star atlas odfnew Z+0.5 <Teff> [ <log(g)> ]\n" );
263 if( lgValidBinFile( "atlas_fp02k2_odfnew.mod", dum, as ) )
264 fprintf( ioQQQ, " table star atlas odfnew Z+0.2 <Teff> [ <log(g)> ]\n" );
265 if( lgValidBinFile( "atlas_fp00k2_odfnew.mod", dum, as ) )
266 fprintf( ioQQQ, " table star atlas odfnew Z+0.0 <Teff> [ <log(g)> ]\n" );
267 if( lgValidBinFile( "atlas_fm05k2_odfnew.mod", dum, as ) )
268 fprintf( ioQQQ, " table star atlas odfnew Z-0.5 <Teff> [ <log(g)> ]\n" );
269 if( lgValidBinFile( "atlas_fm10k2_odfnew.mod", dum, as ) )
270 fprintf( ioQQQ, " table star atlas odfnew Z-1.0 <Teff> [ <log(g)> ]\n" );
271 if( lgValidBinFile( "atlas_fm15k2_odfnew.mod", dum, as ) )
272 fprintf( ioQQQ, " table star atlas odfnew Z-1.5 <Teff> [ <log(g)> ]\n" );
273 if( lgValidBinFile( "atlas_fm20k2_odfnew.mod", dum, as ) )
274 fprintf( ioQQQ, " table star atlas odfnew Z-2.0 <Teff> [ <log(g)> ]\n" );
275 if( lgValidBinFile( "atlas_fm25k2_odfnew.mod", dum, as ) )
276 fprintf( ioQQQ, " table star atlas odfnew Z-2.5 <Teff> [ <log(g)> ]\n" );
277
278 if( lgValidBinFile( "atlas_3d.mod", dum, as ) )
279 fprintf( ioQQQ, " table star atlas 3-dim <Teff> <log(g)> <log(Z)>\n" );
280
281 if( lgValidBinFile( "atlas_3d_odfnew.mod", dum, as ) )
282 fprintf( ioQQQ, " table star atlas odfnew 3-dim <Teff> <log(g)> <log(Z)>\n" );
283
284 if( lgValidBinFile( "Sc1_costar_solar.mod", dum, as ) )
285 fprintf( ioQQQ, " table star costar solar (see Hazy for parameters)\n" );
286 if( lgValidBinFile( "Sc1_costar_halo.mod", dum, as ) )
287 fprintf( ioQQQ, " table star costar halo (see Hazy for parameters)\n" );
288
289 if( lgValidBinFile( "kurucz79.mod", dum, as ) )
290 fprintf( ioQQQ, " table star kurucz79 <Teff>\n" );
291
292 if( lgValidBinFile( "mihalas.mod", dum, as ) )
293 fprintf( ioQQQ, " table star mihalas <Teff>\n" );
294
295 if( lgValidBinFile( "rauch_h-ca_solar.mod", dum, as ) )
296 fprintf( ioQQQ, " table star rauch H-Ca solar <Teff> [ <log(g)> ]\n" );
297 if( lgValidBinFile( "rauch_h-ca_halo.mod", dum, as ) )
298 fprintf( ioQQQ, " table star rauch H-Ca halo <Teff> [ <log(g)> ]\n" );
299 if( lgValidBinFile( "rauch_h-ca_3d.mod", dum, as ) )
300 fprintf( ioQQQ, " table star rauch H-Ca 3-dim <Teff> <log(g)> <log(Z)>\n" );
301
302 if( lgValidBinFile( "rauch_h-ni_solar.mod", dum, as ) )
303 fprintf( ioQQQ, " table star rauch H-Ni solar <Teff> [ <log(g)> ]\n" );
304 if( lgValidBinFile( "rauch_h-ni_halo.mod", dum, as ) )
305 fprintf( ioQQQ, " table star rauch H-Ni halo <Teff> [ <log(g)> ]\n" );
306 if( lgValidBinFile( "rauch_h-ni_3d.mod", dum, as ) )
307 fprintf( ioQQQ, " table star rauch H-Ni 3-dim <Teff> <log(g)> <log(Z)>\n" );
308
309 if( lgValidBinFile( "rauch_pg1159.mod", dum, as ) )
310 fprintf( ioQQQ, " table star rauch pg1159 <Teff> [ <log(g)> ]\n" );
311 if( lgValidBinFile( "rauch_cowd.mod", dum, as ) )
312 fprintf( ioQQQ, " table star rauch co wd <Teff>\n" );
313
314 if( lgValidBinFile( "rauch_hydr.mod", dum, as ) )
315 fprintf( ioQQQ, " table star rauch hydrogen <Teff> [ <log(g)> ]\n" );
316
317 if( lgValidBinFile( "rauch_helium.mod", dum, as ) )
318 fprintf( ioQQQ, " table star rauch helium <Teff> [ <log(g)> ]\n" );
319
320 if( lgValidBinFile( "rauch_h+he_3d.mod", dum, as ) )
321 fprintf( ioQQQ, " table star rauch H+He <Teff> <log(g)> <frac(He)>\n" );
322
323 if( lgValidBinFile( "starburst99.mod", dum, as ) )
324 fprintf( ioQQQ, " table star \"starburst99.mod\" <age>\n" );
325 if( lgValidBinFile( "starburst99_2d.mod", dum, as ) )
326 fprintf( ioQQQ, " table star \"starburst99_2d.mod\" <age> <Z>\n" );
327
328 if( lgValidBinFile( "obstar_merged_p03.mod", dum, as ) )
329 fprintf( ioQQQ, " table star tlusty OBstar Z+0.3 <Teff> [ <log(g)> ]\n" );
330 if( lgValidBinFile( "obstar_merged_p00.mod", dum, as ) )
331 fprintf( ioQQQ, " table star tlusty OBstar Z+0.0 <Teff> [ <log(g)> ]\n" );
332 if( lgValidBinFile( "obstar_merged_m03.mod", dum, as ) )
333 fprintf( ioQQQ, " table star tlusty OBstar Z-0.3 <Teff> [ <log(g)> ]\n" );
334 if( lgValidBinFile( "obstar_merged_m07.mod", dum, as ) )
335 fprintf( ioQQQ, " table star tlusty OBstar Z-0.7 <Teff> [ <log(g)> ]\n" );
336 if( lgValidBinFile( "obstar_merged_m10.mod", dum, as ) )
337 fprintf( ioQQQ, " table star tlusty OBstar Z-1.0 <Teff> [ <log(g)> ]\n" );
338 if( lgValidBinFile( "obstar_merged_m99.mod", dum, as ) )
339 fprintf( ioQQQ, " table star tlusty OBstar Z-inf <Teff> [ <log(g)> ]\n" );
340
341 if( lgValidBinFile( "obstar_merged_3d.mod", dum, as ) )
342 fprintf( ioQQQ, " table star tlusty OBstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
343
344 if( lgValidBinFile( "bstar2006_p03.mod", dum, as ) )
345 fprintf( ioQQQ, " table star tlusty Bstar Z+0.3 <Teff> [ <log(g)> ]\n" );
346 if( lgValidBinFile( "bstar2006_p00.mod", dum, as ) )
347 fprintf( ioQQQ, " table star tlusty Bstar Z+0.0 <Teff> [ <log(g)> ]\n" );
348 if( lgValidBinFile( "bstar2006_m03.mod", dum, as ) )
349 fprintf( ioQQQ, " table star tlusty Bstar Z-0.3 <Teff> [ <log(g)> ]\n" );
350 if( lgValidBinFile( "bstar2006_m07.mod", dum, as ) )
351 fprintf( ioQQQ, " table star tlusty Bstar Z-0.7 <Teff> [ <log(g)> ]\n" );
352 if( lgValidBinFile( "bstar2006_m10.mod", dum, as ) )
353 fprintf( ioQQQ, " table star tlusty Bstar Z-1.0 <Teff> [ <log(g)> ]\n" );
354 if( lgValidBinFile( "bstar2006_m99.mod", dum, as ) )
355 fprintf( ioQQQ, " table star tlusty Bstar Z-inf <Teff> [ <log(g)> ]\n" );
356
357 if( lgValidBinFile( "bstar2006_3d.mod", dum, as ) )
358 fprintf( ioQQQ, " table star tlusty Bstar 3-dim <Teff> <log(g)> <log(Z)>\n" );
359
360 if( lgValidBinFile( "ostar2002_p03.mod", dum, as ) )
361 fprintf( ioQQQ, " table star tlusty Ostar Z+0.3 <Teff> [ <log(g)> ]\n" );
362 if( lgValidBinFile( "ostar2002_p00.mod", dum, as ) )
363 fprintf( ioQQQ, " table star tlusty Ostar Z+0.0 <Teff> [ <log(g)> ]\n" );
364 if( lgValidBinFile( "ostar2002_m03.mod", dum, as ) )
365 fprintf( ioQQQ, " table star tlusty Ostar Z-0.3 <Teff> [ <log(g)> ]\n" );
366 if( lgValidBinFile( "ostar2002_m07.mod", dum, as ) )
367 fprintf( ioQQQ, " table star tlusty Ostar Z-0.7 <Teff> [ <log(g)> ]\n" );
368 if( lgValidBinFile( "ostar2002_m10.mod", dum, as ) )
369 fprintf( ioQQQ, " table star tlusty Ostar Z-1.0 <Teff> [ <log(g)> ]\n" );
370 if( lgValidBinFile( "ostar2002_m15.mod", dum, as ) )
371 fprintf( ioQQQ, " table star tlusty Ostar Z-1.5 <Teff> [ <log(g)> ]\n" );
372 if( lgValidBinFile( "ostar2002_m17.mod", dum, as ) )
373 fprintf( ioQQQ, " table star tlusty Ostar Z-1.7 <Teff> [ <log(g)> ]\n" );
374 if( lgValidBinFile( "ostar2002_m20.mod", dum, as ) )
375 fprintf( ioQQQ, " table star tlusty Ostar Z-2.0 <Teff> [ <log(g)> ]\n" );
376 if( lgValidBinFile( "ostar2002_m30.mod", dum, as ) )
377 fprintf( ioQQQ, " table star tlusty Ostar Z-3.0 <Teff> [ <log(g)> ]\n" );
378 if( lgValidBinFile( "ostar2002_m99.mod", dum, as ) )
379 fprintf( ioQQQ, " table star tlusty Ostar Z-inf <Teff> [ <log(g)> ]\n" );
380
381 if( lgValidBinFile( "ostar2002_3d.mod", dum, as ) )
382 fprintf( ioQQQ, " table star tlusty Ostar 3-dim <Teff> <log(g)> <log(Z)>\n" );
383
384 if( lgValidBinFile( "kwerner.mod", dum, as ) )
385 fprintf( ioQQQ, " table star werner <Teff> [ <log(g)> ]\n" );
386
387 if( lgValidBinFile( "wmbasic.mod", dum, as ) )
388 fprintf( ioQQQ, " table star wmbasic <Teff> <log(g)> <log(Z)>\n" );
389 return;
390}
391
392/* AtlasCompile rebin Kurucz stellar models to match energy grid of code */
393/* >>chng 05 nov 16, added return value to indicate success (0) or failure (1) */
395{
396 /* these contain frequencies for the major absorption edges */
397 realnum Edges[3];
398
399 bool lgFail = false;
400
401 DEBUG_ENTRY( "AtlasCompile()" );
402
403 /* This is a program to re-bin the Kurucz stellar models spectrum to match the
404 * CLOUDY grid. For wavelengths shorter than supplied in the Kurucz files,
405 * the flux will be set to zero. At long wavelengths a Rayleigh-Jeans
406 * extrapolation will be used. */
407
408 /* This version uses power-law interpolation between the points of the stellar
409 * model.*/
410
411 fprintf( ioQQQ, " AtlasCompile on the job.\n" );
412
413 /* define the major absorption edges that require special attention during rebinning
414 *
415 * NB the frequencies should be chosen here such that they are somewhere in between
416 * the two frequency points that straddle the edge in the atmosphere model, the
417 * software in RebinAtmosphere will seek out the exact values of those two points
418 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
419 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
420 *
421 * NB beware not to choose edges too close to one another (i.e. on the order of the
422 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
423 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
424 * almost certainly lead to erroneous behaviour in RebinAtmosphere */
425 Edges[0] = (realnum)(RYDLAM/911.76);
426 Edges[1] = (realnum)(RYDLAM/504.26);
427 Edges[2] = (realnum)(RYDLAM/227.84);
428
430
431 /* >>chng 05 nov 19, add support for non-solar metalicities as well as odfnew models, PvH */
432 if( lgFileReadable( "atlas_fp10k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp10k2.mod", pc, as ) )
433 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp10k2.ascii", "atlas_fp10k2.mod", Edges, 3L, pc );
434 if( lgFileReadable( "atlas_fp05k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp05k2.mod", pc, as ) )
435 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp05k2.ascii", "atlas_fp05k2.mod", Edges, 3L, pc );
436 if( lgFileReadable( "atlas_fp03k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp03k2.mod", pc, as ) )
437 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp03k2.ascii", "atlas_fp03k2.mod", Edges, 3L, pc );
438 if( lgFileReadable( "atlas_fp02k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp02k2.mod", pc, as ) )
439 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp02k2.ascii", "atlas_fp02k2.mod", Edges, 3L, pc );
440 if( lgFileReadable( "atlas_fp01k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp01k2.mod", pc, as ) )
441 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp01k2.ascii", "atlas_fp01k2.mod", Edges, 3L, pc );
442 if( lgFileReadable( "atlas_fp00k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fp00k2.mod", pc, as ) )
443 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp00k2.ascii", "atlas_fp00k2.mod", Edges, 3L, pc );
444 if( lgFileReadable( "atlas_fm01k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm01k2.mod", pc, as ) )
445 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm01k2.ascii", "atlas_fm01k2.mod", Edges, 3L, pc );
446 if( lgFileReadable( "atlas_fm02k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm02k2.mod", pc, as ) )
447 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm02k2.ascii", "atlas_fm02k2.mod", Edges, 3L, pc );
448 if( lgFileReadable( "atlas_fm03k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm03k2.mod", pc, as ) )
449 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm03k2.ascii", "atlas_fm03k2.mod", Edges, 3L, pc );
450 if( lgFileReadable( "atlas_fm05k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm05k2.mod", pc, as ) )
451 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm05k2.ascii", "atlas_fm05k2.mod", Edges, 3L, pc );
452 if( lgFileReadable( "atlas_fm10k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm10k2.mod", pc, as ) )
453 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm10k2.ascii", "atlas_fm10k2.mod", Edges, 3L, pc );
454 if( lgFileReadable( "atlas_fm15k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm15k2.mod", pc, as ) )
455 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm15k2.ascii", "atlas_fm15k2.mod", Edges, 3L, pc );
456 if( lgFileReadable( "atlas_fm20k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm20k2.mod", pc, as ) )
457 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm20k2.ascii", "atlas_fm20k2.mod", Edges, 3L, pc );
458 if( lgFileReadable( "atlas_fm25k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm25k2.mod", pc, as ) )
459 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm25k2.ascii", "atlas_fm25k2.mod", Edges, 3L, pc );
460 if( lgFileReadable( "atlas_fm30k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm30k2.mod", pc, as ) )
461 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm30k2.ascii", "atlas_fm30k2.mod", Edges, 3L, pc );
462 if( lgFileReadable( "atlas_fm35k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm35k2.mod", pc, as ) )
463 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm35k2.ascii", "atlas_fm35k2.mod", Edges, 3L, pc );
464 if( lgFileReadable( "atlas_fm40k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm40k2.mod", pc, as ) )
465 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm40k2.ascii", "atlas_fm40k2.mod", Edges, 3L, pc );
466 if( lgFileReadable( "atlas_fm45k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm45k2.mod", pc, as ) )
467 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm45k2.ascii", "atlas_fm45k2.mod", Edges, 3L, pc );
468 if( lgFileReadable( "atlas_fm50k2.ascii", pc, as ) && !lgValidBinFile( "atlas_fm50k2.mod", pc, as ) )
469 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm50k2.ascii", "atlas_fm50k2.mod", Edges, 3L, pc );
470
471 if( lgFileReadable( "atlas_fp05k2_odfnew.ascii", pc, as ) &&
472 !lgValidBinFile( "atlas_fp05k2_odfnew.mod", pc, as ) )
473
474 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp05k2_odfnew.ascii", "atlas_fp05k2_odfnew.mod",
475 Edges, 3L, pc );
476 if( lgFileReadable( "atlas_fp02k2_odfnew.ascii", pc, as ) &&
477 !lgValidBinFile( "atlas_fp02k2_odfnew.mod", pc, as ) )
478 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp02k2_odfnew.ascii", "atlas_fp02k2_odfnew.mod",
479 Edges, 3L, pc );
480 if( lgFileReadable( "atlas_fp00k2_odfnew.ascii", pc, as ) &&
481 !lgValidBinFile( "atlas_fp00k2_odfnew.mod", pc, as ) )
482 lgFail = lgFail || lgCompileAtmosphere( "atlas_fp00k2_odfnew.ascii", "atlas_fp00k2_odfnew.mod",
483 Edges, 3L, pc );
484 if( lgFileReadable( "atlas_fm05k2_odfnew.ascii", pc, as ) &&
485 !lgValidBinFile( "atlas_fm05k2_odfnew.mod", pc, as ) )
486 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm05k2_odfnew.ascii", "atlas_fm05k2_odfnew.mod",
487 Edges, 3L, pc );
488 if( lgFileReadable( "atlas_fm10k2_odfnew.ascii", pc, as ) &&
489 !lgValidBinFile( "atlas_fm10k2_odfnew.mod", pc, as ) )
490 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm10k2_odfnew.ascii", "atlas_fm10k2_odfnew.mod",
491 Edges, 3L, pc );
492 if( lgFileReadable( "atlas_fm15k2_odfnew.ascii", pc, as ) &&
493 !lgValidBinFile( "atlas_fm15k2_odfnew.mod", pc, as ) )
494 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm15k2_odfnew.ascii", "atlas_fm15k2_odfnew.mod",
495 Edges, 3L, pc );
496 if( lgFileReadable( "atlas_fm20k2_odfnew.ascii", pc, as ) &&
497 !lgValidBinFile( "atlas_fm20k2_odfnew.mod", pc, as ) )
498 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm20k2_odfnew.ascii", "atlas_fm20k2_odfnew.mod",
499 Edges, 3L, pc );
500 if( lgFileReadable( "atlas_fm25k2_odfnew.ascii", pc, as ) &&
501 !lgValidBinFile( "atlas_fm25k2_odfnew.mod", pc, as ) )
502 lgFail = lgFail || lgCompileAtmosphere( "atlas_fm25k2_odfnew.ascii", "atlas_fm25k2_odfnew.mod",
503 Edges, 3L, pc );
504
505 if( lgFileReadable( "atlas_3d.ascii", pc, as ) && !lgValidBinFile( "atlas_3d.mod", pc, as ) )
506 lgFail = lgFail || lgCompileAtmosphere( "atlas_3d.ascii", "atlas_3d.mod", Edges, 3L, pc );
507
508 if( lgFileReadable( "atlas_3d_odfnew.ascii", pc, as ) &&
509 !lgValidBinFile( "atlas_3d_odfnew.mod", pc, as ) )
510 lgFail = lgFail || lgCompileAtmosphere( "atlas_3d_odfnew.ascii", "atlas_3d_odfnew.mod", Edges, 3L, pc );
511 return lgFail;
512}
513
514/* AtlasInterpolate read in and interpolate on Kurucz grid of atmospheres, originally by K Volk */
515long AtlasInterpolate(double val[], /* val[nval] */
516 long *nval,
517 long *ndim,
518 const char *chMetalicity,
519 const char *chODFNew,
520 bool lgList,
521 double *Tlow,
522 double *Thigh)
523{
524 char chIdent[13];
526
527 DEBUG_ENTRY( "AtlasInterpolate()" );
528
529 grid.name = "atlas_";
530 if( *ndim == 3 )
531 grid.name += "3d";
532 else
533 {
534 grid.name += "f";
535 grid.name += chMetalicity;
536 grid.name += "k2";
537 }
538 grid.name += chODFNew;
539 grid.name += ".mod";
540 grid.scheme = AS_DATA_OPTIONAL;
541 /* identification of this atmosphere set, used in
542 * the Cloudy output, *must* be 12 characters long */
543 if( *ndim == 3 )
544 {
545 strcpy( chIdent, "3-dim" );
546 }
547 else
548 {
549 strcpy( chIdent, "Z " );
550 strcat( chIdent, chMetalicity );
551 }
552 strcat( chIdent, ( strlen(chODFNew) == 0 ? " Kurucz" : " ODFNew" ) );
553 grid.ident = chIdent;
554 /* the Cloudy command needed to recompile the binary model file */
555 grid.command = "COMPILE STARS";
556
557 InitGrid( &grid, lgList );
558
559 CheckVal( &grid, val, nval, ndim );
560
561 /* Note on the interpolation (solar abundance grid): 26 October 2000 (Peter van Hoof)
562 *
563 * I computed the effective temperature for a random sample of interpolated
564 * atmospheres by integrating the flux as shown above and compared the results
565 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
566 *
567 * I found that the average discrepancy was:
568 *
569 * DELTA = -0.10% +/- 0.06% (sample size 5000)
570 *
571 * The most extreme discrepancies were
572 * -0.30% <= DELTA <= 0.21%
573 *
574 * The most negative discrepancies were for Teff = 36 - 39 kK, log(g) = 4.5 - 5
575 * The most positive discrepancies were for Teff = 3.5 - 4.0 kK, log(g) = 0 - 1
576 *
577 * The interpolation in the ATLAS grid is clearly very accurate */
578
579 InterpolateRectGrid( &grid, val, Tlow, Thigh );
580
581 FreeGrid( &grid );
582 return rfield.nupper;
583}
584
585/* CoStarCompile rebin costar stellar models to match energy grid of code*/
587{
588 realnum Edges[3];
589 bool lgFail = false;
590
591 DEBUG_ENTRY( "CoStarCompile()" );
592
593 fprintf( ioQQQ, " CoStarCompile on the job.\n" );
594
595 /* define the major absorption edges that require special attention during rebinning
596 *
597 * NB the frequencies should be chosen here such that they are somewhere in between
598 * the two frequency points that straddle the edge in the atmosphere model, the
599 * software in RebinAtmosphere will seek out the exact values of those two points
600 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
601 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
602 *
603 * NB beware not to choose edges too close to one another (i.e. on the order of the
604 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
605 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
606 * almost certainly lead to erroneous behaviour in RebinAtmosphere */
607 Edges[0] = (realnum)(RYDLAM/911.76);
608 Edges[1] = (realnum)(RYDLAM/504.26);
609 Edges[2] = (realnum)(RYDLAM/227.84);
610
612
613 if( lgFileReadable( "Sc1_costar_z020_lb.fluxes", pc, as ) && !lgValidBinFile( "Sc1_costar_solar.mod", pc, as ) )
614 lgFail = lgFail || lgCompileAtmosphereCoStar( "Sc1_costar_z020_lb.fluxes", "Sc1_costar_solar.mod",
615 Edges, 3L, pc );
616 if( lgFileReadable( "Sc1_costar_z004_lb.fluxes", pc, as ) && !lgValidBinFile( "Sc1_costar_halo.mod", pc, as ) )
617 lgFail = lgFail || lgCompileAtmosphereCoStar( "Sc1_costar_z004_lb.fluxes", "Sc1_costar_halo.mod",
618 Edges, 3L, pc );
619 return lgFail;
620}
621
622/* CoStarInterpolate read in and interpolate on CoStar grid of atmospheres */
623long CoStarInterpolate(double val[], /* requested model parameters */
624 long *nval,
625 long *ndim,
626 IntMode imode, /* which interpolation mode is requested */
627 bool lgHalo, /* flag indicating whether solar (==0) or halo (==1) abundances */
628 bool lgList,
629 double *val0_lo,
630 double *val0_hi)
631{
633
634 DEBUG_ENTRY( "CoStarInterpolate()" );
635
636 grid.name = ( lgHalo ? "Sc1_costar_halo.mod" : "Sc1_costar_solar.mod" );
637 grid.scheme = AS_DATA_OPTIONAL;
638 /* identification of this atmosphere set, used in
639 * the Cloudy output, *must* be 12 characters long */
640 grid.ident = " costar";
641 /* the Cloudy command needed to recompile the binary model file */
642 grid.command = "COMPILE STARS";
643
644 /* listing the models in the grid is implemented in CoStarListModels() */
645 InitGrid( &grid, false );
646 /* now sort the models according to track */
648 /* override default interpolation mode */
649 grid.imode = imode;
650
651 if( lgList )
652 {
655 }
656
657 CheckVal( &grid, val, nval, ndim );
658
659 /* Note on the interpolation: 26 October 2000 (Peter van Hoof)
660 *
661 * I computed the effective temperature for a random sample of interpolated
662 * atmospheres by integrating the flux as shown above and compared the results
663 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
664 *
665 * I found that the average discrepancy was:
666 *
667 * DELTA = -1.16% +/- 0.69% (SOLAR models, sample size 4590)
668 * DELTA = -1.17% +/- 0.70% (HALO models, sample size 4828)
669 *
670 * The most extreme discrepancies for the SOLAR models were
671 * -3.18% <= DELTA <= -0.16%
672 *
673 * The most negative discrepancies were for Teff = 35 kK, log(g) = 3.5
674 * The least negative discrepancies were for Teff = 50 kK, log(g) = 4.1
675 *
676 * The most extreme discrepancies for the HALO models were
677 * -2.90% <= DELTA <= -0.13%
678 *
679 * The most negative discrepancies were for Teff = 35 kK, log(g) = 3.5
680 * The least negative discrepancies were for Teff = 50 kK, log(g) = 4.1
681 *
682 * Since Cloudy checks the scaling elsewhere there is no need to re-scale
683 * things here, but this inaccuracy should be kept in mind since it could
684 * indicate problems with the flux distribution */
685
686 InterpolateGridCoStar( &grid, val, val0_lo, val0_hi );
687
688 FreeGrid( &grid );
689 return rfield.nupper;
690}
691
692/* GridCompile rebin user supplied stellar models to match energy grid of code */
693bool GridCompile(const char *InName)
694{
695 bool lgFail = false;
696 realnum Edges[1];
697 string OutName( InName );
698
699 DEBUG_ENTRY( "GridCompile()" );
700
701 fprintf( ioQQQ, " GridCompile on the job.\n" );
702
703 // replace filename extension with ".mod"
704 string::size_type ptr = OutName.find( '.' );
705 ASSERT( ptr != string::npos );
706 OutName.replace( ptr, string::npos, ".mod" );
707
708 process_counter dum;
709 lgFail = lgCompileAtmosphere( InName, OutName.c_str(), Edges, 0L, dum );
710
711 if( !lgFail )
712 {
714
715 /* the file must be local */
716 grid.name = OutName;
717 grid.scheme = AS_LOCAL_ONLY;
718 grid.ident = "bogus ident.";
719 grid.command = "bogus command.";
720
721 InitGrid( &grid, false );
722
723 /* check whether the models in the grid have the correct effective temperature */
724
725 if( strcmp( grid.names[0], "Teff" ) == 0 )
726 {
727 fprintf( ioQQQ, " GridCompile: checking effective temperatures...\n" );
728 ValidateGrid( &grid, 0.02 );
729 }
730
731 FreeGrid( &grid );
732 }
733
734 return lgFail;
735}
736
737/* GridInterpolate read in and interpolate on user supplied grid of atmospheres */
738long GridInterpolate(double val[], /* val[nval] */
739 long *nval,
740 long *ndim,
741 const char *FileName,
742 bool lgList,
743 double *Tlow,
744 double *Thigh)
745{
746 char chIdent[13];
748
749 DEBUG_ENTRY( "GridInterpolate()" );
750
751 // make filename without extension
752 string chTruncName( FileName );
753 string::size_type ptr = chTruncName.find( '.' );
754 if( ptr != string::npos )
755 chTruncName.replace( ptr, string::npos, "" );
756
757 grid.name = FileName;
758 grid.scheme = AS_DATA_OPTIONAL;
759 /* identification of this atmosphere set, used in
760 * the Cloudy output, *must* be 12 characters long */
761 sprintf( chIdent, "%12.12s", chTruncName.c_str() );
762 grid.ident = chIdent;
763 /* the Cloudy command needed to recompile the binary model file */
764 string chString( "COMPILE STARS \"" + chTruncName + ".ascii\"" );
765 grid.command = chString.c_str();
766
767 InitGrid( &grid, lgList );
768
769 CheckVal( &grid, val, nval, ndim );
770
771 InterpolateRectGrid( &grid, val, Tlow, Thigh );
772
773 FreeGrid( &grid );
774 return rfield.nupper;
775}
776
777/* Kurucz79Compile rebin Kurucz 1979 stellar models to match energy grid of code */
779{
780 realnum Edges[1];
781
782 bool lgFail = false;
783
784 DEBUG_ENTRY( "Kurucz79Compile()" );
785
786 fprintf( ioQQQ, " Kurucz79Compile on the job.\n" );
787
788 /* following atmospheres LTE from Kurucz 1979, Ap.J. Sup 40, 1. and
789 * Kurucz (1989) private communication, newer opacities */
790
792
793 if( lgFileReadable( "kurucz79.ascii", pc, as ) && !lgValidBinFile( "kurucz79.mod", pc, as ) )
794 lgFail = lgCompileAtmosphere( "kurucz79.ascii", "kurucz79.mod", Edges, 0L, pc );
795 return lgFail;
796}
797
798/* Kurucz79Interpolate read in and interpolate on Kurucz79 grid of atmospheres */
799long Kurucz79Interpolate(double val[], /* val[nval] */
800 long *nval,
801 long *ndim,
802 bool lgList,
803 double *Tlow,
804 double *Thigh)
805{
807
808 DEBUG_ENTRY( "Kurucz79Interpolate()" );
809
810 grid.name = "kurucz79.mod";
811 grid.scheme = AS_DATA_OPTIONAL;
812 /* identification of this atmosphere set, used in
813 * the Cloudy output, *must* be 12 characters long */
814 grid.ident = " Kurucz 1979";
815 /* the Cloudy command needed to recompile the binary model file */
816 grid.command = "COMPILE STARS";
817
818 InitGrid( &grid, lgList );
819
820 CheckVal( &grid, val, nval, ndim );
821
822 InterpolateRectGrid( &grid, val, Tlow, Thigh );
823
824 FreeGrid( &grid );
825 return rfield.nupper;
826}
827
828/* MihalasCompile rebin Mihalas stellar models to match energy grid of code */
830{
831 realnum Edges[1];
832
833 bool lgFail = false;
834
835 DEBUG_ENTRY( "MihalasCompile()" );
836
837 fprintf( ioQQQ, " MihalasCompile on the job.\n" );
838
839 /* following atmospheres NLTE from Mihalas, NCAR-TN/STR-76 */
840
842
843 if( lgFileReadable( "mihalas.ascii", pc, as ) && !lgValidBinFile( "mihalas.mod", pc, as ) )
844 lgFail = lgCompileAtmosphere( "mihalas.ascii", "mihalas.mod", Edges, 0L, pc );
845 return lgFail;
846}
847
848/* MihalasInterpolate read in and interpolate on Mihalas grid of atmospheres */
849long MihalasInterpolate(double val[], /* val[nval] */
850 long *nval,
851 long *ndim,
852 bool lgList,
853 double *Tlow,
854 double *Thigh)
855{
857
858 DEBUG_ENTRY( "MihalasInterpolate()" );
859
860 grid.name = "mihalas.mod";
861 grid.scheme = AS_DATA_OPTIONAL;
862 /* identification of this atmosphere set, used in
863 * the Cloudy output, *must* be 12 characters long */
864 grid.ident = " Mihalas";
865 /* the Cloudy command needed to recompile the binary model file */
866 grid.command = "COMPILE STARS";
867
868 InitGrid( &grid, lgList );
869
870 CheckVal( &grid, val, nval, ndim );
871
872 InterpolateRectGrid( &grid, val, Tlow, Thigh );
873
874 FreeGrid( &grid );
875 return rfield.nupper;
876}
877
878/* RauchCompile create ascii and mod files for Rauch atmospheres */
880{
881 bool lgFail = false;
882
883 /* these contain frequencies for the major absorption edges */
884 realnum Edges[3];
885
886 /* Before running this program issue the following command where the Rauch
887 * model atmosphere files are kept (0050000_50_solar_bin_0.1 and so on)
888 *
889 * ls *solar_bin_0.1 > rauchmods.list
890 *
891 * and check to see that there are 66 lines in the file.
892 */
893
894 vector<mpp> telg1(NMODS_HCA);
895 vector<mpp> telg2(NMODS_HNI);
896 vector<mpp> telg3(NMODS_PG1159);
897 vector<mpp> telg4(NMODS_HYDR);
898 vector<mpp> telg5(NMODS_HELIUM);
899 vector<mpp> telg6(NMODS_HpHE);
900
901 /* metalicities of the solar and halo grid */
902 static const double par2[2] = { 0., -1. };
903
904 /* Helium fraction by mass */
905 static const double par3[11] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
906
907 DEBUG_ENTRY( "RauchCompile()" );
908
909 fprintf( ioQQQ, " RauchCompile on the job.\n" );
910
911 RauchReadMPP( telg1, telg2, telg3, telg4, telg5, telg6 );
912
913 process_counter dum;
915
916 /* this is the H-Ca grid */
917 if( lgFileReadable( "0050000_50_solar_bin_0.1", dum, as ) && !lgValidAsciiFile( "rauch_h-ca_solar.ascii", as ) )
918 {
919 fprintf( ioQQQ, " Creating rauch_h-ca_solar.ascii....\n" );
920 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_solar.ascii", "_solar_bin_0.1",
921 telg1, NMODS_HCA, 1, 1, par2, 1 );
922 }
923
924 if( lgFileReadable( "0050000_50_halo__bin_0.1", dum, as ) && !lgValidAsciiFile( "rauch_h-ca_halo.ascii", as ) )
925 {
926 fprintf( ioQQQ, " Creating rauch_h-ca_halo.ascii....\n" );
927 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_halo.ascii", "_halo__bin_0.1",
928 telg1, NMODS_HCA, 1, 1, par2, 1 );
929 }
930
931 if( lgFileReadable( "0050000_50_solar_bin_0.1", dum, as ) &&
932 lgFileReadable( "0050000_50_halo__bin_0.1", dum, as ) &&
933 !lgValidAsciiFile( "rauch_h-ca_3d.ascii", as ) )
934 {
935 fprintf( ioQQQ, " Creating rauch_h-ca_3d.ascii....\n" );
936 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_3d.ascii", "_solar_bin_0.1",
937 telg1, NMODS_HCA, 1, 2, par2, 1 );
938 lgFail = lgFail || RauchInitializeSub( "rauch_h-ca_3d.ascii", "_halo__bin_0.1",
939 telg1, NMODS_HCA, 2, 2, par2, 1 );
940 }
941
942 /* this is the H-Ni grid */
943 if( lgFileReadable( "0050000_50_solar_iron.bin_0.1", dum, as ) &&
944 !lgValidAsciiFile( "rauch_h-ni_solar.ascii", as ) )
945 {
946 fprintf( ioQQQ, " Creating rauch_h-ni_solar.ascii....\n" );
947 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_solar.ascii", "_solar_iron.bin_0.1",
948 telg2, NMODS_HNI, 1, 1, par2, 1 );
949 }
950
951 if( lgFileReadable( "0050000_50_halo__iron.bin_0.1", dum, as ) &&
952 !lgValidAsciiFile( "rauch_h-ni_halo.ascii", as ) )
953 {
954 fprintf( ioQQQ, " Creating rauch_h-ni_halo.ascii....\n" );
955 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_halo.ascii", "_halo__iron.bin_0.1",
956 telg2, NMODS_HNI, 1, 1, par2, 1 );
957 }
958
959 if( lgFileReadable( "0050000_50_solar_iron.bin_0.1", dum, as ) &&
960 lgFileReadable( "0050000_50_halo__iron.bin_0.1", dum, as ) &&
961 !lgValidAsciiFile( "rauch_h-ni_3d.ascii", as ) )
962 {
963 fprintf( ioQQQ, " Creating rauch_h-ni_3d.ascii....\n" );
964 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_3d.ascii", "_solar_iron.bin_0.1",
965 telg2, NMODS_HNI, 1, 2, par2, 1 );
966 lgFail = lgFail || RauchInitializeSub( "rauch_h-ni_3d.ascii", "_halo__iron.bin_0.1",
967 telg2, NMODS_HNI, 2, 2, par2, 1 );
968 }
969
970 /* this is the hydrogen deficient PG1159 grid */
971 if( lgFileReadable( "0040000_5.00_33_50_02_15.bin_0.1", dum, as ) &&
972 !lgValidAsciiFile( "rauch_pg1159.ascii", as ) )
973 {
974 fprintf( ioQQQ, " Creating rauch_pg1159.ascii....\n" );
975 lgFail = lgFail || RauchInitializeSub( "rauch_pg1159.ascii", "_33_50_02_15.bin_0.1",
976 telg3, NMODS_PG1159, 1, 1, par2, 2 );
977 }
978
979 /* this is the pure hydrogen grid */
980 if( lgFileReadable( "0020000_4.00_H_00005-02000A.bin_0.1", dum, as ) &&
981 !lgValidAsciiFile( "rauch_hydr.ascii", as ) )
982 {
983 fprintf( ioQQQ, " Creating rauch_hydr.ascii....\n" );
984 lgFail = lgFail || RauchInitializeSub( "rauch_hydr.ascii", "_H_00005-02000A.bin_0.1",
985 telg4, NMODS_HYDR, 1, 1, par2, 2 );
986 }
987
988 /* this is the pure helium grid */
989 if( lgFileReadable( "0050000_5.00_He_00005-02000A.bin_0.1", dum, as ) &&
990 !lgValidAsciiFile( "rauch_helium.ascii", as ) )
991 {
992 fprintf( ioQQQ, " Creating rauch_helium.ascii....\n" );
993 lgFail = lgFail || RauchInitializeSub( "rauch_helium.ascii", "_He_00005-02000A.bin_0.1",
994 telg5, NMODS_HELIUM, 1, 1, par2, 2 );
995 }
996
997 /* this is the 3D grid for arbitrary H+He mixtures */
998 if( lgFileReadable( "0050000_5.00_H+He_1.000_0.000_00005-02000A.bin_0.1", dum, as ) &&
999 !lgValidAsciiFile( "rauch_h+he_3d.ascii", as ) )
1000 {
1001 fprintf( ioQQQ, " Creating rauch_h+he_3d.ascii....\n" );
1002 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_1.000_0.000_00005-02000A.bin_0.1",
1003 telg6, NMODS_HpHE, 1, 11, par3, 2 );
1004 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.900_0.100_00005-02000A.bin_0.1",
1005 telg6, NMODS_HpHE, 2, 11, par3, 2 );
1006 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.800_0.200_00005-02000A.bin_0.1",
1007 telg6, NMODS_HpHE, 3, 11, par3, 2 );
1008 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.700_0.300_00005-02000A.bin_0.1",
1009 telg6, NMODS_HpHE, 4, 11, par3, 2 );
1010 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.600_0.400_00005-02000A.bin_0.1",
1011 telg6, NMODS_HpHE, 5, 11, par3, 2 );
1012 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.500_0.500_00005-02000A.bin_0.1",
1013 telg6, NMODS_HpHE, 6, 11, par3, 2 );
1014 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.400_0.600_00005-02000A.bin_0.1",
1015 telg6, NMODS_HpHE, 7, 11, par3, 2 );
1016 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.300_0.700_00005-02000A.bin_0.1",
1017 telg6, NMODS_HpHE, 8, 11, par3, 2 );
1018 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.200_0.800_00005-02000A.bin_0.1",
1019 telg6, NMODS_HpHE, 9, 11, par3, 2 );
1020 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.100_0.900_00005-02000A.bin_0.1",
1021 telg6, NMODS_HpHE, 10, 11, par3, 2 );
1022 lgFail = lgFail || RauchInitializeSub( "rauch_h+he_3d.ascii", "_H+He_0.000_1.000_00005-02000A.bin_0.1",
1023 telg6, NMODS_HpHE, 11, 11, par3, 2 );
1024 }
1025
1026 /* define the major absorption edges that require special attention during rebinning
1027 *
1028 * NB the frequencies should be chosen here such that they are somewhere in between
1029 * the two frequency points that straddle the edge in the atmosphere model, the
1030 * software in RebinAtmosphere will seek out the exact values of those two points
1031 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
1032 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
1033 *
1034 * NB beware not to choose edges too close to one another (i.e. on the order of the
1035 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
1036 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
1037 * almost certainly lead to erroneous behaviour in RebinAtmosphere */
1038 Edges[0] = 0.99946789f;
1039 Edges[1] = 1.8071406f;
1040 Edges[2] = 3.9996377f;
1041
1042 if( lgFileReadable( "rauch_h-ca_solar.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_solar.mod", pc, as ) )
1043 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_solar.ascii", "rauch_h-ca_solar.mod",Edges,3L, pc );
1044 if( lgFileReadable( "rauch_h-ca_halo.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_halo.mod", pc, as ) )
1045 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_halo.ascii", "rauch_h-ca_halo.mod", Edges, 3L, pc );
1046 if( lgFileReadable( "rauch_h-ca_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ca_3d.mod", pc, as ) )
1047 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ca_3d.ascii", "rauch_h-ca_3d.mod", Edges, 3L, pc );
1048
1049 if( lgFileReadable( "rauch_h-ni_solar.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_solar.mod", pc, as ) )
1050 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_solar.ascii", "rauch_h-ni_solar.mod",Edges,3L, pc );
1051 if( lgFileReadable( "rauch_h-ni_halo.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_halo.mod", pc, as ) )
1052 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_halo.ascii", "rauch_h-ni_halo.mod", Edges, 3L, pc );
1053 if( lgFileReadable( "rauch_h-ni_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h-ni_3d.mod", pc, as ) )
1054 lgFail = lgFail || lgCompileAtmosphere( "rauch_h-ni_3d.ascii", "rauch_h-ni_3d.mod", Edges, 3L, pc );
1055
1056 if( lgFileReadable( "rauch_pg1159.ascii", pc, as ) && !lgValidBinFile( "rauch_pg1159.mod", pc, as ) )
1057 lgFail = lgFail || lgCompileAtmosphere( "rauch_pg1159.ascii", "rauch_pg1159.mod", Edges, 3L, pc );
1058 if( lgFileReadable( "rauch_cowd.ascii", pc, as ) && !lgValidBinFile( "rauch_cowd.mod", pc, as ) )
1059 lgFail = lgFail || lgCompileAtmosphere( "rauch_cowd.ascii", "rauch_cowd.mod", Edges, 3L, pc );
1060
1061 if( lgFileReadable( "rauch_hydr.ascii", pc, as ) && !lgValidBinFile( "rauch_hydr.mod", pc, as ) )
1062 lgFail = lgFail || lgCompileAtmosphere( "rauch_hydr.ascii", "rauch_hydr.mod", Edges, 3L, pc );
1063
1064 if( lgFileReadable( "rauch_helium.ascii", pc, as ) && !lgValidBinFile( "rauch_helium.mod", pc, as ) )
1065 lgFail = lgFail || lgCompileAtmosphere( "rauch_helium.ascii", "rauch_helium.mod", Edges, 3L, pc );
1066
1067 if( lgFileReadable( "rauch_h+he_3d.ascii", pc, as ) && !lgValidBinFile( "rauch_h+he_3d.mod", pc, as ) )
1068 lgFail = lgFail || lgCompileAtmosphere( "rauch_h+he_3d.ascii", "rauch_h+he_3d.mod", Edges, 3L, pc );
1069 return lgFail;
1070}
1071
1072/* RauchInterpolateHCa get one of the Rauch H-Ca model atmospheres, originally by K. Volk */
1073long RauchInterpolateHCa(double val[], /* val[nval] */
1074 long *nval,
1075 long *ndim,
1076 bool lgHalo,
1077 bool lgList,
1078 double *Tlow,
1079 double *Thigh)
1080{
1082
1083 DEBUG_ENTRY( "RauchInterpolateHCa()" );
1084
1085 if( *ndim == 3 )
1086 grid.name = "rauch_h-ca_3d.mod";
1087 else
1088 grid.name = ( lgHalo ? "rauch_h-ca_halo.mod" : "rauch_h-ca_solar.mod" );
1089 grid.scheme = AS_DATA_OPTIONAL;
1090 /* identification of this atmosphere set, used in
1091 * the Cloudy output, *must* be 12 characters long */
1092 grid.ident = " H-Ca Rauch";
1093 /* the Cloudy command needed to recompile the binary model file */
1094 grid.command = "COMPILE STARS";
1095
1096 InitGrid( &grid, lgList );
1097
1098 CheckVal( &grid, val, nval, ndim );
1099
1100 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1101
1102 FreeGrid( &grid );
1103 return rfield.nupper;
1104}
1105
1106/* RauchInterpolateHNi get one of the Rauch H-Ni model atmospheres */
1107long RauchInterpolateHNi(double val[], /* val[nval] */
1108 long *nval,
1109 long *ndim,
1110 bool lgHalo,
1111 bool lgList,
1112 double *Tlow,
1113 double *Thigh)
1114{
1116
1117 DEBUG_ENTRY( "RauchInterpolateHNi()" );
1118
1119 if( *ndim == 3 )
1120 grid.name = "rauch_h-ni_3d.mod";
1121 else
1122 grid.name = ( lgHalo ? "rauch_h-ni_halo.mod" : "rauch_h-ni_solar.mod" );
1123 grid.scheme = AS_DATA_OPTIONAL;
1124 /* identification of this atmosphere set, used in
1125 * the Cloudy output, *must* be 12 characters long */
1126 grid.ident = " H-Ni Rauch";
1127 /* the Cloudy command needed to recompile the binary model file */
1128 grid.command = "COMPILE STARS";
1129
1130 InitGrid( &grid, lgList );
1131
1132 CheckVal( &grid, val, nval, ndim );
1133
1134 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1135
1136 FreeGrid( &grid );
1137 return rfield.nupper;
1138}
1139
1140/* RauchInterpolatePG1159 get one of the Rauch PG1159 model atmospheres */
1141long RauchInterpolatePG1159(double val[], /* val[nval] */
1142 long *nval,
1143 long *ndim,
1144 bool lgList,
1145 double *Tlow,
1146 double *Thigh)
1147{
1149
1150 DEBUG_ENTRY( "RauchInterpolatePG1159()" );
1151
1152 grid.name = "rauch_pg1159.mod";
1153 grid.scheme = AS_DATA_OPTIONAL;
1154 /* identification of this atmosphere set, used in
1155 * the Cloudy output, *must* be 12 characters long */
1156 grid.ident = "PG1159 Rauch";
1157 /* the Cloudy command needed to recompile the binary model file */
1158 grid.command = "COMPILE STARS";
1159
1160 InitGrid( &grid, lgList );
1161
1162 CheckVal( &grid, val, nval, ndim );
1163
1164 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1165
1166 FreeGrid( &grid );
1167 return rfield.nupper;
1168}
1169
1170/* RauchInterpolateCOWD get one of the Rauch C/O white dwarf model atmospheres */
1171long RauchInterpolateCOWD(double val[], /* val[nval] */
1172 long *nval,
1173 long *ndim,
1174 bool lgList,
1175 double *Tlow,
1176 double *Thigh)
1177{
1179
1180 DEBUG_ENTRY( "RauchInterpolateCOWD()" );
1181
1182 grid.name = "rauch_cowd.mod";
1183 grid.scheme = AS_DATA_OPTIONAL;
1184 /* identification of this atmosphere set, used in
1185 * the Cloudy output, *must* be 12 characters long */
1186 grid.ident = "C/O WD Rauch";
1187 /* the Cloudy command needed to recompile the binary model file */
1188 grid.command = "COMPILE STARS";
1189
1190 InitGrid( &grid, lgList );
1191
1192 CheckVal( &grid, val, nval, ndim );
1193
1194 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1195
1196 FreeGrid( &grid );
1197 return rfield.nupper;
1198}
1199
1200/* RauchInterpolateHydr get one of the Rauch pure hydrogen model atmospheres */
1201long RauchInterpolateHydr(double val[], /* val[nval] */
1202 long *nval,
1203 long *ndim,
1204 bool lgList,
1205 double *Tlow,
1206 double *Thigh)
1207{
1209
1210 DEBUG_ENTRY( "RauchInterpolateHydr()" );
1211
1212 grid.name = "rauch_hydr.mod";
1213 grid.scheme = AS_DATA_OPTIONAL;
1214 /* identification of this atmosphere set, used in
1215 * the Cloudy output, *must* be 12 characters long */
1216 grid.ident = " Hydr Rauch";
1217 /* the Cloudy command needed to recompile the binary model file */
1218 grid.command = "COMPILE STARS";
1219
1220 InitGrid( &grid, lgList );
1221
1222 CheckVal( &grid, val, nval, ndim );
1223
1224 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1225
1226 FreeGrid( &grid );
1227 return rfield.nupper;
1228}
1229
1230/* RauchInterpolateHelium get one of the Rauch pure helium model atmospheres */
1231long RauchInterpolateHelium(double val[], /* val[nval] */
1232 long *nval,
1233 long *ndim,
1234 bool lgList,
1235 double *Tlow,
1236 double *Thigh)
1237{
1239
1240 DEBUG_ENTRY( "RauchInterpolateHelium()" );
1241
1242 grid.name = "rauch_helium.mod";
1243 grid.scheme = AS_DATA_OPTIONAL;
1244 /* identification of this atmosphere set, used in
1245 * the Cloudy output, *must* be 12 characters long */
1246 grid.ident = "Helium Rauch";
1247 /* the Cloudy command needed to recompile the binary model file */
1248 grid.command = "COMPILE STARS";
1249
1250 InitGrid( &grid, lgList );
1251
1252 CheckVal( &grid, val, nval, ndim );
1253
1254 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1255
1256 FreeGrid( &grid );
1257 return rfield.nupper;
1258}
1259
1260/* RauchInterpolateHpHe get one of the Rauch hydrogen plus helium model atmospheres */
1261long RauchInterpolateHpHe(double val[], /* val[nval] */
1262 long *nval,
1263 long *ndim,
1264 bool lgList,
1265 double *Tlow,
1266 double *Thigh)
1267{
1269
1270 DEBUG_ENTRY( "RauchInterpolateHpHe()" );
1271
1272 grid.name = "rauch_h+he_3d.mod";
1273 grid.scheme = AS_DATA_OPTIONAL;
1274 /* identification of this atmosphere set, used in
1275 * the Cloudy output, *must* be 12 characters long */
1276 grid.ident = " H+He Rauch";
1277 /* the Cloudy command needed to recompile the binary model file */
1278 grid.command = "COMPILE STARS";
1279
1280 InitGrid( &grid, lgList );
1281
1282 CheckVal( &grid, val, nval, ndim );
1283
1284 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1285
1286 FreeGrid( &grid );
1287 return rfield.nupper;
1288}
1289
1290/* StarburstInitialize does the actual work of preparing the ascii file */
1291bool StarburstInitialize(const char chInName[],
1292 const char chOutName[],
1293 sb_mode mode)
1294{
1295 char chLine[INPUT_LINE_LENGTH]; /* used for getting input lines */
1296
1297 bool lgHeader = true;
1298 long int i, j, nmods, ngp;
1299
1300 size_t nsb_sz = (size_t)NSB99;
1301
1302 double *wavl, *fluxes[MNTS], Age[MNTS], lwavl;
1303
1304 FILE *ioOut, /* pointer to output file we came here to create*/
1305 *ioIn; /* pointer to input files we will read */
1306
1307 DEBUG_ENTRY( "StarburstInitialize()" );
1308
1309 for( i=0; i < MNTS; i++ )
1310 fluxes[i] = NULL;
1311
1312 /* grab some space for the wavelengths and fluxes */
1313 wavl = (double *)MALLOC( sizeof(double)*nsb_sz);
1314
1315 ioIn = open_data( chInName, "r", AS_LOCAL_ONLY );
1316
1317 lwavl = 0.;
1318 nmods = 0;
1319 ngp = 0;
1320
1321 while( read_whole_line( chLine, INPUT_LINE_LENGTH, ioIn ) != NULL )
1322 {
1323 if( !lgHeader )
1324 {
1325 double cage, cwavl, cfl, cfl1, cfl2, cfl3;
1326
1327 /* format: age/yr wavl/Angstrom log10(flux_total) log10(flux_stellar) log10(flux_neb) */
1328 /* we are only interested in the total flux, so we ignore the remaining numbers */
1329 if( sscanf( chLine, " %le %le %le %le %le", &cage, &cwavl, &cfl1, &cfl2, &cfl3 ) != 5 )
1330 {
1331 fprintf( ioQQQ, "syntax error in data of Starburst grid.\n" );
1332 goto error;
1333 }
1334
1335 if( mode == SB_TOTAL )
1336 cfl = cfl1;
1337 else if( mode == SB_STELLAR )
1338 cfl = cfl2;
1339 else if( mode == SB_NEBULAR )
1340 cfl = cfl3;
1341 else
1342 TotalInsanity();
1343
1344 if( cwavl < lwavl )
1345 {
1346 ++nmods;
1347 ngp = 0;
1348
1349 if( nmods >= MNTS )
1350 {
1351 fprintf( ioQQQ, "too many time steps in Starburst grid.\n" );
1352 fprintf( ioQQQ, "please increase MNTS and recompile.\n" );
1353 goto error;
1354 }
1355 }
1356
1357 if( ngp == 0 )
1358 {
1359 fluxes[nmods] = (double *)MALLOC( sizeof(double)*nsb_sz);
1360 Age[nmods] = cage;
1361 }
1362
1363 if( ngp >= (long)nsb_sz )
1364 {
1365 /* this should only be needed when nmods == 0 */
1366 ASSERT( nmods == 0 );
1367
1368 nsb_sz *= 2;
1369 fluxes[0] = (double *)REALLOC(fluxes[0],sizeof(double)*nsb_sz);
1370 wavl = (double *)REALLOC(wavl,sizeof(double)*nsb_sz);
1371 }
1372
1373 if( !fp_equal(Age[nmods],cage,10) )
1374 {
1375 fprintf( ioQQQ, "age error in Starburst grid.\n" );
1376 goto error;
1377 }
1378
1379 if( nmods == 0 )
1380 wavl[ngp] = cwavl;
1381 else
1382 {
1383 if( !fp_equal(wavl[ngp],cwavl,10) )
1384 {
1385 fprintf( ioQQQ, "wavelength error in Starburst grid.\n" );
1386 goto error;
1387 }
1388 }
1389
1390 /* arbitrarily renormalize to flux in erg/cm^2/s/A at 1kpc */
1391 /* constant is log10( 4*pi*(kpc/cm)^2 ) */
1392 fluxes[nmods][ngp] = pow( 10., cfl - 44.077911 );
1393
1394 lwavl = cwavl;
1395 ++ngp;
1396 }
1397
1398 if( lgHeader && strncmp( &chLine[1], "TIME [YR]", 9 ) == 0 )
1399 lgHeader = false;
1400 }
1401
1402 if( lgHeader )
1403 {
1404 /* this happens when the "TIME [YR]" string was not found in column 1 of the file */
1405 fprintf( ioQQQ, "syntax error in header of Starburst grid.\n" );
1406 goto error;
1407 }
1408
1409 ++nmods;
1410
1411 /* finished - close the unit */
1412 fclose(ioIn);
1413
1414 ioOut = open_data( chOutName, "w", AS_LOCAL_ONLY );
1415
1416 fprintf( ioOut, " %ld\n", VERSION_ASCII );
1417 fprintf( ioOut, " %d\n", 1 );
1418 fprintf( ioOut, " %d\n", 1 );
1419 fprintf( ioOut, " Age\n" );
1420 fprintf( ioOut, " %ld\n", nmods );
1421 fprintf( ioOut, " %ld\n", ngp );
1422 /* Starburst99 models give the wavelength in Angstrom */
1423 fprintf( ioOut, " lambda\n" );
1424 /* conversion factor to Angstrom */
1425 fprintf( ioOut, " %.8e\n", 1. );
1426 /* Starburst99 models give the total flux F_lambda in erg/s/A, will be renormalized at 1 kpc */
1427 fprintf( ioOut, " F_lambda\n" );
1428 /* this factor is irrelevant since Teff check will not be carried out */
1429 fprintf( ioOut, " %.8e\n", 1. );
1430 /* write out the Ages */
1431 for( i=0; i < nmods; i++ )
1432 {
1433 fprintf( ioOut, " %.3e", Age[i] );
1434 if( ((i+1)%4) == 0 )
1435 fprintf( ioOut, "\n" );
1436 }
1437 if( (i%4) != 0 )
1438 fprintf( ioOut, "\n" );
1439
1440 fprintf( ioQQQ, " Writing: " );
1441
1442 /* write out the wavelength grid */
1443 for( j=0; j < ngp; j++ )
1444 {
1445 fprintf( ioOut, " %.4e", wavl[j] );
1446 /* want to have 5 numbers per line */
1447 if( ((j+1)%5) == 0 )
1448 fprintf( ioOut, "\n" );
1449 }
1450 /* need to throw a newline if we did not end on an exact line */
1451 if( (j%5) != 0 )
1452 fprintf( ioOut, "\n" );
1453
1454 /* print to screen to show progress */
1455 fprintf( ioQQQ, "." );
1456 fflush( ioQQQ );
1457
1458 for( i=0; i < nmods; i++ )
1459 {
1460 for( j=0; j < ngp; j++ )
1461 {
1462 fprintf( ioOut, " %.4e", fluxes[i][j] );
1463 /* want to have 5 numbers per line */
1464 if( ((j+1)%5) == 0 )
1465 fprintf( ioOut, "\n" );
1466 }
1467 /* need to throw a newline if we did not end on an exact line */
1468 if( (j%5) != 0 )
1469 fprintf( ioOut, "\n" );
1470
1471 /* print to screen to show progress */
1472 fprintf( ioQQQ, "." );
1473 fflush( ioQQQ );
1474 }
1475
1476 fprintf( ioQQQ, " done.\n" );
1477
1478 fclose(ioOut);
1479
1480 /* free the space grabbed above */
1481 for( i=0; i < MNTS; i++ )
1482 FREE_SAFE( fluxes[i] );
1483 FREE_CHECK( wavl );
1484 return false;
1485
1486error:
1487 for( i=0; i < MNTS; i++ )
1488 FREE_SAFE( fluxes[i] );
1489 FREE_CHECK( wavl );
1490 return true;
1491}
1492
1493/* StarburstCompile, rebin Starburst99 model output to match energy grid of code */
1495{
1496 realnum Edges[1];
1497
1498 bool lgFail = false;
1499
1500 DEBUG_ENTRY( "StarburstCompile()" );
1501
1502 fprintf( ioQQQ, " StarburstCompile on the job.\n" );
1503
1504 process_counter dum;
1506
1507 if( lgFileReadable( "starburst99.stb99", dum, as ) && !lgValidAsciiFile( "starburst99.ascii", as ) )
1508 lgFail = lgFail || StarburstInitialize( "starburst99.stb99", "starburst99.ascii", SB_TOTAL );
1509 if( lgFileReadable( "starburst99.ascii", pc, as ) && !lgValidBinFile( "starburst99.mod", pc, as ) )
1510 lgFail = lgFail || lgCompileAtmosphere( "starburst99.ascii", "starburst99.mod", Edges, 0L, pc );
1511
1512 if( lgFileReadable( "starburst99_2d.ascii", pc, as ) && !lgValidBinFile( "starburst99_2d.mod", pc, as ) )
1513 lgFail = lgFail || lgCompileAtmosphere( "starburst99_2d.ascii", "starburst99_2d.mod", Edges, 0L, pc );
1514 return lgFail;
1515}
1516
1517/* TlustyCompile rebin Tlusty BSTAR2006/OSTAR2002 stellar models to match energy grid of code */
1519{
1520 /* these contain frequencies for the major absorption edges */
1521 realnum Edges[1];
1522
1523 bool lgFail = false;
1524
1525 DEBUG_ENTRY( "TlustyCompile()" );
1526
1527 fprintf( ioQQQ, " TlustyCompile on the job.\n" );
1528
1530
1531 if( lgFileReadable( "obstar_merged_p03.ascii", pc, as ) && !lgValidBinFile( "obstar_merged_p03.mod", pc, as ) )
1532 lgFail = lgFail || lgCompileAtmosphere( "obstar_merged_p03.ascii", "obstar_merged_p03.mod", Edges, 0L, pc );
1533 if( lgFileReadable( "obstar_merged_p00.ascii", pc, as ) && !lgValidBinFile( "obstar_merged_p00.mod", pc, as ) )
1534 lgFail = lgFail || lgCompileAtmosphere( "obstar_merged_p00.ascii", "obstar_merged_p00.mod", Edges, 0L, pc );
1535 if( lgFileReadable( "obstar_merged_m03.ascii", pc, as ) && !lgValidBinFile( "obstar_merged_m03.mod", pc, as ) )
1536 lgFail = lgFail || lgCompileAtmosphere( "obstar_merged_m03.ascii", "obstar_merged_m03.mod", Edges, 0L, pc );
1537 if( lgFileReadable( "obstar_merged_m07.ascii", pc, as ) && !lgValidBinFile( "obstar_merged_m07.mod", pc, as ) )
1538 lgFail = lgFail || lgCompileAtmosphere( "obstar_merged_m07.ascii", "obstar_merged_m07.mod", Edges, 0L, pc );
1539 if( lgFileReadable( "obstar_merged_m10.ascii", pc, as ) && !lgValidBinFile( "obstar_merged_m10.mod", pc, as ) )
1540 lgFail = lgFail || lgCompileAtmosphere( "obstar_merged_m10.ascii", "obstar_merged_m10.mod", Edges, 0L, pc );
1541 if( lgFileReadable( "obstar_merged_m99.ascii", pc, as ) && !lgValidBinFile( "obstar_merged_m99.mod", pc, as ) )
1542 lgFail = lgFail || lgCompileAtmosphere( "obstar_merged_m99.ascii", "obstar_merged_m99.mod", Edges, 0L, pc );
1543
1544 if( lgFileReadable( "obstar_merged_3d.ascii", pc, as ) && !lgValidBinFile( "obstar_merged_3d.mod", pc, as ) )
1545 lgFail = lgFail || lgCompileAtmosphere( "obstar_merged_3d.ascii", "obstar_merged_3d.mod", Edges, 0L, pc );
1546
1547 if( lgFileReadable( "bstar2006_p03.ascii", pc, as ) && !lgValidBinFile( "bstar2006_p03.mod", pc, as ) )
1548 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_p03.ascii", "bstar2006_p03.mod", Edges, 0L, pc );
1549 if( lgFileReadable( "bstar2006_p00.ascii", pc, as ) && !lgValidBinFile( "bstar2006_p00.mod", pc, as ) )
1550 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_p00.ascii", "bstar2006_p00.mod", Edges, 0L, pc );
1551 if( lgFileReadable( "bstar2006_m03.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m03.mod", pc, as ) )
1552 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m03.ascii", "bstar2006_m03.mod", Edges, 0L, pc );
1553 if( lgFileReadable( "bstar2006_m07.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m07.mod", pc, as ) )
1554 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m07.ascii", "bstar2006_m07.mod", Edges, 0L, pc );
1555 if( lgFileReadable( "bstar2006_m10.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m10.mod", pc, as ) )
1556 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m10.ascii", "bstar2006_m10.mod", Edges, 0L, pc );
1557 if( lgFileReadable( "bstar2006_m99.ascii", pc, as ) && !lgValidBinFile( "bstar2006_m99.mod", pc, as ) )
1558 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_m99.ascii", "bstar2006_m99.mod", Edges, 0L, pc );
1559
1560 if( lgFileReadable( "bstar2006_3d.ascii", pc, as ) && !lgValidBinFile( "bstar2006_3d.mod", pc, as ) )
1561 lgFail = lgFail || lgCompileAtmosphere( "bstar2006_3d.ascii", "bstar2006_3d.mod", Edges, 0L, pc );
1562
1563 if( lgFileReadable( "ostar2002_p03.ascii", pc, as ) && !lgValidBinFile( "ostar2002_p03.mod", pc, as ) )
1564 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_p03.ascii", "ostar2002_p03.mod", Edges, 0L, pc );
1565 if( lgFileReadable( "ostar2002_p00.ascii", pc, as ) && !lgValidBinFile( "ostar2002_p00.mod", pc, as ) )
1566 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_p00.ascii", "ostar2002_p00.mod", Edges, 0L, pc );
1567 if( lgFileReadable( "ostar2002_m03.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m03.mod", pc, as ) )
1568 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m03.ascii", "ostar2002_m03.mod", Edges, 0L, pc );
1569 if( lgFileReadable( "ostar2002_m07.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m07.mod", pc, as ) )
1570 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m07.ascii", "ostar2002_m07.mod", Edges, 0L, pc );
1571 if( lgFileReadable( "ostar2002_m10.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m10.mod", pc, as ) )
1572 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m10.ascii", "ostar2002_m10.mod", Edges, 0L, pc );
1573 if( lgFileReadable( "ostar2002_m15.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m15.mod", pc, as ) )
1574 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m15.ascii", "ostar2002_m15.mod", Edges, 0L, pc );
1575 if( lgFileReadable( "ostar2002_m17.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m17.mod", pc, as ) )
1576 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m17.ascii", "ostar2002_m17.mod", Edges, 0L, pc );
1577 if( lgFileReadable( "ostar2002_m20.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m20.mod", pc, as ) )
1578 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m20.ascii", "ostar2002_m20.mod", Edges, 0L, pc );
1579 if( lgFileReadable( "ostar2002_m30.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m30.mod", pc, as ) )
1580 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m30.ascii", "ostar2002_m30.mod", Edges, 0L, pc );
1581 if( lgFileReadable( "ostar2002_m99.ascii", pc, as ) && !lgValidBinFile( "ostar2002_m99.mod", pc, as ) )
1582 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_m99.ascii", "ostar2002_m99.mod", Edges, 0L, pc );
1583
1584 if( lgFileReadable( "ostar2002_3d.ascii", pc, as ) && !lgValidBinFile( "ostar2002_3d.mod", pc, as ) )
1585 lgFail = lgFail || lgCompileAtmosphere( "ostar2002_3d.ascii", "ostar2002_3d.mod", Edges, 0L, pc );
1586 return lgFail;
1587}
1588
1589/* TlustyInterpolate get one of the Tlusty OBSTAR_MERGED/BSTAR2006/OSTAR2002 model atmospheres */
1590long TlustyInterpolate(double val[], /* val[nval] */
1591 long *nval,
1592 long *ndim,
1593 tl_grid tlg,
1594 const char *chMetalicity,
1595 bool lgList,
1596 double *Tlow,
1597 double *Thigh)
1598{
1599 char chIdent[13];
1601
1602 DEBUG_ENTRY( "TlustyInterpolate()" );
1603
1604 if( tlg == TL_OBSTAR )
1605 grid.name = "obstar_merged_";
1606 else if( tlg == TL_BSTAR )
1607 grid.name = "bstar2006_";
1608 else if( tlg == TL_OSTAR )
1609 grid.name = "ostar2002_";
1610 else
1611 TotalInsanity();
1612 if( *ndim == 3 )
1613 grid.name += "3d";
1614 else
1615 grid.name += chMetalicity;
1616 grid.name += ".mod";
1617 grid.scheme = AS_DATA_OPTIONAL;
1618 /* identification of this atmosphere set, used in
1619 * the Cloudy output, *must* be 12 characters long */
1620 if( *ndim == 3 )
1621 {
1622 strcpy( chIdent, "3-dim" );
1623 }
1624 else
1625 {
1626 strcpy( chIdent, "Z " );
1627 strcat( chIdent, chMetalicity );
1628 }
1629 if( tlg == TL_OBSTAR )
1630 strcat( chIdent, " OBstar" );
1631 else if( tlg == TL_BSTAR )
1632 strcat( chIdent, " Bstr06" );
1633 else if( tlg == TL_OSTAR )
1634 strcat( chIdent, " Ostr02" );
1635 else
1636 TotalInsanity();
1637 grid.ident = chIdent;
1638 /* the Cloudy command needed to recompile the binary model file */
1639 grid.command = "COMPILE STARS";
1640
1641 InitGrid( &grid, lgList );
1642
1643 CheckVal( &grid, val, nval, ndim );
1644
1645 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1646
1647 FreeGrid( &grid );
1648 return rfield.nupper;
1649}
1650
1651/* WernerCompile rebin Werner stellar models to match energy grid of code */
1652/* >>chng 05 nov 16, added return value to indicate success (0) or failure (1) */
1654{
1655 /* these contain frequencies for the major absorption edges */
1656 realnum Edges[3];
1657
1658 bool lgFail = false;
1659
1660 DEBUG_ENTRY( "WernerCompile()" );
1661
1662 fprintf( ioQQQ, " WernerCompile on the job.\n" );
1663
1664 /* define the major absorption edges that require special attention during rebinning
1665 *
1666 * NB the frequencies should be chosen here such that they are somewhere in between
1667 * the two frequency points that straddle the edge in the atmosphere model, the
1668 * software in RebinAtmosphere will seek out the exact values of those two points
1669 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
1670 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere in between (e.g. at 911.76A).
1671 *
1672 * NB beware not to choose edges too close to one another (i.e. on the order of the
1673 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
1674 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
1675 * almost certainly lead to erroneous behaviour in RebinAtmosphere */
1676 Edges[0] = 0.99946789f;
1677 Edges[1] = 1.8071406f;
1678 Edges[2] = 3.9996377f;
1679
1680 /* The "kwerner.ascii" file is a modified ascii dump of the Klaus Werner
1681 * stellar model files which he gave to me in 1992. The first set of values
1682 * is the frequency grid (in Ryd) followed by the atmosphere models in order
1683 * of increasing temperature and log(g). The following comments are already
1684 * incorporated in the modified kwerner.ascii file that is supplied with Cloudy.
1685 *
1686 * >>chng 00 oct 18, The frequency grid was slightly tweaked compared to the
1687 * original values supplied by Klaus Werner to make it monotonically increasing;
1688 * this is due to there being fluxes above and below certain wavelengths where
1689 * the opacity changes (i.e. the Lyman and Balmer limits for example) which are
1690 * assigned the same wavelength in the original Klaus Werner files. PvH
1691 *
1692 * >>chng 00 oct 20, StarEner[172] is out of sequence. As per the Klaus Werner comment,
1693 * it should be omitted. The energy grid is very dense in this region and was most likely
1694 * intended to sample an absorption line which was not included in this particular grid.
1695 * StarFlux[172] is therefore always equal to the flux in neighbouring points (at least
1696 * those with slightly smaller energies). It is therefore safe to ignore this point. PvH
1697 *
1698 * >>chng 00 oct 20, As per the same comment, StarFlux[172] is also deleted. PvH */
1699
1701
1702 if( lgFileReadable( "kwerner.ascii", pc, as ) && !lgValidBinFile( "kwerner.mod", pc, as ) )
1703 lgFail = lgCompileAtmosphere( "kwerner.ascii", "kwerner.mod", Edges, 3L, pc );
1704 return lgFail;
1705}
1706
1707/* WernerInterpolate read in and interpolate on Werner grid of PN atmospheres, originally by K Volk */
1708long WernerInterpolate(double val[], /* val[nval] */
1709 long *nval,
1710 long *ndim,
1711 bool lgList,
1712 double *Tlow,
1713 double *Thigh)
1714{
1716
1717 DEBUG_ENTRY( "WernerInterpolate()" );
1718
1719 /* This subroutine was added (28 dec 1992) to read from the set of
1720 * hot white dwarf model atmospheres from Klaus Werner at Kiel. The
1721 * values are read in (energy in Rydberg units, f_nu in cgs units)
1722 * for any of the 20 models. Each model had 513 points before rebinning.
1723 * The Rayleigh-Jeans tail was extrapolated. */
1724
1725 grid.name = "kwerner.mod";
1726 grid.scheme = AS_DATA_OPTIONAL;
1727 /* identification of this atmosphere set, used in
1728 * the Cloudy output, *must* be 12 characters long */
1729 grid.ident = "Klaus Werner";
1730 /* the Cloudy command needed to recompile the binary model file */
1731 grid.command = "COMPILE STARS";
1732
1733 InitGrid( &grid, lgList );
1734
1735 CheckVal( &grid, val, nval, ndim );
1736
1737 /* Note on the interpolation: 26 October 2000 (Peter van Hoof)
1738 *
1739 * I computed the effective temperature for a random sample of interpolated
1740 * atmospheres by integrating the flux as shown above and compared the results
1741 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
1742 *
1743 * I found that the average discrepancy was:
1744 *
1745 * DELTA = -0.71% +/- 0.71% (sample size 5000)
1746 *
1747 * The most extreme discrepancies were
1748 * -4.37% <= DELTA <= 0.24%
1749 *
1750 * The most negative discrepancies were for Teff = 95 kK, log(g) = 5
1751 * The most positive discrepancies were for Teff = 160 kK, log(g) = 8
1752 *
1753 * Since Cloudy checks the scaling elsewhere there is no need to re-scale
1754 * things here, but this inaccuracy should be kept in mind since it could
1755 * indicate problems with the flux distribution */
1756
1757 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1758
1759 FreeGrid( &grid );
1760 return rfield.nupper;
1761}
1762
1763/* WMBASICCompile rebin WMBASIC stellar models to match energy grid of code */
1765{
1766 /* these contain frequencies for the major absorption edges */
1767 realnum Edges[3];
1768
1769 bool lgFail = false;
1770
1771 DEBUG_ENTRY( "WMBASICCompile()" );
1772
1773 fprintf( ioQQQ, " WMBASICCompile on the job.\n" );
1774
1775 /* define the major absorption edges that require special attention during rebinning */
1776 Edges[0] = 0.99946789f;
1777 Edges[1] = 1.8071406f;
1778 Edges[2] = 3.9996377f;
1779
1781
1782 if( lgFileReadable( "wmbasic.ascii", pc, as ) && !lgValidBinFile( "wmbasic.mod", pc, as ) )
1783 lgFail = lgCompileAtmosphere( "wmbasic.ascii", "wmbasic.mod", Edges, 3L, pc );
1784 return lgFail;
1785}
1786
1787/* WMBASICInterpolate read in and interpolate on WMBASIC grid of hot star atmospheres */
1788long WMBASICInterpolate(double val[], /* val[nval] */
1789 long *nval,
1790 long *ndim,
1791 bool lgList,
1792 double *Tlow,
1793 double *Thigh)
1794{
1796
1797 DEBUG_ENTRY( "WMBASICInterpolate()" );
1798
1799 grid.name = "wmbasic.mod";
1800 grid.scheme = AS_DATA_OPTIONAL;
1801 /* identification of this atmosphere set, used in
1802 * the Cloudy output, *must* be 12 characters long */
1803 grid.ident = " WMBASIC";
1804 /* the Cloudy command needed to recompile the binary model file */
1805 grid.command = "COMPILE STARS";
1806
1807 InitGrid( &grid, lgList );
1808
1809 CheckVal( &grid, val, nval, ndim );
1810
1811 InterpolateRectGrid( &grid, val, Tlow, Thigh );
1812
1813 FreeGrid( &grid );
1814 return rfield.nupper;
1815}
1816
1817/* CompileAtmosphereCoStar rebin costar stellar atmospheres to match cloudy energy grid,
1818 * called by the compile stars command */
1819STATIC bool lgCompileAtmosphereCoStar(const char chFNameIn[],
1820 const char chFNameOut[],
1821 const realnum Edges[], /* Edges[nedges] */
1822 long nedges,
1823 process_counter& pc)
1824{
1825 char chLine[INPUT_LINE_LENGTH];
1826 char names[MDIM][MNAM+1];
1827 int32 val[7];
1828 uint32 uval[2];
1829 double dval[3];
1830 char md5sum[NMD5];
1831 long int i, j, nskip, nModels, nWL;
1832
1833 /* these will be malloced into large work arrays*/
1834 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL;
1835 /* this will hold all the model parameters */
1836 mpp *telg = NULL;
1837
1838 FILE *ioIN; /* used for input */
1839 FILE *ioOUT; /* used for output */
1840 vector<realnum> SaveAnu(rfield.nupper);
1841
1842 DEBUG_ENTRY( "CompileAtmosphereCoStar()" );
1843
1844 /* This is a program to re-bin the costar stellar model spectra to match the
1845 * Cloudy grid. For short wavelengths I will use a power law extrapolation
1846 * of the model values (which should be falling rapidly) if needed. At long
1847 * wavelengths I will assume Rayleigh-Jeans from the last stellar model point
1848 * to extrapolate to 1 cm wavelength. */
1849
1850 /* This version uses power-law interpolation between the points of the stellar model. */
1851
1852 /* read the original data file obtained off the web,
1853 * open as read only */
1854 try
1855 {
1856 ioIN = open_data( chFNameIn, "r", AS_LOCAL_ONLY );
1857 }
1858 catch( cloudy_exit )
1859 {
1860 goto error;
1861 }
1862 fprintf( ioQQQ, " CompileAtmosphereCoStar got %s.\n", chFNameIn );
1863
1864 /* get first line and see how many more to skip */
1865 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1866 {
1867 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading nskip.\n" );
1868 goto error;
1869 }
1870 sscanf( chLine, "%li", &nskip );
1871
1872 /* now skip the header information */
1873 for( i=0; i < nskip; ++i )
1874 {
1875 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1876 {
1877 fprintf( ioQQQ, " CompileAtmosphereCoStar fails skipping header.\n" );
1878 goto error;
1879 }
1880 }
1881
1882 /* now get number of models and number of wavelengths */
1883 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1884 {
1885 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading nModels, nWL.\n" );
1886 goto error;
1887 }
1888 sscanf( chLine, "%li%li", &nModels, &nWL );
1889
1890 if( nModels <= 0 || nWL <= 0 )
1891 {
1892 fprintf( ioQQQ, " CompileAtmosphereCoStar scanned off impossible values for nModels=%li or nWL=%li\n",
1893 nModels, nWL );
1894 goto error;
1895 }
1896
1897 /* allocate space for the stellar parameters */
1898 telg = (mpp *)CALLOC( (size_t)nModels, sizeof(mpp) );
1899
1900 /* get all model parameters for the atmospheres */
1901 for( i=0; i < nModels; ++i )
1902 {
1903 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1904 {
1905 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading model parameters.\n" );
1906 goto error;
1907 }
1908 /* first letter on line is indicator of grid */
1909 telg[i].chGrid = chLine[0];
1910 /* get the model id number */
1911 sscanf( chLine+1, "%i", &telg[i].modid );
1912 /* get the temperature */
1913 sscanf( chLine+23, "%lg", &telg[i].par[0] );
1914 /* get the surface gravity */
1915 sscanf( chLine+31, "%lg", &telg[i].par[1] );
1916 /* get the ZAMS mass */
1917 sscanf( chLine+7, "%lg", &telg[i].par[2] );
1918 /* get the model age */
1919 sscanf( chLine+15, "%lg", &telg[i].par[3] );
1920
1921 /* the code in parse_table.cpp implicitly depends on this! */
1922 ASSERT( telg[i].par[2] > 10. );
1923 ASSERT( telg[i].par[3] > 10. );
1924
1925 /* convert ZAMS masses to logarithms */
1926 telg[i].par[2] = log10(telg[i].par[2]);
1927 }
1928
1929 /* this will be the file we create, that will be read to compute models,
1930 * open to write binary */
1931 try
1932 {
1933 ioOUT = open_data( chFNameOut, "wb", AS_LOCAL_ONLY );
1934 }
1935 catch( cloudy_exit )
1936 {
1937 goto error;
1938 }
1939
1940 val[0] = (int32)VERSION_BIN;
1941 val[1] = (int32)MDIM;
1942 val[2] = (int32)MNAM;
1943 val[3] = (int32)2;
1944 val[4] = (int32)4;
1945 val[5] = (int32)nModels;
1946 val[6] = (int32)rfield.nupper;
1947 uval[0] = sizeof(val) + sizeof(uval) + sizeof(dval) + sizeof(md5sum) +
1948 sizeof(names) + nModels*sizeof(mpp); /* nOffset */
1949 uval[1] = rfield.nupper*sizeof(realnum); /* nBlocksize */
1950 dval[0] = double(rfield.emm);
1951 dval[1] = double(rfield.egamry);
1952 dval[2] = double(continuum.ResolutionScaleFactor);
1953
1954 strncpy( md5sum, continuum.mesh_md5sum.c_str(), sizeof(md5sum) );
1955
1956 strncpy( names[0], "Teff\0\0", MNAM+1 );
1957 strncpy( names[1], "log(g)", MNAM+1 );
1958 strncpy( names[2], "log(M)", MNAM+1 );
1959 strncpy( names[3], "Age\0\0\0", MNAM+1 );
1960
1961 for (long i=0; i<rfield.nupper; ++i)
1962 SaveAnu[i] = (realnum) rfield.AnuOrg[i];
1963 if( fwrite( val, sizeof(val), 1, ioOUT ) != 1 ||
1964 fwrite( uval, sizeof(uval), 1, ioOUT ) != 1 ||
1965 /* write out the lower, upper bound of the energy mesh, and the res scale factor */
1966 fwrite( dval, sizeof(dval), 1, ioOUT ) != 1 ||
1967 /* write out the (modified) md5 checksum of continuum_mesh.ini */
1968 fwrite( md5sum, sizeof(md5sum), 1, ioOUT ) != 1 ||
1969 fwrite( names, sizeof(names), 1, ioOUT ) != 1 ||
1970 /* write out the array of {Teff,log(g)} pairs */
1971 fwrite( telg, sizeof(mpp), (size_t)nModels, ioOUT ) != (size_t)nModels ||
1972 /* write out the cloudy energy grid for later sanity checks */
1973 fwrite( get_ptr(SaveAnu), (size_t)uval[1], 1, ioOUT ) != 1 )
1974 {
1975 fprintf( ioQQQ, " CompileAtmosphereCoStar failed writing header of output file.\n" );
1976 goto error;
1977 }
1978
1979 /* MALLOC some workspace */
1980 StarEner = (realnum *)MALLOC( sizeof(realnum)*nWL );
1981 StarFlux = (realnum *)MALLOC( sizeof(realnum)*nWL );
1982 CloudyFlux = (realnum *)MALLOC( (size_t)uval[1] );
1983
1984 fprintf( ioQQQ, " Compiling: " );
1985
1986 /* get the star data */
1987 for( i=0; i < nModels; ++i )
1988 {
1989 /* get number to skip */
1990 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
1991 {
1992 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading the skip to next spectrum.\n" );
1993 goto error;
1994 }
1995 sscanf( chLine, "%li", &nskip );
1996
1997 for( j=0; j < nskip; ++j )
1998 {
1999 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
2000 {
2001 fprintf( ioQQQ, " CompileAtmosphereCoStar fails doing the skip.\n" );
2002 goto error;
2003 }
2004 }
2005
2006 /* now read in the wavelength and flux for this star, read in
2007 * backwards since we want to be in increasing energy order rather
2008 * than wavelength */
2009 for( j=nWL-1; j >= 0; --j )
2010 {
2011 if( read_whole_line( chLine, (int)sizeof(chLine), ioIN ) == NULL )
2012 {
2013 fprintf( ioQQQ, " CompileAtmosphereCoStar fails reading the spectral data.\n" );
2014 goto error;
2015 }
2016 double help1, help2;
2017 sscanf( chLine, "%lg %lg", &help1, &help2 );
2018
2019 /* continuum flux was log, convert to linear, also do
2020 * conversion from "astrophysical" flux to F_nu in cgs units */
2021 StarFlux[j] = (realnum)(PI*pow(10.,help2));
2022 /* StarEner was in Angstroms, convert to Ryd */
2023 StarEner[j] = (realnum)(RYDLAM/help1);
2024
2025 /* sanity check */
2026 if( j < nWL-1 )
2027 ASSERT( StarEner[j] < StarEner[j+1] );
2028 }
2029
2030 /* this will do the heavy lifting, and define arrays used below,
2031 * NB the lowest energy point in these grids appears to be bogus.
2032 * tell rebin about nWL-1 */
2033 RebinAtmosphere(nWL-1, StarEner+1, StarFlux+1, CloudyFlux, nedges, Edges );
2034
2035 /* write the continuum out as a binary file */
2036 if( fwrite( CloudyFlux, (size_t)uval[1], 1, ioOUT ) != 1 )
2037 {
2038 fprintf( ioQQQ, " CompileAtmosphereCoStar failed writing star flux.\n" );
2039 goto error;
2040 }
2041
2042 fprintf( ioQQQ, "." );
2043 fflush( ioQQQ );
2044 }
2045
2046 fprintf( ioQQQ, " done.\n" );
2047
2048 fclose( ioIN );
2049 fclose( ioOUT );
2050
2051 FREE_CHECK( telg );
2052 FREE_CHECK( StarEner );
2053 FREE_CHECK( StarFlux );
2054 FREE_CHECK( CloudyFlux );
2055
2056 fprintf( ioQQQ, "\n CompileAtmosphereCoStar completed ok.\n\n" );
2057 ++pc.nOK;
2058 return false;
2059
2060error:
2061 FREE_SAFE( telg );
2062 FREE_SAFE( StarEner );
2063 FREE_SAFE( StarFlux );
2064 FREE_SAFE( CloudyFlux );
2065 ++pc.nFail;
2066 return true;
2067}
2068
2069/* InterpolateGridCoStar read in and interpolate on costar grid of windy O atmospheres */
2070STATIC void InterpolateGridCoStar(const stellar_grid *grid, /* struct with all the grid parameters */
2071 const double val[], /* val[0]: Teff for imode = 1,2; M_ZAMS for imode = 3;
2072 * age for imode = 4 */
2073 /* val[1]: nmodid for imode = 1; log(g) for imode = 2;
2074 * age for imode = 3; M_ZAMS for imode = 4 */
2075 double *val0_lo,
2076 double *val0_hi)
2077{
2078 long i, j, k, nmodid, off, ptr;
2079 long *indloTr, *indhiTr, useTr[2];
2080 long indlo[2], indhi[2], index[2];
2081 realnum *ValTr;
2082 double lval[2], aval[4];
2083
2084 DEBUG_ENTRY( "InterpolateGridCoStar()" );
2085
2086 switch( grid->imode )
2087 {
2090 lval[0] = val[0];
2091 lval[1] = val[1];
2092 off = 0;
2093 break;
2095 lval[0] = log10(val[0]); /* use log10(M_ZAMS) internally */
2096 lval[1] = val[1];
2097 off = 2;
2098 break;
2100 /* swap parameters, hence mimic IM_COSTAR_MZAMS_AGE */
2101 lval[0] = log10(val[1]); /* use log10(M_ZAMS) internally */
2102 lval[1] = val[0];
2103 off = 2;
2104 break;
2105 default:
2106 fprintf( ioQQQ, " InterpolateGridCoStar called with insane value for imode: %d.\n", grid->imode );
2108 }
2109
2110 nmodid = (long)(lval[1]+0.5);
2111
2112 ASSERT( rfield.lgContMalloc[rfield.nShape] );
2113
2114 /* read in the saved cloudy energy scale so we can confirm this is a good image */
2115 GetBins( grid, rfield.tNu[rfield.nShape] );
2116
2117# if DEBUGPRT
2118 /* check whether the models in the grid have the correct effective temperature */
2119 ValidateGrid( grid, 0.005 );
2120# endif
2121
2122 /* now allocate some temp workspace */
2123 ValTr = (realnum *)MALLOC( sizeof(realnum)*grid->nTracks );
2124 indloTr = (long *)MALLOC( sizeof(long)*grid->nTracks );
2125 indhiTr = (long *)MALLOC( sizeof(long)*grid->nTracks );
2126
2127 /* first do horizontal search, i.e. search along individual tracks */
2128 for( j=0; j < grid->nTracks; j++ )
2129 {
2130 if( grid->imode == IM_COSTAR_TEFF_MODID )
2131 {
2132 if( grid->trackLen[j] >= nmodid ) {
2133 index[0] = nmodid - 1;
2134 index[1] = j;
2135 ptr = grid->jval[JIndex(grid,index)];
2136 indloTr[j] = ptr;
2137 indhiTr[j] = ptr;
2138 ValTr[j] = (realnum)grid->telg[ptr].par[off];
2139 }
2140 else
2141 {
2142 indloTr[j] = -2;
2143 indhiTr[j] = -2;
2144 ValTr[j] = -FLT_MAX;
2145 }
2146 }
2147 else
2148 {
2149 FindHCoStar( grid, j, lval[1], off, ValTr, indloTr, indhiTr );
2150 }
2151 }
2152
2153# if DEBUGPRT
2154 for( j=0; j < grid->nTracks; j++ )
2155 {
2156 if( indloTr[j] >= 0 )
2157 printf( "track %c: models %c%d, %c%d, val %g\n",
2158 (char)('A'+j), grid->telg[indloTr[j]].chGrid, grid->telg[indloTr[j]].modid,
2159 grid->telg[indhiTr[j]].chGrid, grid->telg[indhiTr[j]].modid, ValTr[j]);
2160 }
2161# endif
2162
2163 /* now do vertical search, i.e. interpolate between tracks */
2164 FindVCoStar( grid, lval[0], ValTr, useTr );
2165
2166 /* This should only happen when InterpolateGridCoStar is called in non-optimizing mode,
2167 * when optimizing InterpolateGridCoStar should report back to optimize_func()...
2168 * The fact that FindVCoStar allows interpolation between non-adjoining tracks
2169 * should guarantee that this will not happen. */
2170 if( useTr[0] < 0 )
2171 {
2172 fprintf( ioQQQ, " The parameters for the requested CoStar model are out of range.\n" );
2174 }
2175
2176 ASSERT( useTr[0] >= 0 && useTr[0] < grid->nTracks );
2177 ASSERT( useTr[1] >= 0 && useTr[1] < grid->nTracks );
2178 ASSERT( indloTr[useTr[0]] >= 0 && indloTr[useTr[0]] < (int)grid->nmods );
2179 ASSERT( indhiTr[useTr[0]] >= 0 && indhiTr[useTr[0]] < (int)grid->nmods );
2180 ASSERT( indloTr[useTr[1]] >= 0 && indloTr[useTr[1]] < (int)grid->nmods );
2181 ASSERT( indhiTr[useTr[1]] >= 0 && indhiTr[useTr[1]] < (int)grid->nmods );
2182
2183# if DEBUGPRT
2184 printf( "interpolate between tracks %c and %c\n", (char)('A'+useTr[0]), (char)('A'+useTr[1]) );
2185# endif
2186
2187 indlo[0] = indloTr[useTr[0]];
2188 indhi[0] = indhiTr[useTr[0]];
2189 indlo[1] = indloTr[useTr[1]];
2190 indhi[1] = indhiTr[useTr[1]];
2191
2192 InterpolateModelCoStar( grid, lval, aval, indlo, indhi, index, 0, off, rfield.tslop[rfield.nShape] );
2193
2194 for( i=0; i < rfield.nupper; i++ )
2195 {
2196 rfield.tslop[rfield.nShape][i] = (realnum)pow((realnum)10.f,rfield.tslop[rfield.nShape][i]);
2197 if( rfield.tslop[rfield.nShape][i] < 1e-37 )
2198 rfield.tslop[rfield.nShape][i] = 0.;
2199 }
2200
2201 if( false )
2202 {
2203 FILE *ioBUG = fopen( "interpolated.txt", "w" );
2204 for( k=0; k < rfield.nupper; ++k )
2205 fprintf( ioBUG, "%e %e\n", rfield.tNu[rfield.nShape][k].Ryd(), rfield.tslop[rfield.nShape][k] );
2206 fclose( ioBUG );
2207 }
2208
2209 /* sanity check: see whether this model has the correct effective temperature */
2210 if( ! lgValidModel( rfield.tNu[rfield.nShape], rfield.tslop[rfield.nShape], aval[0], 0.05 ) )
2211 TotalInsanity();
2212
2213 /* set limits for optimizer */
2214 SetLimits( grid, lval[0], NULL, NULL, useTr, ValTr, val0_lo, val0_hi );
2215
2216 /* now write some final info */
2217 if( called.lgTalk )
2218 {
2219 fprintf( ioQQQ, " * c<< FINAL: T_eff = %7.1f, ", aval[0] );
2220 fprintf( ioQQQ, "log(g) = %4.2f, M(ZAMS) = %5.1f, age = ", aval[1], pow(10.,aval[2]) );
2221 fprintf( ioQQQ, PrintEfmt("%8.2e",aval[3]) );
2222 fprintf( ioQQQ, " >>> *\n" );
2223 }
2224
2225 FREE_CHECK( indhiTr );
2226 FREE_CHECK( indloTr );
2227 FREE_CHECK( ValTr );
2228 return;
2229}
2230
2231/* find which models to use for interpolation along a given evolutionary track */
2233 long track,
2234 double par2, /* requested log(g) or age */
2235 long off, /* determines which parameter to match 0 -> log(g), 2 -> age */
2236 realnum *ValTr,/* ValTr[track]: Teff/log(M) value for interpolated model along track */
2237 long *indloTr, /* indloTr[track]: model number for first model used in interpolation */
2238 long *indhiTr) /* indhiTr[track]: model number for second model used in interpolation */
2239{
2240 long index[2], j, mod1, mod2;
2241
2242 DEBUG_ENTRY( "FindHCoStar()" );
2243
2244 indloTr[track] = -2;
2245 indhiTr[track] = -2;
2246 ValTr[track] = -FLT_MAX;
2247
2248 index[1] = track;
2249
2250 for( j=0; j < grid->trackLen[track]; j++ )
2251 {
2252 index[0] = j;
2253 mod1 = grid->jval[JIndex(grid,index)];
2254
2255 /* do we have an exact match ? */
2256 if( fabs(par2-grid->telg[mod1].par[off+1]) <= 10.*FLT_EPSILON*fabs(grid->telg[mod1].par[off+1]) )
2257 {
2258 indloTr[track] = mod1;
2259 indhiTr[track] = mod1;
2260 ValTr[track] = (realnum)grid->telg[mod1].par[off];
2261 return;
2262 }
2263 }
2264
2265 for( j=0; j < grid->trackLen[track]-1; j++ )
2266 {
2267 index[0] = j;
2268 mod1 = grid->jval[JIndex(grid,index)];
2269 index[0] = j+1;
2270 mod2 = grid->jval[JIndex(grid,index)];
2271
2272 /* do we interpolate ? */
2273 if( (par2 - grid->telg[mod1].par[off+1])*(par2 - grid->telg[mod2].par[off+1]) < 0. )
2274 {
2275 double frac;
2276
2277 indloTr[track] = mod1;
2278 indhiTr[track] = mod2;
2279 frac = (par2 - grid->telg[mod2].par[off+1])/
2280 (grid->telg[mod1].par[off+1] - grid->telg[mod2].par[off+1]);
2281 ValTr[track] = (realnum)(frac*grid->telg[mod1].par[off] +
2282 (1.-frac)*grid->telg[mod2].par[off] );
2283 break;
2284 }
2285 }
2286 return;
2287}
2288
2289/* find which tracks to use for interpolation in between tracks */
2291 double par1, /* requested Teff or ZAMS mass */
2292 realnum *ValTr, /* internal workspace */
2293 long useTr[]) /* useTr[0]: track number for first track to be used in interpolation
2294 * (i.e., 0 means 'A', etc.)
2295 * useTr[1]: track number for second track to be used in interpolation
2296 * NOTE: FindVCoStar raises a flag when interpolating between non-adjoining
2297 * tracks, i.e. when (useTr[1]-useTr[0]) > 1 */
2298{
2299 long j;
2300
2301 DEBUG_ENTRY( "FindVCoStar()" );
2302
2303 useTr[0] = -1;
2304 useTr[1] = -1;
2305
2306 for( j=0; j < grid->nTracks; j++ )
2307 {
2308 /* do we have an exact match ? */
2309 if( ValTr[j] != -FLT_MAX && fabs(par1-(double)ValTr[j]) <= 10.*FLT_EPSILON*fabs(ValTr[j]) )
2310 {
2311 useTr[0] = j;
2312 useTr[1] = j;
2313 break;
2314 }
2315 }
2316
2317 if( useTr[0] >= 0 )
2318 {
2319 return;
2320 }
2321
2322 for( j=0; j < grid->nTracks-1; j++ )
2323 {
2324 if( ValTr[j] != -FLT_MAX )
2325 {
2326 long int i,j2;
2327
2328 /* find next valid track */
2329 j2 = 0;
2330 for( i = j+1; i < grid->nTracks; i++ )
2331 {
2332 if( ValTr[i] != -FLT_MAX )
2333 {
2334 j2 = i;
2335 break;
2336 }
2337 }
2338
2339 /* do we interpolate ? */
2340 if( j2 > 0 && ((realnum)par1-ValTr[j])*((realnum)par1-ValTr[j2]) < 0.f )
2341 {
2342 useTr[0] = j;
2343 useTr[1] = j2;
2344 break;
2345 }
2346 }
2347 }
2348
2349 /* raise caution when we interpolate between non-adjoining tracks */
2350 continuum.lgCoStarInterpolationCaution = ( useTr[1]-useTr[0] > 1 );
2351 return;
2352}
2353
2354/* Make a listing of all the models in the CoStar grid */
2356{
2357 long index[2], maxlen, n;
2358
2359 DEBUG_ENTRY( "CoStarListModels()" );
2360
2361 maxlen = 0;
2362 for( n=0; n < grid->nTracks; n++ )
2363 maxlen = MAX2( maxlen, grid->trackLen[n] );
2364
2365 fprintf( ioQQQ, "\n" );
2366 fprintf( ioQQQ, " Track\\Index |" );
2367 for( n = 0; n < maxlen; n++ )
2368 fprintf( ioQQQ, " %5ld ", n+1 );
2369 fprintf( ioQQQ, "\n" );
2370 fprintf( ioQQQ, "--------------|" );
2371 for( n = 0; n < maxlen; n++ )
2372 fprintf( ioQQQ, "----------------" );
2373 fprintf( ioQQQ, "\n" );
2374
2375 for( index[1]=0; index[1] < grid->nTracks; ++index[1] )
2376 {
2377 long ptr;
2378 double Teff, alogg, Mass;
2379
2380 fprintf( ioQQQ, " %c", (char)('A'+index[1]) );
2381 index[0] = 0;
2382 ptr = grid->jval[JIndex(grid,index)];
2383 Mass = pow(10.,grid->telg[ptr].par[2]);
2384 fprintf( ioQQQ, " (%3.0f Msol) |", Mass );
2385
2386 for( index[0]=0; index[0] < grid->trackLen[index[1]]; ++index[0] )
2387 {
2388 ptr = grid->jval[JIndex(grid,index)];
2389 Teff = grid->telg[ptr].par[0];
2390 alogg = grid->telg[ptr].par[1];
2391 fprintf( ioQQQ, " (%6.1f,%4.2f)", Teff, alogg );
2392 }
2393 fprintf( ioQQQ, "\n" );
2394 }
2395 return;
2396}
2397
2398/* RauchInitializeSub does the actual work of preparing the ascii file */
2399STATIC int RauchInitializeSub(const char chFName[],
2400 const char chSuff[],
2401 const vector<mpp>& telg,
2402 long nmods,
2403 long n,
2404 long ngrids,
2405 const double par2[], /* par2[ngrids] */
2406 int format)
2407{
2408 char chLine[INPUT_LINE_LENGTH]; /* used for getting input lines */
2409
2410 FILE *ioOut, /* pointer to output file we came here to create*/
2411 *ioIn; /* pointer to input files we will read */
2412
2413 long int i, j;
2414
2415 double *wavl, *fluxes;
2416
2417 DEBUG_ENTRY( "RauchInitializeSub()" );
2418
2419 /* grab some space for the wavelengths and fluxes */
2420 wavl = (double *)MALLOC( sizeof(double)*NRAUCH);
2421 fluxes = (double *)MALLOC( sizeof(double)*NRAUCH);
2422
2423 try
2424 {
2425 if( n == 1 )
2426 ioOut = open_data( chFName, "w", AS_LOCAL_ONLY );
2427 else
2428 ioOut = open_data( chFName, "a", AS_LOCAL_ONLY );
2429 }
2430 catch( cloudy_exit )
2431 {
2432 goto error;
2433 }
2434
2435 if( n == 1 )
2436 {
2437 fprintf( ioOut, " %ld\n", VERSION_ASCII );
2438 fprintf( ioOut, " %d\n", ( ngrids == 1 ? 2 : 3 ) );
2439 fprintf( ioOut, " %d\n", ( ngrids == 1 ? 2 : 3 ) );
2440 fprintf( ioOut, " Teff\n" );
2441 fprintf( ioOut, " log(g)\n" );
2442 if( ngrids == 2 )
2443 fprintf( ioOut, " log(Z)\n" );
2444 else if( ngrids == 11 )
2445 fprintf( ioOut, " f(He)\n" );
2446 /* NB - this is based on the assumption that each of the planes in the cubic grid is the same */
2447 fprintf( ioOut, " %ld\n", nmods*ngrids );
2448 fprintf( ioOut, " %d\n", NRAUCH );
2449 /* Rauch models give the wavelength in Angstrom */
2450 fprintf( ioOut, " lambda\n" );
2451 /* conversion factor to Angstrom */
2452 fprintf( ioOut, " %.8e\n", 1. );
2453 /* Rauch models give the "Astrophysical" flux F_lambda in erg/cm^2/s/cm */
2454 fprintf( ioOut, " F_lambda\n" );
2455 /* the factor PI*1e-8 is needed to convert to "regular" flux in erg/cm^2/s/Angstrom */
2456 fprintf( ioOut, " %.8e\n", PI*1.e-8 );
2457 /* NB - this is based on the assumption that each of the planes in the cubic grid is the same */
2458 for( j=0; j < ngrids; j++ )
2459 {
2460 /* write out the {Teff,log(g)} grid */
2461 for( i=0; i < nmods; i++ )
2462 {
2463 if( ngrids == 1 )
2464 fprintf( ioOut, " %.0f %.1f", telg[i].par[0], telg[i].par[1] );
2465 else
2466 fprintf( ioOut, " %.0f %.1f %.1f", telg[i].par[0], telg[i].par[1], par2[j] );
2467 if( ((i+1)%4) == 0 )
2468 fprintf( ioOut, "\n" );
2469 }
2470 if( (i%4) != 0 )
2471 fprintf( ioOut, "\n" );
2472 }
2473
2474 fprintf( ioQQQ, " Writing: " );
2475 }
2476
2477 for( i=0; i < nmods; i++ )
2478 {
2479 /* must create name of next stellar atmosphere */
2480 if( format == 1 )
2481 sprintf( chLine, "%7.7ld_%2ld", (long)(telg[i].par[0]+0.5), (long)(10.*telg[i].par[1]+0.5) );
2482 else if( format == 2 )
2483 sprintf( chLine, "%7.7ld_%.2f", (long)(telg[i].par[0]+0.5), telg[i].par[1] );
2484 else
2485 {
2486 fprintf( ioQQQ, " insanity in RauchInitializeSub\n" );
2487 ShowMe();
2489 }
2490 string chFileName( chLine );
2491 chFileName += chSuff;
2492 /* now open next stellar atmosphere for reading*/
2493 try
2494 {
2495 ioIn = open_data( chFileName.c_str(), "r", AS_LOCAL_ONLY );
2496 }
2497 catch( cloudy_exit )
2498 {
2499 goto error;
2500 }
2501
2502 /* get first line */
2503 j = 0;
2504 if( read_whole_line( chLine, (int)sizeof(chLine), ioIn ) == NULL )
2505 {
2506 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2507 i, j );
2508 goto error;
2509 }
2510 /* >>chng 02 nov 20, now keep reading them until don't hit the *
2511 * since number of comments may change */
2512 while( chLine[0] == '*' )
2513 {
2514 if( read_whole_line( chLine, (int)sizeof(chLine), ioIn ) == NULL )
2515 {
2516 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2517 i, j );
2518 goto error;
2519 }
2520 ++j;
2521 }
2522
2523 for( j=0; j < NRAUCH; j++ )
2524 {
2525 double ttemp, wl;
2526 /* get the input line */
2527 /* >>chng 02 nov 20, don't reread very first line image since we got it above */
2528 if( j > 0 )
2529 {
2530 if(read_whole_line( chLine, (int)sizeof(chLine), ioIn )==NULL )
2531 {
2532 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2533 i, j );
2534 goto error;
2535 }
2536 }
2537
2538 /* scan off wavelength and flux)*/
2539 if( sscanf( chLine, "%lf %le", &wl, &ttemp ) != 2 )
2540 {
2541 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2542 i, j );
2543 goto error;
2544 }
2545
2546 if( i == 0 )
2547 wavl[j] = wl;
2548 else
2549 {
2550 /* check if this model is on the same wavelength grid as the first */
2551 if( !fp_equal(wavl[j],wl,10) )
2552 {
2553 fprintf( ioQQQ, " RauchInitializeSub error in atmosphere file %4ld%4ld\n",
2554 i, j );
2555 goto error;
2556 }
2557 }
2558 fluxes[j] = ttemp;
2559 }
2560
2561 /* finished - close the unit */
2562 fclose(ioIn);
2563
2564 /* now write to output file */
2565 if( i == 0 && n == 1 )
2566 {
2567 /* wavelength grid is the same for all models, so write only once */
2568 for( j=0; j < NRAUCH; j++ )
2569 {
2570 fprintf( ioOut, " %.4e", wavl[j] );
2571 /* want to have 5 numbers per line */
2572 if( ((j+1)%5) == 0 )
2573 fprintf( ioOut, "\n" );
2574 }
2575 /* need to throw a newline if we did not end on an exact line */
2576 if( (j%5) != 0 )
2577 fprintf( ioOut, "\n" );
2578 }
2579
2580 for( j=0; j < NRAUCH; j++ )
2581 {
2582 fprintf( ioOut, " %.4e", fluxes[j] );
2583 /* want to have 5 numbers per line */
2584 if( ((j+1)%5) == 0 )
2585 fprintf( ioOut, "\n" );
2586 }
2587 /* need to throw a newline if we did not end on an exact line */
2588 if( (j%5) != 0 )
2589 fprintf( ioOut, "\n" );
2590
2591 /* print to screen to show progress */
2592 fprintf( ioQQQ, "." );
2593 fflush( ioQQQ );
2594 }
2595
2596 if( n == ngrids )
2597 fprintf( ioQQQ, " done.\n" );
2598
2599 fclose(ioOut);
2600
2601 /* free the space grabbed above */
2602 FREE_CHECK( fluxes );
2603 FREE_CHECK( wavl );
2604 return 0;
2605
2606error:
2607 FREE_CHECK( fluxes );
2608 FREE_CHECK( wavl );
2609 return 1;
2610}
2611
2612STATIC void RauchReadMPP(vector<mpp>& telg1,
2613 vector<mpp>& telg2,
2614 vector<mpp>& telg3,
2615 vector<mpp>& telg4,
2616 vector<mpp>& telg5,
2617 vector<mpp>& telg6)
2618{
2619 DEBUG_ENTRY( "RauchReadMPP()" );
2620
2621 const char fnam[] = "rauch_models.dat";
2622 fstream ioDATA;
2623 open_data( ioDATA, fnam, mode_r );
2624
2625 string line;
2626 getdataline( ioDATA, line );
2627 long version;
2628 istringstream iss( line );
2629 iss >> version;
2630 if( version != VERSION_RAUCH_MPP )
2631 {
2632 fprintf( ioQQQ, " RauchReadMPP: the version of %s is not the current version.\n", fnam );
2633 fprintf( ioQQQ, " Please obtain the current version from the Cloudy web site.\n" );
2634 fprintf( ioQQQ, " I expected to find version %ld and got %ld instead.\n",
2635 VERSION_RAUCH_MPP, version );
2637 }
2638
2639 getdataline( ioDATA, line );
2640 unsigned long ndata;
2641 istringstream iss2( line );
2642 iss2 >> ndata;
2643 ASSERT( ndata == telg1.size() );
2644 // this implicitly assumes there is exactly one comment line between
2645 // the number of data points and the start of the data
2646 getline( ioDATA, line );
2647 // read data for H-Ca grid
2648 for( unsigned long i=0; i < ndata; ++i )
2649 ioDATA >> telg1[i].par[0] >> telg1[i].par[1];
2650 getline( ioDATA, line );
2651
2652 getdataline( ioDATA, line );
2653 istringstream iss3( line );
2654 iss3 >> ndata;
2655 ASSERT( ndata == telg2.size() );
2656 getline( ioDATA, line );
2657 // read data for H-Ni grid
2658 for( unsigned long i=0; i < ndata; ++i )
2659 ioDATA >> telg2[i].par[0] >> telg2[i].par[1];
2660 getline( ioDATA, line );
2661
2662 getdataline( ioDATA, line );
2663 istringstream iss4( line );
2664 iss4 >> ndata;
2665 ASSERT( ndata == telg3.size() );
2666 getline( ioDATA, line );
2667 // read data for PG1159 grid
2668 for( unsigned long i=0; i < ndata; ++i )
2669 ioDATA >> telg3[i].par[0] >> telg3[i].par[1];
2670 getline( ioDATA, line );
2671
2672 getdataline( ioDATA, line );
2673 istringstream iss5( line );
2674 iss5 >> ndata;
2675 ASSERT( ndata == telg4.size() );
2676 getline( ioDATA, line );
2677 // read data for pure H grid
2678 for( unsigned long i=0; i < ndata; ++i )
2679 ioDATA >> telg4[i].par[0] >> telg4[i].par[1];
2680 getline( ioDATA, line );
2681
2682 getdataline( ioDATA, line );
2683 istringstream iss6( line );
2684 iss6 >> ndata;
2685 ASSERT( ndata == telg5.size() );
2686 getline( ioDATA, line );
2687 // read data for pure He grid
2688 for( unsigned long i=0; i < ndata; ++i )
2689 ioDATA >> telg5[i].par[0] >> telg5[i].par[1];
2690 getline( ioDATA, line );
2691
2692 getdataline( ioDATA, line );
2693 istringstream iss7( line );
2694 iss7 >> ndata;
2695 ASSERT( ndata == telg6.size() );
2696 getline( ioDATA, line );
2697 // read data for pure H+He grid
2698 for( unsigned long i=0; i < ndata; ++i )
2699 ioDATA >> telg6[i].par[0] >> telg6[i].par[1];
2700 getline( ioDATA, line );
2701
2702 getdataline( ioDATA, line );
2703 istringstream iss8( line );
2704 iss8 >> version;
2705 ASSERT( version == VERSION_RAUCH_MPP );
2706
2707 return;
2708}
2709
2710inline void getdataline(fstream& ioDATA,
2711 string& line)
2712{
2713 do
2714 {
2715 getline( ioDATA, line );
2716 }
2717 while( line[0] == '#' );
2718 return;
2719}
2720
2721/* lgCompileAtmosphere does the actual rebinning onto the Cloudy grid and writes the binary file */
2722/* >>chng 01 feb 12, added return value to indicate success (0) or failure (1) */
2723STATIC bool lgCompileAtmosphere(const char chFNameIn[],
2724 const char chFNameOut[],
2725 const realnum Edges[], /* Edges[nedges] */
2726 long nedges,
2727 process_counter& pc)
2728{
2729 FILE *ioIN; /* used for input */
2730 FILE *ioOUT; /* used for output */
2731
2732 char chDataType[11];
2733 char names[MDIM][MNAM+1];
2734
2735 bool lgFreqX, lgFreqY, lgFlip;
2736 int32 val[7];
2737 uint32 uval[2];
2738 double dval[3];
2739 char md5sum[NMD5];
2740 long int i, imod, version, nd, ndim, npar, nmods, ngrid;
2741
2742 /* these will be malloced into large work arrays */
2743 realnum *StarEner = NULL, *StarFlux = NULL, *CloudyFlux = NULL, *scratch = NULL;
2744 vector<realnum> SaveAnu(rfield.nupper);
2745
2746 double convert_wavl, convert_flux;
2747
2748 mpp *telg = NULL;
2749
2750 DEBUG_ENTRY( "lgCompileAtmosphere()" );
2751
2752 try
2753 {
2754 ioIN = open_data( chFNameIn, "r", AS_LOCAL_ONLY );
2755 }
2756 catch( cloudy_exit )
2757 {
2758 goto error;
2759 }
2760 fprintf( ioQQQ, " lgCompileAtmosphere got %s.\n", chFNameIn );
2761
2762 /* read version number */
2763 if( fscanf( ioIN, "%ld", &version ) != 1 )
2764 {
2765 fprintf( ioQQQ, " lgCompileAtmosphere failed reading VERSION.\n" );
2766 goto error;
2767 }
2768
2769 if( version != VERSION_ASCII )
2770 {
2771 fprintf( ioQQQ, " lgCompileAtmosphere: there is a version number mismatch in"
2772 " the ascii atmosphere file: %s.\n", chFNameIn );
2773 fprintf( ioQQQ, " lgCompileAtmosphere: Please recreate this file or download the"
2774 " latest version following the instructions on the Cloudy website.\n" );
2775 goto error;
2776 }
2777
2778 /* >>chng 06 jun 10, read the dimension of the grid, PvH */
2779 if( fscanf( ioIN, "%ld", &ndim ) != 1 || ndim <= 0 || ndim > MDIM )
2780 {
2781 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid dimension of grid.\n" );
2782 goto error;
2783 }
2784
2785 /* >>chng 06 jun 12, read the number of model parameters, PvH */
2786 if( fscanf( ioIN, "%ld", &npar ) != 1 || npar <= 0 || npar < ndim || npar > MDIM )
2787 {
2788 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid no. of model parameters.\n" );
2789 goto error;
2790 }
2791
2792 /* make sure valgrind doesn't trip on the binary write of this array */
2793 memset( names, '\0', MDIM*(MNAM+1) );
2794
2795 for( nd=0; nd < npar; nd++ )
2796 {
2797 if( fscanf( ioIN, "%6s", names[nd] ) != 1 )
2798 {
2799 fprintf( ioQQQ, " lgCompileAtmosphere failed reading parameter label.\n" );
2800 goto error;
2801 }
2802 }
2803
2804 /* >>chng 05 nov 18, read the following extra parameters from the ascii file, PvH */
2805 if( fscanf( ioIN, "%ld", &nmods ) != 1 || nmods <= 0 )
2806 {
2807 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid number of models.\n" );
2808 goto error;
2809 }
2810
2811 if( fscanf( ioIN, "%ld", &ngrid ) != 1 || ngrid <= 1 )
2812 {
2813 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid number of grid points.\n" );
2814 goto error;
2815 }
2816
2817 /* read data type for wavelengths, allowed values are lambda, nu */
2818 if( fscanf( ioIN, "%10s", chDataType ) != 1 )
2819 {
2820 fprintf( ioQQQ, " lgCompileAtmosphere failed reading wavl DataType string.\n" );
2821 goto error;
2822 }
2823
2824 if( strcmp( chDataType, "lambda" ) == 0 )
2825 lgFreqX = false;
2826 else if( strcmp( chDataType, "nu" ) == 0 )
2827 lgFreqX = true;
2828 else {
2829 fprintf( ioQQQ, " lgCompileAtmosphere found illegal wavl DataType: %s.\n", chDataType );
2830 goto error;
2831 }
2832
2833 if( fscanf( ioIN, "%le", &convert_wavl ) != 1 || convert_wavl <= 0. )
2834 {
2835 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid wavl conversion factor.\n" );
2836 goto error;
2837 }
2838
2839 /* read data type for flux, allowed values F_lambda, H_lambda, F_nu, H_nu */
2840 if( fscanf( ioIN, "%10s", chDataType ) != 1 )
2841 {
2842 fprintf( ioQQQ, " lgCompileAtmosphere failed reading flux DataType string.\n" );
2843 goto error;
2844 }
2845
2846 if( strcmp( chDataType, "F_lambda" ) == 0 || strcmp( chDataType, "H_lambda" ) == 0 )
2847 lgFreqY = false;
2848 else if( strcmp( chDataType, "F_nu" ) == 0 || strcmp( chDataType, "H_nu" ) == 0 )
2849 lgFreqY = true;
2850 else {
2851 fprintf( ioQQQ, " lgCompileAtmosphere found illegal flux DataType: %s.\n", chDataType );
2852 goto error;
2853 }
2854
2855 if( fscanf( ioIN, "%le", &convert_flux ) != 1 || convert_flux <= 0. )
2856 {
2857 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid flux conversion factor.\n" );
2858 goto error;
2859 }
2860
2861 telg = (mpp *)CALLOC( (size_t)nmods, sizeof(mpp) );
2862
2863 for( i=0; i < nmods; i++ )
2864 {
2865 for( nd=0; nd < npar; nd++ )
2866 {
2867 if( fscanf( ioIN, "%le", &telg[i].par[nd] ) != 1 )
2868 {
2869 fprintf( ioQQQ, " lgCompileAtmosphere failed reading valid model parameter.\n" );
2870 goto error;
2871 }
2872 }
2873 }
2874
2875 try
2876 {
2877 ioOUT = open_data( chFNameOut, "wb", AS_LOCAL_ONLY );
2878 }
2879 catch( cloudy_exit )
2880 {
2881 goto error;
2882 }
2883
2884 val[0] = (int32)VERSION_BIN;
2885 val[1] = (int32)MDIM;
2886 val[2] = (int32)MNAM;
2887 val[3] = (int32)ndim;
2888 val[4] = (int32)npar;
2889 val[5] = (int32)nmods;
2890 val[6] = (int32)rfield.nupper;
2891 uval[0] = sizeof(val) + sizeof(uval) + sizeof(dval) + sizeof(md5sum) +
2892 sizeof(names) + nmods*sizeof(mpp); /* nOffset */
2893 uval[1] = rfield.nupper*sizeof(realnum); /* nBlocksize */
2894 dval[0] = double(rfield.emm);
2895 dval[1] = double(rfield.egamry);
2896 dval[2] = double(continuum.ResolutionScaleFactor);
2897
2898 strncpy( md5sum, continuum.mesh_md5sum.c_str(), sizeof(md5sum) );
2899
2900 for (long i=0; i<rfield.nupper; ++i)
2901 SaveAnu[i] = (realnum) rfield.AnuOrg[i];
2902
2903 if( fwrite( val, sizeof(val), 1, ioOUT ) != 1 ||
2904 fwrite( uval, sizeof(uval), 1, ioOUT ) != 1 ||
2905 /* write out the lower, upper bound of the energy mesh, and the res scale factor */
2906 fwrite( dval, sizeof(dval), 1, ioOUT ) != 1 ||
2907 /* write out the (modified) md5 checksum of continuum_mesh.ini */
2908 fwrite( md5sum, sizeof(md5sum), 1, ioOUT ) != 1 ||
2909 fwrite( names, sizeof(names), 1, ioOUT ) != 1 ||
2910 /* write out the array of {Teff,log(g)} pairs */
2911 fwrite( telg, sizeof(mpp), (size_t)nmods, ioOUT ) != (size_t)nmods ||
2912 /* write out the cloudy energy grid for later sanity checks */
2913 fwrite( get_ptr(SaveAnu), (size_t)uval[1], 1, ioOUT ) != 1 )
2914 {
2915 fprintf( ioQQQ, " lgCompileAtmosphere failed writing header of output file.\n" );
2916 goto error;
2917 }
2918
2919 /* MALLOC some workspace */
2920 StarEner = (realnum *)MALLOC( sizeof(realnum)*ngrid );
2921 scratch = (realnum *)MALLOC( sizeof(realnum)*ngrid );
2922 StarFlux = (realnum *)MALLOC( sizeof(realnum)*ngrid );
2923 CloudyFlux = (realnum *)MALLOC( (size_t)uval[1] );
2924
2925 /* read wavelength grid */
2926 for( i=0; i < ngrid; i++ )
2927 {
2928 double help;
2929 if( fscanf( ioIN, "%lg", &help ) != 1 )
2930 {
2931 fprintf( ioQQQ, " lgCompileAtmosphere failed reading wavelength.\n" );
2932 goto error;
2933 }
2934 /* this conversion makes sure that scratch[i] is
2935 * either wavelength in Angstrom or frequency in Hz */
2936 scratch[i] = (realnum)(help*convert_wavl);
2937
2938 if( scratch[i] <= 0.f )
2939 {
2940 fprintf( ioQQQ, " PROBLEM: a non-positive %s was found, value: %e\n",
2941 lgFreqX ? "frequency" : "wavelength", help );
2943 }
2944 }
2945
2946 lgFlip = ( !lgFreqX && scratch[0] < scratch[1] ) || ( lgFreqX && scratch[0] > scratch[1] );
2947
2948 /* convert continuum over to increasing frequency in Ryd */
2949 for( i=0; i < ngrid; i++ )
2950 {
2951 /* convert scratch[i] to frequency in Ryd */
2952 if( lgFreqX )
2953 scratch[i] /= (realnum)FR1RYD;
2954 else
2955 scratch[i] = (realnum)(RYDLAM/scratch[i]);
2956
2957 if( lgFlip )
2958 StarEner[ngrid-i-1] = scratch[i];
2959 else
2960 StarEner[i] = scratch[i];
2961 }
2962
2963 ASSERT( StarEner[0] > 0.f );
2964 /* make sure the array is in ascending order */
2965 for( i=1; i < ngrid; i++ )
2966 {
2967 if( StarEner[i] <= StarEner[i-1] )
2968 {
2969 fprintf( ioQQQ, " PROBLEM: the %s grid is not strictly monotonically increasing/decreasing\n",
2970 lgFreqX ? "frequency" : "wavelength" );
2972 }
2973 }
2974
2975 fprintf( ioQQQ, " Compiling: " );
2976
2977 for( imod=0; imod < nmods; imod++ )
2978 {
2979 const realnum CONVERT_FNU = (realnum)(1.e8*SPEEDLIGHT/POW2(FR1RYD));
2980
2981 /* now read the stellar fluxes */
2982 for( i=0; i < ngrid; i++ )
2983 {
2984 double help;
2985 if( fscanf( ioIN, "%lg", &help ) != 1 )
2986 {
2987 fprintf( ioQQQ, " lgCompileAtmosphere failed reading star flux.\n" );
2988 goto error;
2989 }
2990 /* this conversion makes sure that scratch[i] is either
2991 * F_nu in erg/cm^2/s/Hz or F_lambda in erg/cm^2/s/A */
2992 scratch[i] = (realnum)(help*convert_flux);
2993
2994 /* this can underflow on the Wien tail */
2995 if( scratch[i] < 0.f )
2996 {
2997 fprintf( ioQQQ, "\n PROBLEM: a negative flux was found, model number %ld, value: %e\n",
2998 imod+1, help );
3000 }
3001 }
3002
3003 for( i=0; i < ngrid; i++ )
3004 {
3005 if( lgFlip )
3006 StarFlux[ngrid-i-1] = scratch[i];
3007 else
3008 StarFlux[i] = scratch[i];
3009 }
3010
3011 for( i=0; i < ngrid; i++ )
3012 {
3013 /* this converts to F_nu in erg/cm^2/s/Hz */
3014 if( !lgFreqY )
3015 StarFlux[i] *= CONVERT_FNU/POW2(StarEner[i]);
3016 ASSERT( StarFlux[i] >= 0.f );
3017 }
3018
3019 /* the re-binned values are returned in the "CloudyFlux" array */
3020 RebinAtmosphere( ngrid, StarEner, StarFlux, CloudyFlux, nedges, Edges );
3021
3022 /* write the continuum out as a binary file */
3023 if( fwrite( CloudyFlux, (size_t)uval[1], 1, ioOUT ) != 1 )
3024 {
3025 fprintf( ioQQQ, " lgCompileAtmosphere failed writing star flux.\n" );
3026 goto error;
3027 }
3028
3029 fprintf( ioQQQ, "." );
3030 fflush( ioQQQ );
3031 }
3032
3033 fprintf( ioQQQ, " done.\n" );
3034
3035 fclose(ioIN);
3036 fclose(ioOUT);
3037
3038 /* now free up the memory we claimed */
3039 FREE_CHECK( CloudyFlux );
3040 FREE_CHECK( StarFlux );
3041 FREE_CHECK( StarEner );
3042 FREE_CHECK( scratch );
3043 FREE_CHECK( telg );
3044
3045 fprintf( ioQQQ, " lgCompileAtmosphere completed ok.\n\n" );
3046 ++pc.nOK;
3047 return false;
3048
3049error:
3050 FREE_SAFE( CloudyFlux );
3051 FREE_SAFE( StarFlux );
3052 FREE_SAFE( StarEner );
3053 FREE_SAFE( scratch );
3054 FREE_SAFE( telg );
3055 ++pc.nFail;
3056 return true;
3057}
3058
3060 bool lgList)
3061{
3062 long nd;
3063 int32 version, mdim, mnam;
3064 double mesh_elo, mesh_ehi;
3065 char md5sum[NMD5];
3066
3067 DEBUG_ENTRY( "InitGrid()" );
3068
3069 try
3070 {
3071 grid->ioIN = open_data( grid->name.c_str(), "rb", grid->scheme );
3072 }
3073 catch( cloudy_exit )
3074 {
3075 /* something went wrong */
3076 /* NB NB - DO NOT CHANGE THE FOLLOWING ERROR MESSAGE! checkall.pl picks it up */
3077 fprintf( ioQQQ, " Error: stellar atmosphere file not found.\n" );
3078 fprintf(ioQQQ , "\n\n If the path is set then it is possible that the stellar atmosphere data files do not exist.\n");
3079 fprintf(ioQQQ , " Have the stellar data files been downloaded and compiled with the COMPILE STARS command?\n");
3080 fprintf(ioQQQ , " If you are simply running the test suite and do not need the stellar continua then you should simply ignore this failure\n");
3082 }
3083
3084 /* >>chng 01 oct 17, add version and size to this array */
3085 if( fread( &version, sizeof(version), 1, grid->ioIN ) != 1 ||
3086 fread( &mdim, sizeof(mdim), 1, grid->ioIN ) != 1 ||
3087 fread( &mnam, sizeof(mnam), 1, grid->ioIN ) != 1 ||
3088 fread( &grid->ndim, sizeof(grid->ndim), 1, grid->ioIN ) != 1 ||
3089 fread( &grid->npar, sizeof(grid->npar), 1, grid->ioIN ) != 1 ||
3090 fread( &grid->nmods, sizeof(grid->nmods), 1, grid->ioIN ) != 1 ||
3091 fread( &grid->ngrid, sizeof(grid->ngrid), 1, grid->ioIN ) != 1 ||
3092 fread( &grid->nOffset, sizeof(grid->nOffset), 1, grid->ioIN ) != 1 ||
3093 fread( &grid->nBlocksize, sizeof(grid->nBlocksize), 1, grid->ioIN ) != 1 ||
3094 fread( &mesh_elo, sizeof(mesh_elo), 1, grid->ioIN ) != 1 ||
3095 fread( &mesh_ehi, sizeof(mesh_ehi), 1, grid->ioIN ) != 1 ||
3096 fread( &rfield.RSFCheck[rfield.nShape], sizeof(rfield.RSFCheck[rfield.nShape]), 1, grid->ioIN ) != 1 ||
3097 fread( md5sum, sizeof(md5sum), 1, grid->ioIN ) != 1 )
3098 {
3099 fprintf( ioQQQ, " InitGrid failed reading header.\n" );
3101 }
3102
3103 /* do some sanity checks */
3104 if( version != VERSION_BIN )
3105 {
3106 fprintf( ioQQQ, " InitGrid: there is a version mismatch between"
3107 " the compiled atmospheres file I expected and the one I found.\n" );
3108 fprintf( ioQQQ, " InitGrid: Please recompile the stellar"
3109 " atmospheres file with the command: %s.\n", grid->command );
3111 }
3112
3113 if( mdim != MDIM || mnam != MNAM )
3114 {
3115 fprintf( ioQQQ, " InitGrid: the compiled atmospheres file is produced"
3116 " with an incompatible version of Cloudy.\n" );
3117 fprintf( ioQQQ, " InitGrid: Please recompile the stellar"
3118 " atmospheres file with the command: %s.\n", grid->command );
3120 }
3121
3122 if( !fp_equal( double(rfield.emm), mesh_elo ) ||
3123 !fp_equal( double(rfield.egamry), mesh_ehi ) ||
3124 strncmp( continuum.mesh_md5sum.c_str(), md5sum, NMD5 ) != 0 )
3125 {
3126 fprintf( ioQQQ, " InitGrid: the compiled atmospheres file is produced"
3127 " with an incompatible frequency grid.\n" );
3128 fprintf( ioQQQ, " InitGrid: Please recompile the stellar"
3129 " atmospheres file with the command: %s.\n", grid->command );
3131 }
3132
3133 ASSERT( grid->ndim > 0 && grid->ndim <= MDIM );
3134 ASSERT( grid->npar >= grid->ndim && grid->npar <= MDIM );
3135 ASSERT( grid->nmods > 0 );
3136 ASSERT( grid->ngrid > 0 );
3137 ASSERT( grid->nOffset > 0 );
3138 ASSERT( grid->nBlocksize > 0 );
3139
3140 rfield.nupper = grid->ngrid;
3141
3142 if( fread( &grid->names, sizeof(grid->names), 1, grid->ioIN ) != 1 )
3143 {
3144 fprintf( ioQQQ, " InitGrid failed reading names array.\n" );
3146 }
3147
3148 grid->lgIsTeffLoggGrid = ( grid->ndim >= 2 &&
3149 strcmp( grid->names[0], "Teff" ) == 0 &&
3150 strcmp( grid->names[1], "log(g)" ) == 0 );
3151
3152 grid->telg = (mpp *)MALLOC( sizeof(mpp)*grid->nmods );
3153 grid->val = (double **)MALLOC( sizeof(double*)*grid->ndim );
3154 for( nd=0; nd < grid->ndim; nd++ )
3155 {
3156 grid->val[nd] = (double *)MALLOC( sizeof(double)*grid->nmods );
3157 }
3158 grid->nval = (long *)MALLOC( sizeof(long)*grid->ndim );
3159
3160 if( fread( grid->telg, sizeof(mpp), grid->nmods, grid->ioIN ) != (size_t)grid->nmods )
3161 {
3162 fprintf( ioQQQ, " InitGrid failed reading model parameter block.\n" );
3164 }
3165
3166# ifdef SEEK_END
3167 /* sanity check: does the file have the correct length ? */
3168 /* NOTE: this operation is not necessarily supported by all operating systems
3169 * but if the preprocessor symbol SEEK_END exists it is assumed to be supported */
3170 int res = fseek( grid->ioIN, 0, SEEK_END );
3171 if( res == 0 )
3172 {
3173 long End = ftell( grid->ioIN );
3174 long Expected = grid->nOffset + (grid->nmods+1)*grid->nBlocksize;
3175 if( End != Expected )
3176 {
3177 fprintf( ioQQQ, " InitGrid: Problem performing sanity check for size of binary file.\n" );
3178 fprintf( ioQQQ, " InitGrid: I expected to find %ld bytes, but actually found %ld bytes.\n",
3179 Expected, End );
3180 fprintf( ioQQQ, " InitGrid: Please recompile the stellar"
3181 " atmospheres file with the command: %s.\n", grid->command );
3183 }
3184 }
3185# endif
3186
3187 InitIndexArrays( grid, lgList );
3188
3189 /* set default interpolation mode */
3190 grid->imode = IM_RECT_GRID;
3191 /* these are only used by CoStar grids */
3192 grid->trackLen = NULL;
3193 grid->nTracks = 0;
3194 grid->jval = NULL;
3195 return;
3196}
3197
3198/* check whether a binary atmosphere exists and is up-to-date */
3199STATIC bool lgValidBinFile(const char *binName, process_counter& pc, access_scheme scheme)
3200{
3201 int32 version, mdim, mnam;
3202 double mesh_elo, mesh_ehi, mesh_res_factor;
3203 char md5sum[NMD5];
3205
3206 DEBUG_ENTRY( "lgValidBinFile()" );
3207
3208 //
3209 // NB NB NB
3210 //
3211 // this routine is called when either of these two commands is issued:
3212 //
3213 // TABLE STAR AVAIL
3214 // COMPILE STAR [ additional parameters ]
3215 //
3216 // both these commands execute as soon as they are parsed and then terminate
3217 // hence it is safe to assume that no SET CONTINUUM RESOLUTION command will follow!
3218 //
3219 // !!! THE TEST BELOW FOR VALIDITY OF THE FILE DEPENDS ON THAT ASSUMPTION !!!
3220 //
3221
3222 grid.name = binName;
3223
3224 if( (grid.ioIN = open_data( grid.name.c_str(), "rb", scheme )) == NULL )
3225 return false;
3226
3227 if( fread( &version, sizeof(version), 1, grid.ioIN ) != 1 ||
3228 fread( &mdim, sizeof(mdim), 1, grid.ioIN ) != 1 ||
3229 fread( &mnam, sizeof(mnam), 1, grid.ioIN ) != 1 ||
3230 fread( &grid.ndim, sizeof(grid.ndim), 1, grid.ioIN ) != 1 ||
3231 fread( &grid.npar, sizeof(grid.npar), 1, grid.ioIN ) != 1 ||
3232 fread( &grid.nmods, sizeof(grid.nmods), 1, grid.ioIN ) != 1 ||
3233 fread( &grid.ngrid, sizeof(grid.ngrid), 1, grid.ioIN ) != 1 ||
3234 fread( &grid.nOffset, sizeof(grid.nOffset), 1, grid.ioIN ) != 1 ||
3235 fread( &grid.nBlocksize, sizeof(grid.nBlocksize), 1, grid.ioIN ) != 1 ||
3236 fread( &mesh_elo, sizeof(mesh_elo), 1, grid.ioIN ) != 1 ||
3237 fread( &mesh_ehi, sizeof(mesh_ehi), 1, grid.ioIN ) != 1 ||
3238 fread( &mesh_res_factor, sizeof(mesh_res_factor), 1, grid.ioIN ) != 1 ||
3239 fread( md5sum, sizeof(md5sum), 1, grid.ioIN ) != 1 )
3240 {
3241 fclose( grid.ioIN );
3242 return false;
3243 }
3244
3245 /* do some sanity checks */
3246 if( version != VERSION_BIN || mdim != MDIM || mnam != MNAM ||
3247 !fp_equal( double(rfield.emm), mesh_elo ) ||
3248 !fp_equal( double(rfield.egamry), mesh_ehi ) ||
3249 !fp_equal( double(continuum.ResolutionScaleFactor), mesh_res_factor ) ||
3250 strncmp( continuum.mesh_md5sum.c_str(), md5sum, NMD5 ) != 0 )
3251 {
3252 fclose( grid.ioIN );
3253 return false;
3254 }
3255
3256# ifdef SEEK_END
3257 /* sanity check: does the file have the correct length ? */
3258 /* NOTE: this operation is not necessarily supported by all operating systems
3259 * but if the preprocessor symbol SEEK_END exists it is assumed to be supported */
3260 int res = fseek( grid.ioIN, 0, SEEK_END );
3261 if( res == 0 )
3262 {
3263 long End = ftell( grid.ioIN );
3264 long Expected = grid.nOffset + (grid.nmods+1)*grid.nBlocksize;
3265 if( End != Expected )
3266 {
3267 fclose( grid.ioIN );
3268 return false;
3269 }
3270 }
3271# endif
3272
3273 fclose( grid.ioIN );
3274 ++pc.notProcessed; // the file is up-to-date -> no processing
3275 return true;
3276}
3277
3278/* check whether a ascii atmosphere file exists and is up-to-date */
3279STATIC bool lgValidAsciiFile(const char *ascName, access_scheme scheme)
3280{
3281 long version;
3282 FILE *ioIN;
3283
3284 DEBUG_ENTRY( "lgValidAsciiFile()" );
3285
3286 /* can we read the file? */
3287 if( (ioIN = open_data( ascName, "r", scheme )) == NULL )
3288 return false;
3289
3290 /* check version number */
3291 if( fscanf( ioIN, "%ld", &version ) != 1 || version != VERSION_ASCII )
3292 {
3293 fclose( ioIN );
3294 return false;
3295 }
3296
3297 fclose( ioIN );
3298 return true;
3299}
3300
3301/* sort CoStar models according to track and index number, store indices in grid->jval[] */
3302STATIC void InitGridCoStar(stellar_grid *grid) /* the grid parameters */
3303{
3304 char track;
3305 bool lgFound;
3306 long i, index[2];
3307
3308 DEBUG_ENTRY( "InitGridCoStar()" );
3309
3310 ASSERT( grid->ndim == 2 );
3311 ASSERT( grid->jlo != NULL && grid->jhi != NULL );
3312
3313 grid->jval = grid->jlo;
3314 FREE_CHECK( grid->jhi );
3315 grid->jlo = grid->jhi = NULL;
3316
3317 /* invalidate contents set by InitGrid first */
3318 memset( grid->jval, 0xff, (size_t)(grid->nval[0]*grid->nval[1]*sizeof(long)) );
3319
3320 grid->trackLen = (long *)CALLOC( (size_t)grid->nmods, sizeof(long) );
3321
3322 index[1] = 0;
3323 while( true )
3324 {
3325 index[0] = 0;
3326 track = (char)('A'+index[1]);
3327 do
3328 {
3329 lgFound = false;
3330 for( i=0; i < grid->nmods; i++ )
3331 {
3332 if( grid->telg[i].chGrid == track && grid->telg[i].modid == index[0]+1 )
3333 {
3334 grid->jval[JIndex(grid,index)] = i;
3335 ++index[0];
3336 lgFound = true;
3337 break;
3338 }
3339 }
3340 }
3341 while( lgFound );
3342
3343 if( index[0] == 0 )
3344 break;
3345
3346 grid->trackLen[index[1]] = index[0];
3347 ++index[1];
3348 }
3349
3350 grid->nTracks = index[1];
3351 return;
3352}
3353
3355 double val[], /* val[ndim] */
3356 long *nval,
3357 long *ndim)
3358{
3359 DEBUG_ENTRY( "CheckVal()" );
3360
3361 if( *ndim == 0 )
3362 *ndim = (long)grid->ndim;
3363 if( *ndim == 2 && *nval == 1 && grid->lgIsTeffLoggGrid )
3364 {
3365 /* default gravity is maximum gravity */
3366 val[*nval] = grid->val[1][grid->nval[1]-1];
3367 ++(*nval);
3368 }
3369 if( *ndim != (long)grid->ndim )
3370 {
3371 fprintf( ioQQQ, " A %ld-dim grid was requested, but a %ld-dim grid was found.\n",
3372 *ndim, (long)grid->ndim );
3374 }
3375 if( *nval < *ndim )
3376 {
3377 fprintf( ioQQQ, " A %ld-dim grid was requested, but only %ld parameters were entered.\n",
3378 *ndim, *nval );
3380 }
3381}
3382
3384 const double val[], /* val[ndim] */
3385 double *Tlow,
3386 double *Thigh)
3387{
3388 bool lgInvalid;
3389 long int i,
3390 *indlo,
3391 *indhi,
3392 *index,
3393 k,
3394 nd;
3395 double *aval;
3396
3397 DEBUG_ENTRY( "InterpolateRectGrid()" );
3398
3399 /* create some space */
3400 indlo = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) );
3401 indhi = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) );
3402 index = (long *)MALLOC((size_t)(grid->ndim*sizeof(long)) );
3403 aval = (double *)MALLOC((size_t)(grid->npar*sizeof(double)) );
3404
3405 ASSERT( rfield.lgContMalloc[rfield.nShape] );
3406 ASSERT( grid->nBlocksize == rfield.nupper*sizeof(realnum) );
3407
3408 /* save energy scale for check against code's in conorm (scale not yet defined when this routine called) */
3409 GetBins( grid, rfield.tNu[rfield.nShape] );
3410
3411# if DEBUGPRT
3412 /* check whether the models have the correct effective temperature, for debugging only */
3413 ValidateGrid( grid, 0.02 );
3414# endif
3415
3416 /* now generate pointers for models to use */
3417 for( nd=0; nd < grid->ndim; nd++ )
3418 {
3419 FindIndex( grid->val[nd], grid->nval[nd], val[nd], &indlo[nd], &indhi[nd], &lgInvalid );
3420 if( lgInvalid )
3421 {
3422 fprintf( ioQQQ,
3423 " Requested parameter %s = %.2f is not within the range %.2f to %.2f\n",
3424 grid->names[nd], val[nd], grid->val[nd][0], grid->val[nd][grid->nval[nd]-1] );
3426 }
3427 }
3428
3429 InterpolateModel( grid, val, aval, indlo, indhi, index, grid->ndim, rfield.tslop[rfield.nShape], IS_UNDEFINED );
3430
3431 /* print the parameters of the interpolated model */
3432 if( called.lgTalk )
3433 {
3434 if( grid->npar == 1 )
3435 fprintf( ioQQQ,
3436 " * c<< FINAL: %6s = %13.2f"
3437 " >>> *\n",
3438 grid->names[0], aval[0] );
3439 else if( grid->npar == 2 )
3440 fprintf( ioQQQ,
3441 " * c<< FINAL: %6s = %10.2f"
3442 " %6s = %8.5f >>> *\n",
3443 grid->names[0], aval[0], grid->names[1], aval[1] );
3444 else if( grid->npar == 3 )
3445 fprintf( ioQQQ,
3446 " * c<< FINAL: %6s = %7.0f"
3447 " %6s = %5.2f %6s = %5.2f >>> *\n",
3448 grid->names[0], aval[0], grid->names[1], aval[1],
3449 grid->names[2], aval[2] );
3450 else if( grid->npar >= 4 )
3451 {
3452 fprintf( ioQQQ,
3453 " * c<< FINAL: %4s = %7.0f"
3454 " %6s = %4.2f %6s = %5.2f %6s = ",
3455 grid->names[0], aval[0], grid->names[1], aval[1],
3456 grid->names[2], aval[2], grid->names[3] );
3457 fprintf( ioQQQ, PrintEfmt( "%9.2e", aval[3] ) );
3458 fprintf( ioQQQ, " >>> *\n" );
3459 }
3460 }
3461
3462 for( i=0; i < rfield.nupper; i++ )
3463 {
3464 rfield.tslop[rfield.nShape][i] = (realnum)pow((realnum)10.f,rfield.tslop[rfield.nShape][i]);
3465 if( rfield.tslop[rfield.nShape][i] < 1e-37 )
3466 rfield.tslop[rfield.nShape][i] = 0.;
3467 }
3468
3469 if( false )
3470 {
3471 FILE *ioBUG = fopen( "interpolated.txt", "w" );
3472 for( k=0; k < rfield.nupper; ++k )
3473 fprintf( ioBUG, "%e %e\n", rfield.tNu[rfield.nShape][k].Ryd(), rfield.tslop[rfield.nShape][k] );
3474 fclose( ioBUG );
3475 }
3476
3477 if( strcmp( grid->names[0], "Teff" ) == 0 )
3478 {
3479 if( ! lgValidModel( rfield.tNu[rfield.nShape], rfield.tslop[rfield.nShape], val[0], 0.10 ) )
3480 TotalInsanity();
3481 }
3482
3483 /* set limits for optimizer */
3484 SetLimits( grid, val[0], indlo, indhi, NULL, NULL, Tlow, Thigh );
3485
3486 FREE_CHECK( aval );
3487 FREE_CHECK( index );
3488 FREE_CHECK( indhi );
3489 FREE_CHECK( indlo );
3490 return;
3491}
3492
3494{
3495 long i;
3496
3497 DEBUG_ENTRY( "FreeGrid()" );
3498
3499 /* this was opened/allocated in InitGrid and subsidiaries,
3500 * this should become a destructor in C++ */
3501 fclose( grid->ioIN );
3502 FREE_CHECK( grid->telg );
3503 for( i = 0; i < grid->ndim; i++ )
3504 FREE_CHECK( grid->val[i] );
3505 FREE_CHECK( grid->val );
3506 FREE_CHECK( grid->nval );
3507 FREE_SAFE( grid->jlo );
3508 FREE_SAFE( grid->jhi );
3509 FREE_SAFE( grid->trackLen )
3510 FREE_SAFE( grid->jval );
3511 return;
3512}
3513
3515 const double val[],
3516 double aval[],
3517 const long indlo[],
3518 const long indhi[],
3519 long index[],
3520 long nd,
3521 vector<realnum>& flux1,
3522 IntStage stage)
3523{
3524 bool lgDryRun;
3525 long i, ind, j;
3526
3527 DEBUG_ENTRY( "InterpolateModel()" );
3528
3529 --nd;
3530
3531 lgDryRun = ( flux1.size() == 0 );
3532
3533 if( nd < 0 )
3534 {
3535 long n = JIndex(grid,index);
3536 if( stage == IS_FIRST )
3537 ind = ( grid->jlo[n] >= 0 ) ? grid->jlo[n] : grid->jhi[n];
3538 else if( stage == IS_SECOND )
3539 ind = ( grid->jhi[n] >= 0 ) ? grid->jhi[n] : grid->jlo[n];
3540 else if( grid->ndim == 1 )
3541 /* in this case grid->jlo[n] and grid->jhi[n] should be identical */
3542 ind = grid->jlo[n];
3543 else
3544 TotalInsanity();
3545
3546 if( ind < 0 )
3547 {
3548 fprintf( ioQQQ, " The requested interpolation could not be completed, sorry.\n" );
3549 fprintf( ioQQQ, " No suitable match was found for a model with" );
3550 for( i=0; i < grid->ndim; i++ )
3551 fprintf( ioQQQ, " %s=%.6g ", grid->names[i], grid->val[i][index[i]] );
3552 fprintf( ioQQQ, "\n" );
3554 }
3555
3556 for( i=0; i < grid->npar; i++ )
3557 aval[i] = grid->telg[ind].par[i];
3558
3559 if( !lgDryRun )
3560 {
3561 for( i=0; i < grid->ndim && called.lgTalk; i++ )
3562 {
3563 if( !fp_equal(grid->val[i][index[i]],aval[i],10) )
3564 {
3565 fprintf( ioQQQ, " No exact match was found for a model with" );
3566 for( j=0; j < grid->ndim; j++ )
3567 fprintf( ioQQQ, " %s=%.6g ", grid->names[j], grid->val[j][index[j]] );
3568 fprintf( ioQQQ, "- using the following model instead:\n" );
3569 break;
3570 }
3571 }
3572
3573 GetModel( grid, ind, flux1, lgVERBOSE, lgTAKELOG );
3574 }
3575 }
3576 else
3577 {
3578 vector<realnum> flux2(rfield.nupper);
3579 double *aval2;
3580
3581# if !defined NDEBUG
3582 const realnum SECURE = 10.f*FLT_EPSILON;
3583# endif
3584
3585 aval2 = (double*)MALLOC((size_t)(grid->npar*sizeof(double)) );
3586
3587 /* Interpolation is carried out first in the parameter with nd == 0 (usually
3588 * Teff), then the parameter with nd == 1 (usually log(g)), etc. One or two
3589 * atmosphere models are read depending on whether the parameter was matched
3590 * exactly or not. If needed, logarithmic interpolation is done.
3591 */
3592
3593 if( nd == 1 )
3594 stage = IS_FIRST;
3595
3596 index[nd] = indlo[nd];
3597 InterpolateModel( grid, val, aval, indlo, indhi, index, nd, flux1, stage );
3598
3599 if( nd == 1 )
3600 stage = IS_SECOND;
3601
3602 index[nd] = indhi[nd];
3603 vector<realnum> empty;
3604 InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, empty, stage );
3605
3606 if( !fp_equal(aval2[nd],aval[nd],10) )
3607 {
3608 double fr1, fr2, fc1 = 0., fc2 = 0.;
3609
3610 if( !lgDryRun )
3611 InterpolateModel( grid, val, aval2, indlo, indhi, index, nd, flux2, stage );
3612
3613 fr1 = (aval2[nd]-val[nd])/(aval2[nd]-aval[nd]);
3614 /* when interpolating in log(g) it can happen that fr1 is outside the range 0 .. 1
3615 * this can be the case when the requested log(g) was not present in the grid
3616 * and it had to be approximated by another model. In this case do not extrapolate */
3617 if( nd == 1 )
3618 fr1 = MIN2( MAX2( fr1, 0. ), 1. );
3619 fr2 = 1. - fr1;
3620
3621 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3622
3623 if( !lgDryRun )
3624 {
3625# if DEBUGPRT
3626 fprintf( ioQQQ, "interpolation nd=%ld fr1=%g\n", nd, fr1 );
3627# endif
3628
3629 /* special treatment for high-temperature Rauch models */
3630 if( nd == 0 && strcmp( grid->names[nd], "Teff" ) == 0 )
3631 {
3632 /* The following is an approximate scaling to use for the range of
3633 * temperatures above 200000 K in the H-Ca Rauch models where the
3634 * temperature steps are large and thus the interpolations are over
3635 * large ranges. For the lower temperatures I assume that there is
3636 * no need for this.
3637 *
3638 * It should be remembered that this interpolation is not exact, and
3639 * the possible error at high temperatures might be large enough to
3640 * matter. (Kevin Volk)
3641 */
3642 fc1 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->val[nd][indlo[nd]])*4. : 0.;
3643 fc2 = ( val[nd] > 200000. ) ? log10(val[nd]/grid->val[nd][indhi[nd]])*4. : 0.;
3644 }
3645
3646 for( i=0; i < rfield.nupper; ++i )
3647 flux1[i] = (realnum)(fr1*(flux1[i]+fc1) + fr2*(flux2[i]+fc2));
3648 }
3649
3650 for( i=0; i < grid->npar; i++ )
3651 aval[i] = fr1*aval[i] + fr2*aval2[i];
3652 }
3653
3654 FREE_CHECK( aval2 );
3655 }
3656 return;
3657}
3658
3660 const double val[],
3661 double aval[],
3662 const long indlo[],
3663 const long indhi[],
3664 long index[],
3665 long nd,
3666 long off,
3667 vector<realnum>& flux1)
3668{
3669 long i, ind;
3670
3671 DEBUG_ENTRY( "InterpolateModelCoStar()" );
3672
3673 if( nd == 2 )
3674 {
3675 ind = ( index[1] == 0 ) ? indlo[index[0]] : indhi[index[0]];
3676
3677 GetModel( grid, ind, flux1, lgVERBOSE, lgTAKELOG );
3678
3679 for( i=0; i < grid->npar; i++ )
3680 aval[i] = grid->telg[ind].par[i];
3681 }
3682 else
3683 {
3684 bool lgSkip;
3685# if !defined NDEBUG
3686 const realnum SECURE = 10.f*FLT_EPSILON;
3687# endif
3688
3689 /* Interpolation is carried out first along evolutionary tracks, then
3690 * in between evolutionary tracks. Between 1 and 4 atmosphere models are read
3691 * depending on whether the parameter/track was matched exactly or not.
3692 */
3693
3694 index[nd] = 0;
3695 InterpolateModelCoStar( grid, val, aval, indlo, indhi, index, nd+1, off, flux1 );
3696
3697 lgSkip = ( nd == 1 ) ? ( indhi[index[0]] == indlo[index[0]] ) :
3698 ( indlo[0] == indlo[1] && indhi[0] == indhi[1] );
3699
3700 if( ! lgSkip )
3701 {
3702 vector<realnum> flux2(rfield.nupper);
3703 double fr1, fr2, *aval2;
3704
3705 aval2 = (double*)MALLOC((size_t)(grid->npar*sizeof(double)) );
3706
3707 index[nd] = 1;
3708 InterpolateModelCoStar( grid, val, aval2, indlo, indhi, index, nd+1, off, flux2 );
3709
3710 fr1 = (aval2[nd+off]-val[nd])/(aval2[nd+off]-aval[nd+off]);
3711 fr2 = 1. - fr1;
3712
3713# if DEBUGPRT
3714 fprintf( ioQQQ, "interpolation nd=%ld fr1=%g\n", nd, fr1 );
3715# endif
3716
3717 ASSERT( 0.-SECURE <= fr1 && fr1 <= 1.+SECURE );
3718
3719 for( i=0; i < rfield.nupper; ++i )
3720 flux1[i] = (realnum)(fr1*flux1[i] + fr2*flux2[i]);
3721
3722 for( i=0; i < grid->npar; i++ )
3723 aval[i] = fr1*aval[i] + fr2*aval2[i];
3724
3725 FREE_CHECK( aval2 );
3726 }
3727 }
3728 return;
3729}
3730
3732 vector<Energy>& ener)
3733{
3734 DEBUG_ENTRY( "GetBins()" );
3735
3736 /* make sure ident is exactly 12 characters long, otherwise output won't fit */
3737 ASSERT( strlen(grid->ident) == 12 );
3738
3739 ASSERT( grid->nBlocksize == rfield.nupper*sizeof(realnum) );
3740
3741 /* skip over ind stars */
3742 /* >>chng 01 oct 18, add nOffset */
3743 if( fseek( grid->ioIN, (long)(grid->nOffset), SEEK_SET ) != 0 )
3744 {
3745 fprintf( ioQQQ, " Error finding atmosphere frequency bins\n");
3747 }
3748
3749 vector<realnum> data(rfield.nupper);
3750 if( fread( get_ptr(data), 1, grid->nBlocksize, grid->ioIN ) != grid->nBlocksize )
3751 {
3752 fprintf( ioQQQ, " Error reading atmosphere frequency bins\n" );
3754 }
3755
3756 for( long i=0; i < rfield.nupper; ++i )
3757 ener[i].set(data[i]);
3758
3759 return;
3760}
3761
3763 long ind,
3764 vector<realnum>& flux,
3765 bool lgTalk,
3766 bool lgTakeLog)
3767{
3768 long i;
3769
3770 DEBUG_ENTRY( "GetModel()" );
3771
3772 /* add 1 to account for frequency grid that is stored in front of all the atmospheres */
3773 ind++;
3774
3775 /* make sure ident is exactly 12 characters long, otherwise output won't fit */
3776 ASSERT( strlen(grid->ident) == 12 );
3777 /* ind == 0 is the frequency grid, ind == 1 .. nmods are the atmosphere models */
3778 ASSERT( ind >= 0 && ind <= grid->nmods );
3779
3780 /* skip over ind stars */
3781 /* >>chng 01 oct 18, add nOffset */
3782 if( fseek( grid->ioIN, (long)(ind*grid->nBlocksize+grid->nOffset), SEEK_SET ) != 0 )
3783 {
3784 fprintf( ioQQQ, " Error seeking atmosphere %ld\n", ind );
3786 }
3787
3788 if( fread( get_ptr(flux), 1, grid->nBlocksize, grid->ioIN ) != grid->nBlocksize )
3789 {
3790 fprintf( ioQQQ, " Error trying to read atmosphere %ld\n", ind );
3792 }
3793
3794 /* print the parameters of the atmosphere model */
3795 if( called.lgTalk && lgTalk )
3796 {
3797 /* ind-1 below since telg doesn't have the entry for the frequency grid */
3798 if( grid->npar == 1 )
3799 fprintf( ioQQQ,
3800 " * c<< %s model%5ld read. "
3801 " %6s = %13.2f >>> *\n",
3802 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0] );
3803 else if( grid->npar == 2 )
3804 fprintf( ioQQQ,
3805 " * c<< %s model%5ld read. "
3806 " %6s = %10.2f %6s = %8.5f >>> *\n",
3807 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0],
3808 grid->names[1], grid->telg[ind-1].par[1] );
3809 else if( grid->npar == 3 )
3810 fprintf( ioQQQ,
3811 " * c<< %s model%5ld read. "
3812 " %6s=%7.0f %6s=%5.2f %6s=%5.2f >>> *\n",
3813 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0],
3814 grid->names[1], grid->telg[ind-1].par[1],
3815 grid->names[2], grid->telg[ind-1].par[2] );
3816 else if( grid->npar >= 4 )
3817 {
3818 fprintf( ioQQQ,
3819 " * c< %s mdl%4ld"
3820 " %4s=%5.0f %6s=%4.2f %6s=%5.2f %6s=",
3821 grid->ident, ind, grid->names[0], grid->telg[ind-1].par[0],
3822 grid->names[1], grid->telg[ind-1].par[1],
3823 grid->names[2], grid->telg[ind-1].par[2], grid->names[3] );
3824 fprintf( ioQQQ, PrintEfmt( "%9.2e", grid->telg[ind-1].par[3] ) );
3825 fprintf( ioQQQ, " >> *\n" );
3826 }
3827 }
3828
3829 if( lgTakeLog )
3830 {
3831 /* convert to logs since we will interpolate in log flux */
3832 for( i=0; i < rfield.nupper; ++i )
3833 {
3834 // the keyword volatile is needed to work around a
3835 // compiler bug in g++ versions 4.7.0 and later
3836 // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65425
3837 volatile double help = flux[i];
3838 if( help > 0. )
3839 help = log10(help);
3840 else
3841 help = -99999.;
3842 flux[i] = realnum(help);
3843 }
3844 }
3845 return;
3846}
3847
3849 double val,
3850 const long indlo[],
3851 const long indhi[],
3852 const long useTr[],
3853 const realnum ValTr[],
3854 double *loLim,
3855 double *hiLim)
3856{
3857 DEBUG_ENTRY( "SetLimits()" );
3858
3859 if( optimize.lgVarOn )
3860 {
3861 int ptr0,ptr1;
3862 long index[MDIM], j;
3863 const double SECURE = (1. + 20.*(double)FLT_EPSILON);
3864
3865 *loLim = +DBL_MAX;
3866 *hiLim = -DBL_MAX;
3867
3868 switch( grid->imode )
3869 {
3870 case IM_RECT_GRID:
3871 *loLim = -DBL_MAX;
3872 *hiLim = +DBL_MAX;
3873 SetLimitsSub( grid, val, indlo, indhi, index, grid->ndim, loLim, hiLim );
3874 break;
3878 for( j=0; j < grid->nTracks; j++ )
3879 {
3880 if( ValTr[j] != -FLT_MAX )
3881 {
3882 /* M_ZAMS is already logarithm, Teff is linear */
3883 double temp = ( grid->imode == IM_COSTAR_MZAMS_AGE ) ?
3884 pow(10.,(double)ValTr[j]) : ValTr[j];
3885 *loLim = MIN2(*loLim,temp);
3886 *hiLim = MAX2(*hiLim,temp);
3887 }
3888 }
3889 break;
3891 index[0] = 0;
3892 index[1] = useTr[0];
3893 ptr0 = grid->jval[JIndex(grid,index)];
3894 index[1] = useTr[1];
3895 ptr1 = grid->jval[JIndex(grid,index)];
3896 *loLim = MAX2(grid->telg[ptr0].par[3],grid->telg[ptr1].par[3]);
3897# if DEBUGPRT
3898 printf( "set limit 0: (models %d, %d) %f %f\n",
3899 ptr0+1, ptr1+1, grid->telg[ptr0].par[3], grid->telg[ptr1].par[3] );
3900# endif
3901 index[0] = grid->trackLen[useTr[0]]-1;
3902 index[1] = useTr[0];
3903 ptr0 = grid->jval[JIndex(grid,index)];
3904 index[0] = grid->trackLen[useTr[1]]-1;
3905 index[1] = useTr[1];
3906 ptr1 = grid->jval[JIndex(grid,index)];
3907 *hiLim = MIN2(grid->telg[ptr0].par[3],grid->telg[ptr1].par[3]);
3908# if DEBUGPRT
3909 printf( "set limit 1: (models %d, %d) %f %f\n",
3910 ptr0+1, ptr1+1, grid->telg[ptr0].par[3], grid->telg[ptr1].par[3] );
3911# endif
3912 break;
3913 default:
3914 fprintf( ioQQQ, " SetLimits called with insane value for imode: %d.\n", grid->imode );
3916 }
3917
3918 ASSERT( fabs(*loLim) < DBL_MAX && fabs(*hiLim) < DBL_MAX );
3919
3920 /* check sanity of optimization limits */
3921 if( *hiLim <= *loLim )
3922 {
3923 fprintf( ioQQQ, " no room to optimize: lower limit %.4f, upper limit %.4f.\n",
3924 *loLim,*hiLim );
3926 }
3927
3928 /* make a bit of room for round-off errors */
3929 *loLim *= SECURE;
3930 *hiLim /= SECURE;
3931
3932# if DEBUGPRT
3933 printf("set limits: %g %g\n",*loLim,*hiLim);
3934# endif
3935 }
3936 else
3937 {
3938 *loLim = 0.;
3939 *hiLim = 0.;
3940 }
3941 return;
3942}
3943
3945 double val,
3946 const long indlo[],
3947 const long indhi[],
3948 long index[],
3949 long nd,
3950 double *loLim,
3951 double *hiLim)
3952{
3953 long n;
3954
3955 DEBUG_ENTRY( "SetLimitsSub()" );
3956
3957 --nd;
3958
3959 if( nd < 1 )
3960 {
3961 double loLoc = +DBL_MAX;
3962 double hiLoc = -DBL_MAX;
3963
3964 for( index[0]=0; index[0] < grid->nval[0]; ++index[0] )
3965 {
3966 /* grid->val[0][i] is the array of Par0 values (Teff/Age/...) in the
3967 * grid, which it is sorted in strict monotonically increasing order.
3968 * This routine searches for the largest range [loLoc,hiLoc] in Par0
3969 * such that loLoc <= val <= hiLoc, and at least one model exists for
3970 * each Par0 value in this range. This assures that interpolation is
3971 * safe and the optimizer will not trip... */
3972 n = JIndex(grid,index);
3973 if( grid->jlo[n] < 0 && grid->jhi[n] < 0 )
3974 {
3975 /* there are no models with this value of Par0 */
3976 /* this value of Par0 should be outside of allowed range */
3977 if( grid->val[0][index[0]] < val )
3978 loLoc = DBL_MAX;
3979 /* this is beyond the legal range, so terminate the search */
3980 if( grid->val[0][index[0]] > val )
3981 break;
3982 }
3983 else
3984 {
3985 /* there are models with this value of Par0 */
3986 /* update range to include this value of Par0 */
3987 if( grid->val[0][index[0]] <= val )
3988 {
3989 /* remember lowest legal value of loLoc
3990 * -> only update if previous value was illegal */
3991 if( loLoc == DBL_MAX )
3992 loLoc = grid->val[0][index[0]];
3993 }
3994 if( grid->val[0][index[0]] >= val )
3995 {
3996 /* remember highest legal value of hiLoc
3997 * -> always update */
3998 hiLoc = grid->val[0][index[0]];
3999 }
4000 }
4001 }
4002
4003 ASSERT( fabs(loLoc) < DBL_MAX && fabs(hiLoc) < DBL_MAX && loLoc <= hiLoc );
4004
4005 *loLim = MAX2(*loLim,loLoc);
4006 *hiLim = MIN2(*hiLim,hiLoc);
4007 }
4008 else
4009 {
4010 index[nd] = indlo[nd];
4011 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
4012
4013 if( indhi[nd] != indlo[nd] )
4014 {
4015 index[nd] = indhi[nd];
4016 SetLimitsSub( grid, val, indlo, indhi, index, nd, loLim, hiLim );
4017 }
4018 }
4019 return;
4020}
4021
4023 bool lgList)
4024{
4025 long i, *index, j, jsize, nd;
4026 double *val;
4027
4028 DEBUG_ENTRY( "InitIndexArrays()" );
4029
4030 ASSERT( grid->telg != NULL );
4031 ASSERT( grid->nmods > 0 );
4032
4033 jsize = 1;
4034
4035 /* this loop creates a list of all unique model parameter values in increasing order */
4036 for( nd=0; nd < grid->ndim; nd++ )
4037 {
4038 double pval = grid->telg[0].par[nd];
4039 grid->val[nd][0] = pval;
4040 grid->nval[nd] = 1;
4041
4042 for( i=1; i < grid->nmods; i++ )
4043 {
4044 bool lgOutOfRange;
4045 long i1, i2;
4046
4047 pval = grid->telg[i].par[nd];
4048 FindIndex( grid->val[nd], grid->nval[nd], pval, &i1, &i2, &lgOutOfRange );
4049 /* if i1 < i2, the new parameter value was not present yet and
4050 * it needs to be inserted in between i1 and i2 --> first move
4051 * all entries from i2 to grid->nval[nd]-1 one slot upward and
4052 * then insert the new value at i2; this also works correctly
4053 * if lgOutOfRange is set, hence no special check is needed */
4054 if( i1 < i2 )
4055 {
4056 /* val[nd] has grid->nmods entries, so cannot overflow */
4057 for( j = grid->nval[nd]-1; j >= i2; j-- )
4058 grid->val[nd][j+1] = grid->val[nd][j];
4059 grid->val[nd][i2] = pval;
4060 grid->nval[nd]++;
4061 }
4062 }
4063
4064 jsize *= grid->nval[nd];
4065
4066# if DEBUGPRT
4067 printf( "%s[%ld]:", grid->names[nd], grid->nval[nd] );
4068 for( i=0; i < grid->nval[nd]; i++ )
4069 printf( " %g", grid->val[nd][i] );
4070 printf( "\n" );
4071# endif
4072 }
4073
4074 index = (long *)MALLOC( sizeof(long)*grid->ndim );
4075 val = (double *)MALLOC( sizeof(double)*grid->ndim );
4076
4077 /* this memory will be freed in the calling function */
4078 grid->jlo = (long *)MALLOC( sizeof(long)*jsize );
4079 grid->jhi = (long *)MALLOC( sizeof(long)*jsize );
4080
4081 /* set up square array of model indices; this will be used to
4082 * choose the correct models for the interpolation process */
4083 FillJ( grid, index, val, grid->ndim, lgList );
4084
4085 FREE_CHECK( val );
4086 FREE_CHECK( index );
4087
4088 if( lgList )
4090 return;
4091}
4092
4094 long index[], /* index[grid->ndim] */
4095 double val[], /* val[grid->ndim] */
4096 long nd,
4097 bool lgList)
4098{
4099 DEBUG_ENTRY( "FillJ()" );
4100
4101 --nd;
4102
4103 if( nd < 0 )
4104 {
4105 long n = JIndex(grid,index);
4106 SearchModel( grid->telg, grid->lgIsTeffLoggGrid, grid->nmods, val, grid->ndim,
4107 &grid->jlo[n], &grid->jhi[n] );
4108 }
4109 else
4110 {
4111 for( index[nd]=0; index[nd] < grid->nval[nd]; index[nd]++ )
4112 {
4113 val[nd] = grid->val[nd][index[nd]];
4114 FillJ( grid, index, val, nd, lgList );
4115 }
4116 }
4117
4118 if( lgList && nd == MIN2(grid->ndim-1,1) )
4119 {
4120 fprintf( ioQQQ, "\n" );
4121 if( grid->ndim > 2 )
4122 {
4123 fprintf( ioQQQ, "subgrid for" );
4124 for( long n = nd+1; n < grid->ndim; n++ )
4125 fprintf( ioQQQ, " %s=%g", grid->names[n], val[n] );
4126 fprintf( ioQQQ, ":\n\n" );
4127 }
4128 if( grid->ndim > 1 )
4129 {
4130 fprintf( ioQQQ, "%6.6s\\%6.6s |", grid->names[0], grid->names[1] );
4131 for( long n = 0; n < grid->nval[1]; n++ )
4132 fprintf( ioQQQ, " %9.3g", grid->val[1][n] );
4133 fprintf( ioQQQ, "\n" );
4134 fprintf( ioQQQ, "--------------|" );
4135 for( long n = 0; n < grid->nval[1]; n++ )
4136 fprintf( ioQQQ, "----------" );
4137 }
4138 else
4139 {
4140 fprintf( ioQQQ, "%13.13s |\n", grid->names[0] );
4141 fprintf( ioQQQ, "--------------|----------" );
4142 }
4143 fprintf( ioQQQ, "\n" );
4144 for( index[0]=0; index[0] < grid->nval[0]; index[0]++ )
4145 {
4146 fprintf( ioQQQ, "%13.7g |", grid->val[0][index[0]] );
4147 if( grid->ndim > 1 )
4148 {
4149 for( index[1]=0; index[1] < grid->nval[1]; index[1]++ )
4150 if( grid->jlo[JIndex(grid,index)] == grid->jhi[JIndex(grid,index)] &&
4151 grid->jlo[JIndex(grid,index)] >= 0 )
4152 fprintf( ioQQQ, " %9ld", grid->jlo[JIndex(grid,index)]+1 );
4153 else
4154 fprintf( ioQQQ, " --" );
4155 }
4156 else
4157 {
4158 fprintf( ioQQQ, " %9ld", grid->jlo[JIndex(grid,index)]+1 );
4159 }
4160 fprintf( ioQQQ, "\n" );
4161 }
4162 fprintf( ioQQQ, "\n" );
4163 }
4164 return;
4165}
4166
4168 const long index[]) /* index[grid->ndim] */
4169{
4170 long i, ind, mul;
4171
4172 DEBUG_ENTRY( "JIndex()" );
4173
4174 ind = 0;
4175 mul = 1;
4176 for( i=0; i < grid->ndim; i++ )
4177 {
4178 ind += index[i]*mul;
4179 mul *= grid->nval[i];
4180 }
4181 return ind;
4182}
4183
4184STATIC void SearchModel(const mpp telg[], /* telg[nmods] */
4185 bool lgIsTeffLoggGrid,
4186 long nmods,
4187 const double val[], /* val[ndim] */
4188 long ndim,
4189 long *index_low,
4190 long *index_high)
4191{
4192 long i, nd;
4193 double alogg_low = -DBL_MAX, alogg_high = DBL_MAX;
4194
4195 DEBUG_ENTRY( "SearchModel()" );
4196
4197 /* given values for the model parameters, this routine searches for the atmosphere
4198 * model that is the best match. If all parameters can be matched simultaneously the
4199 * choice is obvious, but this cannot always be achieved (typically for high Teff, the
4200 * low log(g) models will be missing). If lgIsTeffLoggGrid is true, the rule is that
4201 * all parameters except log(g) must always be matched (such a model is not always
4202 * guaranteed to exist). If all requested parameters can be matched exactly, both
4203 * index_low and index_high will point to that model. If all parameters except log(g)
4204 * can be matched exactly, it will return the model with the lowest log(g) value larger
4205 * than the requested value in index_high, and the model with the highest log(g) value
4206 * lower than the requested value in index_low. If either requirement cannot be
4207 * fulfilled, -2 will be returned. When lgIsTeffLoggGrid is false, all parameters must
4208 * be matched and both index_low and index_high will point to that model. If no such
4209 * model can be found, -2 will be returned. */
4210
4211 *index_low = *index_high = -2;
4212 for( i=0; i < nmods; i++ )
4213 {
4214 bool lgNext = false;
4215 /* ignore models with different parameters */
4216 for( nd=0; nd < ndim; nd++ )
4217 {
4218 if( nd != 1 && !fp_equal(telg[i].par[nd],val[nd],10) )
4219 {
4220 lgNext = true;
4221 break;
4222 }
4223 }
4224 if( lgNext )
4225 continue;
4226
4227 /* an exact match is found */
4228 if( ndim == 1 || fp_equal(telg[i].par[1],val[1],10) )
4229 {
4230 *index_low = i;
4231 *index_high = i;
4232 return;
4233 }
4234 if( lgIsTeffLoggGrid )
4235 {
4236 /* keep a record of the highest log(g) model smaller than alogg */
4237 if( telg[i].par[1] < val[1] && telg[i].par[1] > alogg_low )
4238 {
4239 *index_low = i;
4240 alogg_low = telg[i].par[1];
4241 }
4242 /* also keep a record of the lowest log(g) model greater than alogg */
4243 if( telg[i].par[1] > val[1] && telg[i].par[1] < alogg_high )
4244 {
4245 *index_high = i;
4246 alogg_high = telg[i].par[1];
4247 }
4248 }
4249 }
4250 return;
4251}
4252
4253STATIC void FindIndex(const double xval[], /* xval[NVAL] */
4254 long NVAL,
4255 double x,
4256 long *ind1,
4257 long *ind2,
4258 bool *lgInvalid)
4259{
4260 bool lgOutLo, lgOutHi;
4261 long i;
4262
4263 DEBUG_ENTRY( "FindIndex()" );
4264
4265 /* this routine searches for indices ind1, ind2 such that
4266 * xval[ind1] < x < xval[ind2]
4267 * if x is equal to one of the values in xval, then
4268 * ind1 == ind2 and xval[ind1] == x
4269 *
4270 * if x is outside the range xval[0] ... xval[NVAL-1]
4271 * then lgInvalid will be set to true
4272 *
4273 * NB NB -- this routine implicitly assumes that xval is
4274 * strictly monotonically increasing!
4275 */
4276
4277 ASSERT( NVAL > 0 );
4278
4279 /* is x outside of range xval[0] ... xval[NVAL-1]? */
4280 lgOutLo = ( x-xval[0] < -10.*DBL_EPSILON*fabs(xval[0]) );
4281 lgOutHi = ( x-xval[NVAL-1] > 10.*DBL_EPSILON*fabs(xval[NVAL-1]) );
4282
4283 if( lgOutLo || lgOutHi )
4284 {
4285 /* pretend there are two fictitious array elements
4286 * xval[-1] = -Inf and xval[NVAL] = +Inf,
4287 * and return ind1 and ind2 accordingly. This behavior
4288 * is needed for InitIndexArrays() to work correctly */
4289 *ind1 = lgOutLo ? -1 : NVAL-1;
4290 *ind2 = lgOutLo ? 0 : NVAL;
4291 *lgInvalid = true;
4292 return;
4293 }
4294
4295 *lgInvalid = false;
4296
4297 /* there are more efficient ways of doing this, e.g. a binary search.
4298 * However, the xval arrays typically only have 1 or 2 dozen elements,
4299 * so the overhead is negligible and the clarity of this code is preferred */
4300
4301 /* first look for an "exact" match */
4302 for( i=0; i < NVAL; i++ )
4303 {
4304 if( fp_equal(xval[i],x,10) )
4305 {
4306 *ind1 = i;
4307 *ind2 = i;
4308 return;
4309 }
4310 }
4311
4312 /* no match was found -> bracket the x value */
4313 for( i=0; i < NVAL-1; i++ )
4314 {
4315 if( xval[i] < x && x < xval[i+1] )
4316 {
4317 *ind1 = i;
4318 *ind2 = i+1;
4319 return;
4320 }
4321 }
4322
4323 /* this should never be reached ! */
4324 fprintf( ioQQQ, " insanity in FindIndex\n" );
4325 ShowMe();
4327}
4328
4329STATIC bool lgFileReadable(const char *chFnam, process_counter& pc, access_scheme scheme)
4330{
4331 DEBUG_ENTRY( "lgFileReadable()" );
4332
4333 FILE *ioIN;
4334
4335 ioIN = open_data( chFnam, "r", scheme );
4336 if( ioIN != NULL )
4337 {
4338 fclose( ioIN );
4339 ++pc.nFound;
4340 return true;
4341 }
4342 else
4343 {
4344 return false;
4345 }
4346}
4347
4348/*ValidateGrid: check each model in the grid to see if it has the correct Teff */
4350 double toler)
4351{
4352 long i, k, nd;
4353 vector<Energy> anu(rfield.nupper);
4354 vector<realnum> flux(rfield.nupper);
4355
4356 DEBUG_ENTRY( "ValidateGrid()" );
4357
4358 if( strcmp( grid->names[0], "Teff" ) != 0 )
4359 {
4360 return;
4361 }
4362
4363 GetBins( grid, anu );
4364
4365 for( i=0; i < grid->nmods; i++ )
4366 {
4367 fprintf( ioQQQ, "testing model %ld ", i+1 );
4368 for( nd=0; nd < grid->npar; nd++ )
4369 fprintf( ioQQQ, " %s %g", grid->names[nd], grid->telg[i].par[nd] );
4370
4371 GetModel( grid, i, flux, lgSILENT, lgLINEAR );
4372
4373 if( lgValidModel( anu, flux, grid->telg[i].par[0], toler ) )
4374 fprintf( ioQQQ, " OK\n" );
4375
4376 if( false )
4377 {
4378 FILE *ioBUG = fopen( "atmosphere_dump.txt", ( i == 0 ) ? "w" : "a" );
4379
4380 fprintf( ioBUG, "######## MODEL %ld", i+1 );
4381 for( nd=0; nd < grid->npar; nd++ )
4382 fprintf( ioBUG, " %s %g", grid->names[nd], grid->telg[i].par[nd] );
4383 fprintf( ioBUG, "####################\n" );
4384
4385 for( k=0; k < rfield.nupper; ++k )
4386 fprintf( ioBUG, "%e %e\n", anu[k].Ryd(), flux[k] );
4387
4388 fclose( ioBUG );
4389 }
4390 }
4391 return;
4392}
4393
4394STATIC bool lgValidModel(const vector<Energy>& anu,
4395 const vector<realnum>& flux,
4396 double Teff,
4397 double toler)
4398{
4399 bool lgPassed = true;
4400 long k;
4401 double chk, lumi;
4402
4403 DEBUG_ENTRY( "lgValidModel()" );
4404
4405 ASSERT( Teff > 0. );
4406
4407 lumi = 0.;
4408 /* rebinned models are in cgs F_nu units */
4409 for( k=1; k < rfield.nupper; k++ )
4410 lumi += (anu[k].Ryd() - anu[k-1].Ryd())*(flux[k] + flux[k-1])/2.;
4411
4412 /* now convert luminosity to effective temperature */
4413 chk = pow(lumi*FR1RYD/STEFAN_BOLTZ,0.25);
4414 /* the allowed tolerance is set by the caller in toler */
4415 if( fabs(Teff - chk) > toler*Teff ) {
4416 fprintf( ioQQQ, "\n*** WARNING, Teff discrepancy for this model, expected Teff %.2f, ", Teff);
4417 fprintf( ioQQQ, "integration yielded Teff %.2f, delta %.2f%%\n", chk, (chk/Teff-1.)*100. );
4418 lgPassed = false;
4419 }
4420 return lgPassed;
4421}
4422
4423/*RebinAtmosphere: generic routine for rebinning atmospheres onto Cloudy grid */
4424STATIC void RebinAtmosphere(long nCont, /* the number of points in the incident continuum*/
4425 const realnum StarEner[], /* StarEner[nCont], the freq grid for the model, in Ryd*/
4426 const realnum StarFlux[], /* StarFlux[nCont], the original model flux */
4427 realnum CloudyFlux[], /* CloudyFlux[NC_ELL], the model flux on the cloudy grid */
4428 long nEdge, /* the number of bound-free continuum edges in AbsorbEdge */
4429 const realnum AbsorbEdge[]) /* AbsorbEdge[nEdge], energies of the edges */
4430{
4431 bool lgDone;
4432 long int ind,
4433 j,
4434 k;
4435 /* >>chng 00 jun 02, demoted next two to realnum, PvH */
4436 realnum BinHigh,
4437 BinLow,
4438 BinMid,
4439 BinNext,
4440 *EdgeLow=NULL,
4441 *EdgeHigh=NULL,
4442 *StarPower;
4443
4444 DEBUG_ENTRY( "RebinAtmosphere()" );
4445
4446 if( nEdge > 0 )
4447 {
4448 EdgeLow = (realnum*)MALLOC( sizeof(realnum)*(unsigned)nEdge );
4449 EdgeHigh = (realnum*)MALLOC( sizeof(realnum)*(unsigned)nEdge );
4450 }
4451
4452 /* this loop should be before the next loop, otherwise models with a
4453 * very strong He II edge (i.e. no flux beyond that edge) will fail */
4454 for( j=0; j < nEdge; j++ )
4455 {
4456 ind = RebinFind(StarEner,nCont,AbsorbEdge[j]);
4457
4458 /* sanity check */
4459 ASSERT( ind >= 0 && ind+1 < nCont );
4460
4461 EdgeLow[j] = StarEner[ind];
4462 EdgeHigh[j] = StarEner[ind+1];
4463 }
4464
4465 /* cut off that part of the Wien tail that evaluated to zero */
4466 /* >> chng 05 nov 22, inverted loop, slightly faster PvH */
4467 /*for( j=nCont-1; j >= 0; j-- )*/
4468 for( j=0; j < nCont; j++ )
4469 {
4470 if( StarFlux[j] == 0.f )
4471 {
4472 nCont = j;
4473 break;
4474 }
4475 }
4476 ASSERT( nCont > 0 );
4477
4478 StarPower = (realnum *)MALLOC( sizeof(realnum)*(unsigned)(nCont-1) );
4479
4480 for( j=0; j < nCont-1; j++ )
4481 {
4482 double ratio_x, ratio_y;
4483
4484 /* >>chng 05 nov 22, add sanity check to prevent invalid fp operations */
4485 ASSERT( StarEner[j+1] > StarEner[j] );
4486
4487 /* >>chng 06 aug 11, on some systems (e.g., macbook pro) y/x can get evaluated as y*(1/x);
4488 * this causes overflows if x is a denormalized number, hence we force a cast to double, PvH */
4489 ratio_x = (double)StarEner[j+1]/(double)StarEner[j];
4490 ratio_y = (double)StarFlux[j+1]/(double)StarFlux[j];
4491 StarPower[j] = (realnum)(log(ratio_y)/log(ratio_x));
4492 }
4493
4494 for( j=0; j < rfield.nupper; j++ )
4495 {
4496 /* >>chng 05 nov 22, modified BinLow, BinHigh, BinNext to make boundaries match exactly, PvH */
4497 /* BinLow is lower bound of this continuum cell */
4498 BinLow = ( j > 0 ) ?
4499 (realnum)sqrt(rfield.anu[j-1]*rfield.anu[j]) : (realnum)sqrt(POW3(rfield.anu[0])/rfield.anu[1]);
4500
4501 /* BinHigh is upper bound of this continuum cell */
4502 BinHigh = ( j+1 < rfield.nupper ) ?
4503 (realnum)sqrt(rfield.anu[j]*rfield.anu[j+1]) : rfield.anu[rfield.nupper-1];
4504
4505 /* BinNext is upper bound of next continuum cell */
4506 BinNext = ( j+2 < rfield.nupper ) ?
4507 (realnum)sqrt(rfield.anu[j+1]*rfield.anu[j+2]) : rfield.anu[rfield.nupper-1];
4508
4509 lgDone = false;
4510
4511 /* >>chng 00 aug 14, take special care not to interpolate over major edges,
4512 * the region in between EdgeLow and EdgeHigh should be avoided,
4513 * the spectrum is extremely steep there, leading to significant roundoff error, PvH */
4514 for( k=0; k < nEdge; k++ )
4515 {
4516 if( BinLow < EdgeLow[k] && BinNext > EdgeHigh[k] )
4517 {
4518 BinMid = 0.99999f*EdgeLow[k];
4519 CloudyFlux[j] = RebinSingleCell(BinLow,BinMid,StarEner,StarFlux,StarPower,nCont);
4520 j++;
4521
4522 /* sanity check */
4523 ASSERT( j < rfield.nupper );
4524
4525 BinMid = 1.00001f*EdgeHigh[k];
4526 CloudyFlux[j] = RebinSingleCell(BinMid,BinNext,StarEner,StarFlux,StarPower,nCont);
4527 lgDone = true;
4528 break;
4529 }
4530 }
4531
4532 /* default case when we are not close to an edge */
4533 if( !lgDone )
4534 {
4535 CloudyFlux[j] = RebinSingleCell(BinLow,BinHigh,StarEner,StarFlux,StarPower,nCont);
4536 }
4537 }
4538
4539 FREE_CHECK( StarPower );
4540 FREE_SAFE( EdgeHigh );
4541 FREE_SAFE( EdgeLow );
4542 return;
4543}
4544
4546 realnum BinHigh,
4547 const realnum StarEner[], /* StarEner[nCont] */
4548 const realnum StarFlux[], /* StarFlux[nCont] */
4549 const realnum StarPower[], /* StarPower[nCont-1] */
4550 long nCont)
4551{
4552 long int i,
4553 ipHi,
4554 ipLo;
4555 double anu,
4556 retval,
4557 widflx;
4558 double sum,
4559 v1,
4560 val,
4561 x1,
4562 x2;
4563
4564 DEBUG_ENTRY( "RebinSingleCell()" );
4565
4566 /* >>chng 05 nov 22, use geometric mean instead of arithmetic mean, PvH */
4567 anu = sqrt(BinLow*BinHigh);
4568 /* >>chng 05 nov 22, reduce widflx if cell sticks out above highest frequency in model, PvH */
4569 widflx = MIN2(BinHigh,StarEner[nCont-1])-BinLow;
4570
4571 if( BinLow < StarEner[0] )
4572 {
4573 /* this is case where Cloudy's continuum is below stellar continuum,
4574 * (at least for part of the cell), so we do Rayleigh Jeans extrapolation */
4575 retval = (realnum)(StarFlux[0]*pow(anu/StarEner[0],2.));
4576 }
4577 else if( BinLow > StarEner[nCont-1] )
4578 {
4579 /* case where cloudy continuum is entirely above highest stellar point */
4580 retval = 0.0e00;
4581 }
4582 else
4583 {
4584 /* now go through stellar continuum to find bins corresponding to
4585 * this cloudy cell, stellar continuum defined through nCont cells */
4586 ipLo = RebinFind(StarEner,nCont,BinLow);
4587 ipHi = RebinFind(StarEner,nCont,BinHigh);
4588 /* sanity check */
4589 ASSERT( ipLo >= 0 && ipLo < nCont-1 && ipHi >= ipLo );
4590
4591 if( ipHi == ipLo )
4592 {
4593 /* Do the case where the cloudy cell and its edges are between
4594 * two adjacent stellar model points: do power-law interpolation */
4595 retval = (realnum)(StarFlux[ipLo]*pow(anu/StarEner[ipLo],(double)StarPower[ipLo]));
4596 }
4597 else
4598 {
4599 /* Do the case where the cloudy cell and its edges span two or more
4600 * stellar model cells: add segments with power-law interpolation up to
4601 * do the averaging.*/
4602
4603 sum = 0.;
4604
4605 /* ipHi points to stellar point at high end of cloudy continuum cell,
4606 * if the Cloudy cell extends beyond the stellar grid, ipHi == nCont-1
4607 * and the MIN2(ipHi,nCont-2) prevents access beyond allocated memory
4608 * ipLo points to low end, above we asserted that 0 <= ipLo < nCont-1 */
4609 for( i=ipLo; i <= MIN2(ipHi,nCont-2); i++ )
4610 {
4611 double pp1 = StarPower[i] + 1.;
4612
4613 if( i == ipLo )
4614 {
4615 x1 = BinLow;
4616 x2 = StarEner[i+1];
4617 v1 = StarFlux[i]*pow(x1/StarEner[i],(double)StarPower[i]);
4618 /*v2 = StarFlux[i+1];*/
4619 }
4620
4621 else if( i == ipHi )
4622 {
4623 x2 = BinHigh;
4624 x1 = StarEner[i];
4625 /*v2 = StarFlux[i]*pow(x2/StarEner[i],StarPower[i]);*/
4626 v1 = StarFlux[i];
4627 }
4628
4629 /*if( i > ipLo && i < ipHi )*/
4630 else
4631 {
4632 x1 = StarEner[i];
4633 x2 = StarEner[i+1];
4634 v1 = StarFlux[i];
4635 /*v2 = StarFlux[i+1];*/
4636 }
4637
4638 if( fabs(pp1) < 0.001 )
4639 {
4640 val = x1*v1*log(x2/x1);
4641 }
4642 else
4643 {
4644 val = pow(x2/x1,pp1) - 1.;
4645 val = val*x1*v1/pp1;
4646 }
4647 sum += val;
4648 }
4649
4650 retval = sum/widflx;
4651 }
4652 }
4653 return (realnum)retval;
4654}
4655
4656STATIC long RebinFind(const realnum array[], /* array[nArr] */
4657 long nArr,
4658 realnum val)
4659{
4660 long i1,
4661 i2,
4662 i3,
4663 ind = -2,
4664 sgn;
4665
4666 DEBUG_ENTRY( "RebinFind()" );
4667
4668 /* sanity check */
4669 ASSERT( nArr > 1 );
4670
4671 /* return ind(val) such that array[ind] <= val <= array[ind+1],
4672 *
4673 * NB NB: this routine assumes that array[] increases monotonically !
4674 *
4675 * the first two clauses indicate out-of-bounds conditions and
4676 * guarantee that when val1 <= val2, also ind(val1) <= ind(val2) */
4677
4678 if( val < array[0] )
4679 {
4680 ind = -1;
4681 }
4682 else if( val > array[nArr-1] )
4683 {
4684 ind = nArr-1;
4685 }
4686 else
4687 {
4688 /* do a binary search for ind */
4689 i1 = 0;
4690 i3 = nArr-1;
4691 while( i3-i1 > 1 )
4692 {
4693 i2 = (i1+i3)/2;
4694 sgn = sign3(val-array[i2]);
4695
4696 switch(sgn)
4697 {
4698 case -1:
4699 i3 = i2;
4700 break;
4701 case 0:
4702 ind = i2;
4703 return( ind );
4704 case 1:
4705 i1 = i2;
4706 break;
4707 }
4708 }
4709 ind = i1;
4710 }
4711
4712 /* sanity check */
4713 ASSERT( ind > -2 );
4714 return ind;
4715}
4716/*lint +e785 too few initializers */
4717/*lint +e801 use of go to depreciated */
static double x2[63]
static double x1[83]
t_called called
Definition called.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define POW3
Definition cddefines.h:936
#define REALLOC
Definition cddefines.h:519
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
int sign3(T x)
Definition cddefines.h:808
#define CALLOC
Definition cddefines.h:510
#define POW2
Definition cddefines.h:929
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
T * get_ptr(T *v)
Definition cddefines.h:1079
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
t_continuum continuum
Definition continuum.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const ios_base::openmode mode_r
Definition cpu.h:212
access_scheme
Definition cpu.h:207
@ AS_LOCAL_ONLY
Definition cpu.h:208
@ AS_DATA_OPTIONAL
Definition cpu.h:208
@ AS_DATA_ONLY_TRY
Definition cpu.h:207
@ AS_LOCAL_ONLY_TRY
Definition cpu.h:207
t_grid grid
Definition grid.cpp:5
t_optimize optimize
Definition optimize.cpp:5
UNUSED const double FR1RYD
Definition physconst.h:195
UNUSED const double PI
Definition physconst.h:29
UNUSED const double STEFAN_BOLTZ
Definition physconst.h:210
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double RYDLAM
Definition physconst.h:176
t_rfield rfield
Definition rfield.cpp:8
static const bool lgLINEAR
Definition stars.cpp:43
STATIC void InitGridCoStar(stellar_grid *)
Definition stars.cpp:3302
#define FREE_CHECK(PTR)
Definition stars.cpp:37
STATIC void InterpolateModel(const stellar_grid *, const double[], double[], const long[], const long[], long[], long, vector< realnum > &, IntStage)
Definition stars.cpp:3514
int Kurucz79Compile(process_counter &pc)
Definition stars.cpp:778
bool StarburstCompile(process_counter &pc)
Definition stars.cpp:1494
static const bool lgTAKELOG
Definition stars.cpp:44
bool StarburstInitialize(const char chInName[], const char chOutName[], sb_mode mode)
Definition stars.cpp:1291
STATIC void FindVCoStar(const stellar_grid *, double, realnum *, long[])
Definition stars.cpp:2290
STATIC void InterpolateRectGrid(const stellar_grid *, const double[], double *, double *)
Definition stars.cpp:3383
STATIC void FindHCoStar(const stellar_grid *, long, double, long, realnum *, long *, long *)
Definition stars.cpp:2232
STATIC bool lgValidBinFile(const char *, process_counter &, access_scheme)
Definition stars.cpp:3199
STATIC void InterpolateGridCoStar(const stellar_grid *, const double[], double *, double *)
Definition stars.cpp:2070
long RauchInterpolateHydr(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1201
STATIC void SetLimits(const stellar_grid *, double, const long[], const long[], const long[], const realnum[], double *, double *)
Definition stars.cpp:3848
static const long int VERSION_BIN
Definition stars.cpp:197
STATIC bool lgValidModel(const vector< Energy > &, const vector< realnum > &, double, double)
Definition stars.cpp:4394
long Kurucz79Interpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:799
static const bool lgSILENT
Definition stars.cpp:40
int WMBASICCompile(process_counter &pc)
Definition stars.cpp:1764
long WMBASICInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1788
STATIC void SetLimitsSub(const stellar_grid *, double, const long[], const long[], long[], long, double *, double *)
Definition stars.cpp:3944
STATIC void RebinAtmosphere(long, const realnum[], const realnum[], realnum[], long, const realnum[])
Definition stars.cpp:4424
int AtlasCompile(process_counter &pc)
Definition stars.cpp:394
STATIC void CoStarListModels(const stellar_grid *)
Definition stars.cpp:2355
long AtlasInterpolate(double val[], long *nval, long *ndim, const char *chMetalicity, const char *chODFNew, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:515
int RauchCompile(process_counter &pc)
Definition stars.cpp:879
STATIC void FindIndex(const double[], long, double, long *, long *, bool *)
Definition stars.cpp:4253
static const int NMODS_HpHE
Definition stars.cpp:32
bool GridCompile(const char *InName)
Definition stars.cpp:693
long WernerInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1708
STATIC int RauchInitializeSub(const char[], const char[], const vector< mpp > &, long, long, long, const double[], int)
Definition stars.cpp:2399
void getdataline(fstream &, string &)
Definition stars.cpp:2710
STATIC void FillJ(const stellar_grid *, long[], double[], long, bool)
Definition stars.cpp:4093
STATIC bool lgValidAsciiFile(const char *, access_scheme)
Definition stars.cpp:3279
IntStage
Definition stars.cpp:46
@ IS_FIRST
Definition stars.cpp:47
@ IS_SECOND
Definition stars.cpp:47
@ IS_UNDEFINED
Definition stars.cpp:47
static const int NMODS_HYDR
Definition stars.cpp:28
int MihalasCompile(process_counter &pc)
Definition stars.cpp:829
STATIC void GetBins(const stellar_grid *, vector< Energy > &)
Definition stars.cpp:3731
STATIC bool lgCompileAtmosphereCoStar(const char[], const char[], const realnum[], long, process_counter &)
Definition stars.cpp:1819
static const int NMODS_PG1159
Definition stars.cpp:26
static const long int VERSION_ASCII
Definition stars.cpp:192
#define FREE_SAFE(PTR)
Definition stars.cpp:38
long RauchInterpolateHNi(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1107
static const int NRAUCH
Definition stars.cpp:20
STATIC bool lgFileReadable(const char *, process_counter &, access_scheme)
Definition stars.cpp:4329
STATIC void ValidateGrid(const stellar_grid *, double)
Definition stars.cpp:4349
STATIC void RauchReadMPP(vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &)
Definition stars.cpp:2612
int CoStarCompile(process_counter &pc)
Definition stars.cpp:586
long CoStarInterpolate(double val[], long *nval, long *ndim, IntMode imode, bool lgHalo, bool lgList, double *val0_lo, double *val0_hi)
Definition stars.cpp:623
static const long int VERSION_RAUCH_MPP
Definition stars.cpp:199
STATIC void InitGrid(stellar_grid *, bool)
Definition stars.cpp:3059
static const int NMODS_HCA
Definition stars.cpp:22
static const int MNTS
Definition stars.cpp:17
static const int NMODS_HNI
Definition stars.cpp:24
static const int NSB99
Definition stars.cpp:15
STATIC long JIndex(const stellar_grid *, const long[])
Definition stars.cpp:4167
long GridInterpolate(double val[], long *nval, long *ndim, const char *FileName, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:738
STATIC bool lgCompileAtmosphere(const char[], const char[], const realnum[], long, process_counter &)
Definition stars.cpp:2723
long RauchInterpolateHpHe(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1261
int TlustyCompile(process_counter &pc)
Definition stars.cpp:1518
long RauchInterpolateHelium(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1231
long TlustyInterpolate(double val[], long *nval, long *ndim, tl_grid tlg, const char *chMetalicity, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1590
STATIC realnum RebinSingleCell(realnum, realnum, const realnum[], const realnum[], const realnum[], long)
Definition stars.cpp:4545
long RauchInterpolateHCa(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1073
int WernerCompile(process_counter &pc)
Definition stars.cpp:1653
STATIC void CheckVal(const stellar_grid *, double[], long *, long *)
Definition stars.cpp:3354
static const int NMODS_HELIUM
Definition stars.cpp:30
STATIC void InterpolateModelCoStar(const stellar_grid *, const double[], double[], const long[], const long[], long[], long, long, vector< realnum > &)
Definition stars.cpp:3659
STATIC void GetModel(const stellar_grid *, long, vector< realnum > &, bool, bool)
Definition stars.cpp:3762
void AtmospheresAvail(void)
Definition stars.cpp:202
long RauchInterpolateCOWD(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1171
static const bool lgVERBOSE
Definition stars.cpp:41
STATIC void SearchModel(const mpp[], bool, long, const double[], long, long *, long *)
Definition stars.cpp:4184
STATIC void InitIndexArrays(stellar_grid *, bool)
Definition stars.cpp:4022
long MihalasInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:849
STATIC long RebinFind(const realnum[], long, realnum)
Definition stars.cpp:4656
STATIC void FreeGrid(stellar_grid *)
Definition stars.cpp:3493
long RauchInterpolatePG1159(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition stars.cpp:1141
IntMode
Definition stars.h:16
@ IM_COSTAR_MZAMS_AGE
Definition stars.h:18
@ IM_RECT_GRID
Definition stars.h:17
@ IM_COSTAR_AGE_MZAMS
Definition stars.h:18
@ IM_COSTAR_TEFF_LOGG
Definition stars.h:18
@ IM_COSTAR_TEFF_MODID
Definition stars.h:17
sb_mode
Definition stars.h:25
@ SB_STELLAR
Definition stars.h:26
@ SB_NEBULAR
Definition stars.h:26
@ SB_TOTAL
Definition stars.h:26
tl_grid
Definition stars.h:21
@ TL_OBSTAR
Definition stars.h:22
@ TL_OSTAR
Definition stars.h:22
@ TL_BSTAR
Definition stars.h:22
#define MNAM
Definition stars.h:12
#define MDIM
Definition stars.h:9
Definition stars.cpp:52
char chGrid
Definition stars.cpp:55
double par[MDIM]
Definition stars.cpp:53
int modid
Definition stars.cpp:54
int notProcessed
Definition stars.h:32
const char * command
Definition stars.cpp:109
long * nval
Definition stars.cpp:130
mpp * telg
Definition stars.cpp:126
bool lgIsTeffLoggGrid
Definition stars.cpp:100
const char * ident
Definition stars.cpp:107
int32 ngrid
Definition stars.cpp:119
IntMode imode
Definition stars.cpp:111
FILE * ioIN
Definition stars.cpp:104
long * jval
Definition stars.cpp:147
uint32 nOffset
Definition stars.cpp:121
long * jlo
Definition stars.cpp:138
uint32 nBlocksize
Definition stars.cpp:123
int32 ndim
Definition stars.cpp:113
long * jhi
Definition stars.cpp:139
access_scheme scheme
Definition stars.cpp:102
string name
Definition stars.cpp:98
int32 nmods
Definition stars.cpp:117
long * trackLen
Definition stars.cpp:143
int32 npar
Definition stars.cpp:115
char names[MDIM][MNAM+1]
Definition stars.cpp:141
double ** val
Definition stars.cpp:128
long nTracks
Definition stars.cpp:145
static const unsigned int NMD5
Definition thirdparty.h:384