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