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