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