cloudy  trunk
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 
15 static const int NSB99 = 1250;
17 static const int MNTS = 200;
18 
20 static const int NRAUCH = 19951;
22 static const int NMODS_HCA = 66;
24 static const int NMODS_HNI = 51;
26 static const int NMODS_PG1159 = 71;
28 static const int NMODS_HYDR = 100;
30 static const int NMODS_HELIUM = 81;
32 static 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 
40 static const bool lgSILENT = false;
41 static const bool lgVERBOSE = true;
42 
43 static const bool lgLINEAR = false;
44 static const bool lgTAKELOG = true;
45 
46 typedef enum {
47  IS_UNDEFINED, IS_FIRST, IS_SECOND
48 } IntStage;
49 
51 typedef struct
52 {
53  double par[MDIM];
54  int modid;
55  char chGrid;
56 } mpp;
57 
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 
95 typedef struct
96 {
98  string name;
100  bool lgIsTeffLoggGrid;
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;
123  uint32 nBlocksize;
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] */
145  long nTracks;
147  long *jval;
148 } stellar_grid;
149 
150 /* internal routines */
151 STATIC bool lgCompileAtmosphereCoStar(const char[],const char[],const realnum[],long,process_counter&);
152 STATIC void InterpolateGridCoStar(const stellar_grid*,const double[],double*,double*);
153 STATIC void FindHCoStar(const stellar_grid*,long,double,long,realnum*,long*,long*);
154 STATIC void FindVCoStar(const stellar_grid*,double,realnum*,long[]);
155 STATIC void CoStarListModels(const stellar_grid*);
156 STATIC int RauchInitializeSub(const char[],const char[],const vector<mpp>&,long,long,
157  long,const double[],int);
158 STATIC void RauchReadMPP(vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&,vector<mpp>&);
159 inline void getdataline(fstream&,string&);
160 STATIC bool lgCompileAtmosphere(const char[],const char[],const realnum[],long,process_counter&);
161 STATIC void InitGrid(stellar_grid*,bool);
162 STATIC bool lgValidBinFile(const char*,process_counter&,access_scheme);
163 STATIC bool lgValidAsciiFile(const char*,access_scheme);
164 STATIC void InitGridCoStar(stellar_grid*);
165 STATIC void CheckVal(const stellar_grid*,double[],long*,long*);
166 STATIC void InterpolateRectGrid(const stellar_grid*,const double[],double*,double*);
167 STATIC void FreeGrid(stellar_grid*);
168 STATIC void InterpolateModel(const stellar_grid*,const double[],double[],const long[],
169  const long[],long[],long,vector<realnum>&,IntStage);
170 STATIC void InterpolateModelCoStar(const stellar_grid*,const double[],double[],const long[],
171  const long[],long[],long,long,vector<realnum>&);
172 STATIC void GetBins(const stellar_grid*,vector<Energy>&);
173 STATIC void GetModel(const stellar_grid*,long,vector<realnum>&,bool,bool);
174 STATIC void SetLimits(const stellar_grid*,double,const long[],const long[],const long[],
175  const realnum[],double*,double*);
176 STATIC void SetLimitsSub(const stellar_grid*,double,const long[],const long[],long[],long,
177  double*,double*);
178 STATIC void InitIndexArrays(stellar_grid*,bool);
179 STATIC void FillJ(const stellar_grid*,long[],double[],long,bool);
180 STATIC long JIndex(const stellar_grid*,const long[]);
181 STATIC void SearchModel(const mpp[],bool,long,const double[],long,long*,long*);
182 STATIC void FindIndex(const double[],long,double,long*,long*,bool*);
183 STATIC bool lgFileReadable(const char*, process_counter&,access_scheme);
184 STATIC void ValidateGrid(const stellar_grid*,double);
185 STATIC bool lgValidModel(const vector<Energy>&,const vector<realnum>&,double,double);
186 STATIC void RebinAtmosphere(long,const realnum[],const realnum[],realnum[],long,const realnum[]);
187 STATIC realnum RebinSingleCell(realnum,realnum,const realnum[],const realnum[],const realnum[],long);
188 STATIC long RebinFind(const realnum[],long,realnum);
189 
190 
191 /* the version number for the ascii/binary atmosphere files */
192 static const long int VERSION_ASCII = 20060612L;
193 /* binary files are incompatible when floats are converted to doubles */
194 #ifdef FLT_IS_DBL
195 static const long int VERSION_BIN = 201009020L;
196 #else
197 static const long int VERSION_BIN = 201009021L;
198 #endif
199 static const long int VERSION_RAUCH_MPP = 20090324;
200 
202 void AtmospheresAvail( void )
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... */
220  access_scheme as = AS_DATA_ONLY_TRY;
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) */
394 int AtlasCompile(process_counter& pc)
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 
429  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
515 long 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];
525  stellar_grid grid;
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*/
586 int CoStarCompile(process_counter& pc)
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 
611  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
623 long 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 {
632  stellar_grid grid;
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 */
647  InitGridCoStar( &grid );
648  /* override default interpolation mode */
649  grid.imode = imode;
650 
651  if( lgList )
652  {
653  CoStarListModels( &grid );
654  cdEXIT(EXIT_SUCCESS);
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 */
693 bool 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  {
713  stellar_grid grid;
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 */
738 long 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];
747  stellar_grid grid;
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 */
778 int Kurucz79Compile(process_counter& pc)
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 
791  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
799 long Kurucz79Interpolate(double val[], /* val[nval] */
800  long *nval,
801  long *ndim,
802  bool lgList,
803  double *Tlow,
804  double *Thigh)
805 {
806  stellar_grid grid;
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 */
829 int MihalasCompile(process_counter& pc)
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 
841  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
849 long MihalasInterpolate(double val[], /* val[nval] */
850  long *nval,
851  long *ndim,
852  bool lgList,
853  double *Tlow,
854  double *Thigh)
855 {
856  stellar_grid grid;
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 */
879 int RauchCompile(process_counter& pc)
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;
914  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
1073 long RauchInterpolateHCa(double val[], /* val[nval] */
1074  long *nval,
1075  long *ndim,
1076  bool lgHalo,
1077  bool lgList,
1078  double *Tlow,
1079  double *Thigh)
1080 {
1081  stellar_grid grid;
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 */
1107 long RauchInterpolateHNi(double val[], /* val[nval] */
1108  long *nval,
1109  long *ndim,
1110  bool lgHalo,
1111  bool lgList,
1112  double *Tlow,
1113  double *Thigh)
1114 {
1115  stellar_grid grid;
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 */
1141 long RauchInterpolatePG1159(double val[], /* val[nval] */
1142  long *nval,
1143  long *ndim,
1144  bool lgList,
1145  double *Tlow,
1146  double *Thigh)
1147 {
1148  stellar_grid grid;
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 */
1171 long RauchInterpolateCOWD(double val[], /* val[nval] */
1172  long *nval,
1173  long *ndim,
1174  bool lgList,
1175  double *Tlow,
1176  double *Thigh)
1177 {
1178  stellar_grid grid;
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 */
1201 long RauchInterpolateHydr(double val[], /* val[nval] */
1202  long *nval,
1203  long *ndim,
1204  bool lgList,
1205  double *Tlow,
1206  double *Thigh)
1207 {
1208  stellar_grid grid;
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 */
1231 long RauchInterpolateHelium(double val[], /* val[nval] */
1232  long *nval,
1233  long *ndim,
1234  bool lgList,
1235  double *Tlow,
1236  double *Thigh)
1237 {
1238  stellar_grid grid;
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 */
1261 long RauchInterpolateHpHe(double val[], /* val[nval] */
1262  long *nval,
1263  long *ndim,
1264  bool lgList,
1265  double *Tlow,
1266  double *Thigh)
1267 {
1268  stellar_grid grid;
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 */
1291 bool 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 
1486 error:
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 */
1494 bool StarburstCompile(process_counter& pc)
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;
1505  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
1518 int TlustyCompile(process_counter& pc)
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 
1529  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
1590 long 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];
1600  stellar_grid grid;
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) */
1653 int WernerCompile(process_counter& pc)
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 
1700  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
1708 long WernerInterpolate(double val[], /* val[nval] */
1709  long *nval,
1710  long *ndim,
1711  bool lgList,
1712  double *Tlow,
1713  double *Thigh)
1714 {
1715  stellar_grid grid;
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 */
1764 int WMBASICCompile(process_counter& pc)
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 
1780  access_scheme as = AS_LOCAL_ONLY_TRY;
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 */
1788 long WMBASICInterpolate(double val[], /* val[nval] */
1789  long *nval,
1790  long *ndim,
1791  bool lgList,
1792  double *Tlow,
1793  double *Thigh)
1794 {
1795  stellar_grid grid;
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 */
1819 STATIC 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 
2060 error:
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 */
2070 STATIC 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  {
2088  case IM_COSTAR_TEFF_MODID:
2089  case IM_COSTAR_TEFF_LOGG:
2090  lval[0] = val[0];
2091  lval[1] = val[1];
2092  off = 0;
2093  break;
2094  case IM_COSTAR_MZAMS_AGE:
2095  lval[0] = log10(val[0]); /* use log10(M_ZAMS) internally */
2096  lval[1] = val[1];
2097  off = 2;
2098  break;
2099  case IM_COSTAR_AGE_MZAMS:
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 );
2107  cdEXIT(EXIT_FAILURE);
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" );
2173  cdEXIT(EXIT_FAILURE);
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 */
2232 STATIC void FindHCoStar(const stellar_grid *grid,
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 */
2290 STATIC void FindVCoStar(const stellar_grid *grid,
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 */
2355 STATIC void CoStarListModels(const stellar_grid *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 */
2399 STATIC 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();
2488  cdEXIT(EXIT_FAILURE);
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 
2606 error:
2607  FREE_CHECK( fluxes );
2608  FREE_CHECK( wavl );
2609  return 1;
2610 }
2611 
2612 STATIC 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 );
2636  cdEXIT(EXIT_FAILURE);
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 
2710 inline 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) */
2723 STATIC 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 );
2942  cdEXIT(EXIT_FAILURE);
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" );
2971  cdEXIT(EXIT_FAILURE);
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 );
2999  cdEXIT(EXIT_FAILURE);
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 
3049 error:
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 
3059 STATIC void InitGrid(stellar_grid *grid,
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");
3081  cdEXIT(EXIT_FAILURE);
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" );
3100  cdEXIT(EXIT_FAILURE);
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 );
3110  cdEXIT(EXIT_FAILURE);
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 );
3119  cdEXIT(EXIT_FAILURE);
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 );
3130  cdEXIT(EXIT_FAILURE);
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" );
3145  cdEXIT(EXIT_FAILURE);
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" );
3163  cdEXIT(EXIT_FAILURE);
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 );
3182  cdEXIT(EXIT_FAILURE);
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 */
3199 STATIC 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];
3204  stellar_grid grid;
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 */
3279 STATIC 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[] */
3302 STATIC 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 
3354 STATIC void CheckVal(const stellar_grid *grid,
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 );
3373  cdEXIT(EXIT_FAILURE);
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 );
3379  cdEXIT(EXIT_FAILURE);
3380  }
3381 }
3382 
3383 STATIC void InterpolateRectGrid(const stellar_grid *grid,
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] );
3425  cdEXIT(EXIT_FAILURE);
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 
3493 STATIC void FreeGrid(stellar_grid *grid)
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 
3514 STATIC void InterpolateModel(const stellar_grid *grid,
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" );
3553  cdEXIT(EXIT_FAILURE);
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 
3659 STATIC void InterpolateModelCoStar(const stellar_grid *grid,
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 
3731 STATIC void GetBins(const stellar_grid *grid,
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");
3746  cdEXIT(EXIT_FAILURE);
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" );
3753  cdEXIT(EXIT_FAILURE);
3754  }
3755 
3756  for( long i=0; i < rfield.nupper; ++i )
3757  ener[i].set(data[i]);
3758 
3759  return;
3760 }
3761 
3762 STATIC void GetModel(const stellar_grid *grid,
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 );
3785  cdEXIT(EXIT_FAILURE);
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 );
3791  cdEXIT(EXIT_FAILURE);
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 
3848 STATIC void SetLimits(const stellar_grid *grid,
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;
3875  case IM_COSTAR_TEFF_MODID:
3876  case IM_COSTAR_TEFF_LOGG:
3877  case IM_COSTAR_MZAMS_AGE:
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;
3890  case IM_COSTAR_AGE_MZAMS:
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 );
3915  cdEXIT(EXIT_FAILURE);
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 );
3925  cdEXIT(EXIT_FAILURE);
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 
3944 STATIC void SetLimitsSub(const stellar_grid *grid,
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 
4022 STATIC void InitIndexArrays(stellar_grid *grid,
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 )
4089  cdEXIT(EXIT_SUCCESS);
4090  return;
4091 }
4092 
4093 STATIC void FillJ(const stellar_grid *grid,
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 
4167 STATIC long JIndex(const stellar_grid *grid,
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 
4184 STATIC 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 
4253 STATIC 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();
4326  cdEXIT(EXIT_FAILURE);
4327 }
4328 
4329 STATIC 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 */
4349 STATIC void ValidateGrid(const stellar_grid *grid,
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 
4394 STATIC 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 */
4424 STATIC 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 
4545 STATIC realnum RebinSingleCell(realnum BinLow,
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 
4656 STATIC 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 */
RauchInterpolateHydr
long RauchInterpolateHydr(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1201
TL_BSTAR
@ TL_BSTAR
Definition: stars.h:22
FREE_CHECK
#define FREE_CHECK(PTR)
Definition: stars.cpp:37
FR1RYD
const UNUSED double FR1RYD
Definition: physconst.h:195
lgValidModel
STATIC bool lgValidModel(const vector< Energy > &, const vector< realnum > &, double, double)
Definition: stars.cpp:4394
Kurucz79Interpolate
long Kurucz79Interpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:799
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
RauchReadMPP
STATIC void RauchReadMPP(vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &, vector< mpp > &)
Definition: stars.cpp:2612
IM_RECT_GRID
@ IM_RECT_GRID
Definition: stars.h:17
stellar_grid::ngrid
int32 ngrid
Definition: stars.cpp:119
stellar_grid::name
string name
Definition: stars.cpp:98
RauchInterpolateHCa
long RauchInterpolateHCa(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1073
CoStarListModels
STATIC void CoStarListModels(const stellar_grid *)
Definition: stars.cpp:2355
stellar_grid::jval
long * jval
Definition: stars.cpp:147
FindIndex
STATIC void FindIndex(const double[], long, double, long *, long *, bool *)
Definition: stars.cpp:4253
NMODS_PG1159
static const int NMODS_PG1159
Definition: stars.cpp:26
GridCompile
bool GridCompile(const char *InName)
Definition: stars.cpp:693
rfield
t_rfield rfield
Definition: rfield.cpp:8
lgFileReadable
STATIC bool lgFileReadable(const char *, process_counter &, access_scheme)
Definition: stars.cpp:4329
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_rfield::tNu
vector< Energy > tNu[LIMSPC]
Definition: rfield.h:330
STEFAN_BOLTZ
const UNUSED double STEFAN_BOLTZ
Definition: physconst.h:210
IM_COSTAR_TEFF_LOGG
@ IM_COSTAR_TEFF_LOGG
Definition: stars.h:18
lgValidAsciiFile
STATIC bool lgValidAsciiFile(const char *, access_scheme)
Definition: stars.cpp:3279
t_continuum::ResolutionScaleFactor
double ResolutionScaleFactor
Definition: continuum.h:90
stellar_grid::trackLen
long * trackLen
Definition: stars.cpp:143
mpp::par
double par[MDIM]
Definition: stars.cpp:53
realnum
float realnum
Definition: cddefines.h:103
rfield.h
NMODS_HpHE
static const int NMODS_HpHE
Definition: stars.cpp:32
STATIC
#define STATIC
Definition: cddefines.h:97
SB_STELLAR
@ SB_STELLAR
Definition: stars.h:26
TlustyCompile
int TlustyCompile(process_counter &pc)
Definition: stars.cpp:1518
grid
t_grid grid
Definition: grid.cpp:5
InterpolateModelCoStar
STATIC void InterpolateModelCoStar(const stellar_grid *, const double[], double[], const long[], const long[], long[], long, long, vector< realnum > &)
Definition: stars.cpp:3659
AS_DATA_ONLY_TRY
@ AS_DATA_ONLY_TRY
Definition: cpu.h:207
AS_LOCAL_ONLY
@ AS_LOCAL_ONLY
Definition: cpu.h:208
t_optimize::lgVarOn
bool lgVarOn
Definition: optimize.h:203
IS_UNDEFINED
@ IS_UNDEFINED
Definition: stars.cpp:47
MNAM
#define MNAM
Definition: stars.h:12
thirdparty.h
process_counter::nFound
int nFound
Definition: stars.h:31
mpp::chGrid
char chGrid
Definition: stars.cpp:55
t_rfield::RSFCheck
double RSFCheck[LIMSPC]
Definition: rfield.h:339
cloudy_exit
Definition: cddefines.h:396
NMODS_HYDR
static const int NMODS_HYDR
Definition: stars.cpp:28
stellar_grid::jhi
long * jhi
Definition: stars.cpp:139
stars.h
stellar_grid::ndim
int32 ndim
Definition: stars.cpp:113
WMBASICInterpolate
long WMBASICInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1788
NSB99
static const int NSB99
Definition: stars.cpp:15
t_rfield::lgContMalloc
bool lgContMalloc[LIMSPC]
Definition: rfield.h:343
stellar_grid::ident
const char * ident
Definition: stars.cpp:107
IntStage
IntStage
Definition: stars.cpp:46
MDIM
#define MDIM
Definition: stars.h:9
InitGridCoStar
STATIC void InitGridCoStar(stellar_grid *)
Definition: stars.cpp:3302
stellar_grid::ioIN
FILE * ioIN
Definition: stars.cpp:104
AS_DATA_OPTIONAL
@ AS_DATA_OPTIONAL
Definition: cpu.h:208
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
stellar_grid::telg
mpp * telg
Definition: stars.cpp:126
GetModel
STATIC void GetModel(const stellar_grid *, long, vector< realnum > &, bool, bool)
Definition: stars.cpp:3762
t_continuum::mesh_md5sum
string mesh_md5sum
Definition: continuum.h:127
MihalasInterpolate
long MihalasInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:849
TlustyInterpolate
long TlustyInterpolate(double val[], long *nval, long *ndim, tl_grid tlg, const char *chMetalicity, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1590
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
AtlasCompile
int AtlasCompile(process_counter &pc)
Definition: stars.cpp:394
mpp::modid
int modid
Definition: stars.cpp:54
t_rfield::tslop
vector< realnum > tslop[LIMSPC]
Definition: rfield.h:331
x1
static double x1[83]
Definition: atmdat_3body.cpp:28
POW2
#define POW2
Definition: cddefines.h:929
NRAUCH
static const int NRAUCH
Definition: stars.cpp:20
process_counter::notProcessed
int notProcessed
Definition: stars.h:32
MIN2
#define MIN2
Definition: cddefines.h:761
MNTS
static const int MNTS
Definition: stars.cpp:17
CoStarInterpolate
long CoStarInterpolate(double val[], long *nval, long *ndim, IntMode imode, bool lgHalo, bool lgList, double *val0_lo, double *val0_hi)
Definition: stars.cpp:623
lgValidBinFile
STATIC bool lgValidBinFile(const char *, process_counter &, access_scheme)
Definition: stars.cpp:3199
StarburstCompile
bool StarburstCompile(process_counter &pc)
Definition: stars.cpp:1494
NMD5
static const unsigned int NMD5
Definition: thirdparty.h:384
InterpolateGridCoStar
STATIC void InterpolateGridCoStar(const stellar_grid *, const double[], double *, double *)
Definition: stars.cpp:2070
SetLimits
STATIC void SetLimits(const stellar_grid *, double, const long[], const long[], const long[], const realnum[], double *, double *)
Definition: stars.cpp:3848
PI
const UNUSED double PI
Definition: physconst.h:29
stellar_grid::nOffset
uint32 nOffset
Definition: stars.cpp:121
REALLOC
#define REALLOC
Definition: cddefines.h:519
mode_r
const ios_base::openmode mode_r
Definition: cpu.h:212
IM_COSTAR_TEFF_MODID
@ IM_COSTAR_TEFF_MODID
Definition: stars.h:17
Kurucz79Compile
int Kurucz79Compile(process_counter &pc)
Definition: stars.cpp:778
lgCompileAtmosphere
STATIC bool lgCompileAtmosphere(const char[], const char[], const realnum[], long, process_counter &)
Definition: stars.cpp:2723
CheckVal
STATIC void CheckVal(const stellar_grid *, double[], long *, long *)
Definition: stars.cpp:3354
optimize
t_optimize optimize
Definition: optimize.cpp:5
RauchCompile
int RauchCompile(process_counter &pc)
Definition: stars.cpp:879
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_rfield::egamry
realnum egamry
Definition: rfield.h:52
VERSION_RAUCH_MPP
static const long int VERSION_RAUCH_MPP
Definition: stars.cpp:199
stellar_grid::nTracks
long nTracks
Definition: stars.cpp:145
InitIndexArrays
STATIC void InitIndexArrays(stellar_grid *, bool)
Definition: stars.cpp:4022
cddefines.h
lgCompileAtmosphereCoStar
STATIC bool lgCompileAtmosphereCoStar(const char[], const char[], const realnum[], long, process_counter &)
Definition: stars.cpp:1819
IM_COSTAR_AGE_MZAMS
@ IM_COSTAR_AGE_MZAMS
Definition: stars.h:18
FREE_SAFE
#define FREE_SAFE(PTR)
Definition: stars.cpp:38
IS_FIRST
@ IS_FIRST
Definition: stars.cpp:47
InterpolateModel
STATIC void InterpolateModel(const stellar_grid *, const double[], double[], const long[], const long[], long[], long, vector< realnum > &, IntStage)
Definition: stars.cpp:3514
process_counter::nFail
int nFail
Definition: stars.h:34
stellar_grid
Definition: stars.cpp:95
POW3
#define POW3
Definition: cddefines.h:936
GetBins
STATIC void GetBins(const stellar_grid *, vector< Energy > &)
Definition: stars.cpp:3731
optimize.h
stellar_grid::command
const char * command
Definition: stars.cpp:109
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
getdataline
void getdataline(fstream &, string &)
Definition: stars.cpp:2710
SB_NEBULAR
@ SB_NEBULAR
Definition: stars.h:26
t_called::lgTalk
bool lgTalk
Definition: called.h:12
VERSION_ASCII
static const long int VERSION_ASCII
Definition: stars.cpp:192
t_rfield::AnuOrg
double * AnuOrg
Definition: rfield.h:62
access_scheme
access_scheme
Definition: cpu.h:207
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
InterpolateRectGrid
STATIC void InterpolateRectGrid(const stellar_grid *, const double[], double *, double *)
Definition: stars.cpp:3383
stellar_grid::jlo
long * jlo
Definition: stars.cpp:138
SPEEDLIGHT
const UNUSED double SPEEDLIGHT
Definition: physconst.h:100
stellar_grid::nBlocksize
uint32 nBlocksize
Definition: stars.cpp:123
stellar_grid::scheme
access_scheme scheme
Definition: stars.cpp:102
FindHCoStar
STATIC void FindHCoStar(const stellar_grid *, long, double, long, realnum *, long *, long *)
Definition: stars.cpp:2232
MAX2
#define MAX2
Definition: cddefines.h:782
RauchInterpolateHpHe
long RauchInterpolateHpHe(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1261
RYDLAM
const UNUSED double RYDLAM
Definition: physconst.h:176
RauchInterpolateCOWD
long RauchInterpolateCOWD(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1171
FillJ
STATIC void FillJ(const stellar_grid *, long[], double[], long, bool)
Definition: stars.cpp:4093
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
sb_mode
sb_mode
Definition: stars.h:25
NMODS_HELIUM
static const int NMODS_HELIUM
Definition: stars.cpp:30
x2
static double x2[63]
Definition: atmdat_3body.cpp:20
t_rfield::nupper
long int nupper
Definition: rfield.h:46
TL_OBSTAR
@ TL_OBSTAR
Definition: stars.h:22
RebinAtmosphere
STATIC void RebinAtmosphere(long, const realnum[], const realnum[], realnum[], long, const realnum[])
Definition: stars.cpp:4424
stellar_grid::lgIsTeffLoggGrid
bool lgIsTeffLoggGrid
Definition: stars.cpp:100
IM_COSTAR_MZAMS_AGE
@ IM_COSTAR_MZAMS_AGE
Definition: stars.h:18
CoStarCompile
int CoStarCompile(process_counter &pc)
Definition: stars.cpp:586
t_continuum::lgCoStarInterpolationCaution
bool lgCoStarInterpolationCaution
Definition: continuum.h:94
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
SB_TOTAL
@ SB_TOTAL
Definition: stars.h:26
IS_SECOND
@ IS_SECOND
Definition: stars.cpp:47
t_rfield::anu
double * anu
Definition: rfield.h:58
JIndex
STATIC long JIndex(const stellar_grid *, const long[])
Definition: stars.cpp:4167
TL_OSTAR
@ TL_OSTAR
Definition: stars.h:22
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
RauchInitializeSub
STATIC int RauchInitializeSub(const char[], const char[], const vector< mpp > &, long, long, long, const double[], int)
Definition: stars.cpp:2399
InitGrid
STATIC void InitGrid(stellar_grid *, bool)
Definition: stars.cpp:3059
physconst.h
lgVERBOSE
static const bool lgVERBOSE
Definition: stars.cpp:41
lgSILENT
static const bool lgSILENT
Definition: stars.cpp:40
called
t_called called
Definition: called.cpp:5
NMODS_HCA
static const int NMODS_HCA
Definition: stars.cpp:22
AtmospheresAvail
void AtmospheresAvail(void)
Definition: stars.cpp:202
PrintEfmt
#define PrintEfmt(F, V)
Definition: cddefines.h:1472
get_ptr
T * get_ptr(T *v)
Definition: cddefines.h:1079
SetLimitsSub
STATIC void SetLimitsSub(const stellar_grid *, double, const long[], const long[], long[], long, double *, double *)
Definition: stars.cpp:3944
tl_grid
tl_grid
Definition: stars.h:21
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
StarburstInitialize
bool StarburstInitialize(const char chInName[], const char chOutName[], sb_mode mode)
Definition: stars.cpp:1291
MihalasCompile
int MihalasCompile(process_counter &pc)
Definition: stars.cpp:829
t_rfield::nShape
long int nShape
Definition: rfield.h:322
FreeGrid
STATIC void FreeGrid(stellar_grid *)
Definition: stars.cpp:3493
stellar_grid::imode
IntMode imode
Definition: stars.cpp:111
RauchInterpolateHelium
long RauchInterpolateHelium(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1231
AtlasInterpolate
long AtlasInterpolate(double val[], long *nval, long *ndim, const char *chMetalicity, const char *chODFNew, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:515
ShowMe
void ShowMe(void)
Definition: service.cpp:181
RauchInterpolatePG1159
long RauchInterpolatePG1159(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1141
VERSION_BIN
static const long int VERSION_BIN
Definition: stars.cpp:197
WernerCompile
int WernerCompile(process_counter &pc)
Definition: stars.cpp:1653
continuum
t_continuum continuum
Definition: continuum.cpp:5
process_counter::nOK
int nOK
Definition: stars.h:33
t_rfield::emm
realnum emm
Definition: rfield.h:49
stellar_grid::val
double ** val
Definition: stars.cpp:128
ValidateGrid
STATIC void ValidateGrid(const stellar_grid *, double)
Definition: stars.cpp:4349
RebinSingleCell
STATIC realnum RebinSingleCell(realnum, realnum, const realnum[], const realnum[], const realnum[], long)
Definition: stars.cpp:4545
lgTAKELOG
static const bool lgTAKELOG
Definition: stars.cpp:44
called.h
IntMode
IntMode
Definition: stars.h:16
sign3
int sign3(T x)
Definition: cddefines.h:808
continuum.h
SearchModel
STATIC void SearchModel(const mpp[], bool, long, const double[], long, long *, long *)
Definition: stars.cpp:4184
WernerInterpolate
long WernerInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1708
NMODS_HNI
static const int NMODS_HNI
Definition: stars.cpp:24
stellar_grid::npar
int32 npar
Definition: stars.cpp:115
stellar_grid::nval
long * nval
Definition: stars.cpp:130
CALLOC
#define CALLOC
Definition: cddefines.h:510
lgLINEAR
static const bool lgLINEAR
Definition: stars.cpp:43
GridInterpolate
long GridInterpolate(double val[], long *nval, long *ndim, const char *FileName, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:738
stellar_grid::nmods
int32 nmods
Definition: stars.cpp:117
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
mpp
Definition: stars.cpp:51
RebinFind
STATIC long RebinFind(const realnum[], long, realnum)
Definition: stars.cpp:4656
process_counter
Definition: stars.h:29
WMBASICCompile
int WMBASICCompile(process_counter &pc)
Definition: stars.cpp:1764
FindVCoStar
STATIC void FindVCoStar(const stellar_grid *, double, realnum *, long[])
Definition: stars.cpp:2290
AS_LOCAL_ONLY_TRY
@ AS_LOCAL_ONLY_TRY
Definition: cpu.h:207
RauchInterpolateHNi
long RauchInterpolateHNi(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1107
Generated by   doxygen 1.8.17