cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_table.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ParseTable parse the table read command */
4 /*lines_table invoked by table lines command, check if we can find all lines in a given list */
5 #include "cddefines.h"
6 #include "cddrive.h"
7 #include "optimize.h"
8 #include "rfield.h"
9 #include "trace.h"
10 #include "lines.h"
11 #include "radius.h"
12 #include "input.h"
13 #include "stars.h"
14 #include "prt.h"
15 #include "parser.h"
16 #include "save.h"
17 
18 /*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
19 STATIC void ReadTable(const char * fnam);
20 
21 /* these will become the label and wl for a possible list of lines,
22  * obtained when tables lines used */
23 static string chLINE_LIST;
24 
25 /* Black's ISM continuum, with He hole filled in */
26 static const int NISM = 23;
27 static double tnuism[NISM],
28  fnuism[NISM];
29 
30 /* z=2 background,
31  * >>refer continuum background Haardt, Francesco, & Madau, Piero, 1996,
32  * >>refercon ApJ, 461, 20 */
33 static const int NHM96 = 14;
34 /* log energy in Ryd */
35 static const double tnuHM96[NHM96]={-8,-1.722735683,-0.351545683,-0.222905683,-0.133385683,
36 /* changeg these two energies to prevent degeneracy */
37 -0.127655683,-0.004575683,0.297544317,0.476753,0.476756,0.588704317,
38 0.661374317,1.500814317,2.245164317};
39 /*-0.127655683,-0.004575683,0.297544317,0.476754317,0.476754317,0.588704317,*/
40 /*log J in the units of (erg cm^{-2} s^{-1} Hz^{-1} sr^{-1})*/
41 static const double fnuHM96[NHM96]={-32.53342863,-19.9789,-20.4204,-20.4443,-20.5756,-20.7546,
42 -21.2796,-21.6256,-21.8404,-21.4823,-22.2102,-22.9263,-23.32,-24.2865};
43 
44 /* Mathews and Ferland generic AGN continuum */
45 static const int NAGN = 8;
46 static Energy tnuagn[NAGN];
48 
49 /* table Draine ISM continuum */
50 static const int NDRAINE = 15;
51 static double tnudrn[NDRAINE] , tsldrn[NDRAINE];
52 
53 /* routine that stores values for above vectors */
54 STATIC void ZeroContin(void);
55 
56 /* this allows the low energy point of any built in array to be reset to the
57  * current low energy point in the code - nothing need be done if this is reset
58  * tnu is array of energies, [0] is first, and we want it to be lower
59  * fluxlog is flux at tnu, and may or may not be log
60  * lgLog says whether it is */
61 STATIC void resetBltin( Energy *tnu , realnum *fluxlog , bool lgLog )
62 {
63  /* this will multiply low-energy bounds of code and go into element[0]
64  * ensures that energy range is fully covered */
65  const double RESETFACTOR = 0.98;
66  double power;
67  /* this makes sure we are called after emm is defined */
68  ASSERT( rfield.emm() > 0. );
69 
70  if( lgLog )
71  {
72  /* continuum comes in as log of flux */
73  /* this is current power-law slope of low-energy continuum */
74  power = (fluxlog[1] - fluxlog[0] ) / log10( tnu[1].Ryd()/tnu[0].Ryd() );
75  /* this will be new low energy bounds to this continuum */
76  tnu[0] = rfield.emm()*RESETFACTOR;
77  fluxlog[0] = fluxlog[1] + power * log10( tnu[0].Ryd()/tnu[1].Ryd() );
78  }
79  else
80  {
81  /* continuum comes in as linear flux */
82  /* this is current power-law slope of low-energy continuum */
83  power = log10( fluxlog[1]/fluxlog[0]) / log10( tnu[1].Ryd()/tnu[0].Ryd() );
84  /* this will be new low energy bounds to this continuum */
85  tnu[0] = rfield.emm()*RESETFACTOR;
86  fluxlog[0] = log10(fluxlog[1]) + power * log10( tnu[0].Ryd()/tnu[1].Ryd() );
87  /* flux is not really log, we want linear */
88  fluxlog[0] = exp10(fluxlog[0]);
89  }
90  /*fprintf(ioQQQ," power is %f lgLog is %i\n", power, lgLog );*/
91  return;
92 }
93 
95 {
96  string chFile; /*file name for table read */
97 
98  IntMode imode = IM_ILLEGAL_MODE;
99  bool lgHit,
100  lgLogSet;
101  static bool lgCalled=false;
102 
103  long int i,
104  j,
105  nstar;
106 
107  double alpha,
108  brakmm,
109  brakxr,
110  ConBreak,
111  fac,
112  scale,
113  slopir,
114  slopxr;
115 
116  bool lgNoContinuum = false,
117  lgQuoteFound;
118 
119  DEBUG_ENTRY( "ParseTable()" );
120 
121  /* if first call then set up values for table */
122  if( !lgCalled )
123  {
124  ZeroContin();
125  lgCalled = true;
126  }
127 
128  if( rfield.nShape >= LIMSPC )
129  {
130  fprintf( ioQQQ, " %ld is too many spectra entered. Increase LIMSPC\n Sorry.\n",
131  rfield.nShape );
133  }
134 
135  /* four commands, tables line, read, SED, and star, have quotes on the
136  * lines giving file names. must get quotes first so that filename
137  * does not confuse parser */
138  lgQuoteFound = true;
139  if( p.GetQuote( chFile ) )
140  lgQuoteFound = false;
141 
142  // Backwards compatibility with older versions of build in SEDs,
143  // no SED and no file name in quotes, just a keyword
144  bool lgKeyword = false;
145  if(!lgQuoteFound)
146  {
147  if(p.nMatch("AKN1"))
148  {
149  chFile = "akn120.sed";
150  lgKeyword = true;
151  lgQuoteFound = true;
152  }
153  else if(p.nMatch("CRAB"))
154  {
155  if(p.nMatch("DAVIDSON"))
156  {
157  chFile = "CrabDavidson.sed";
158  }
159  else
160  {
161  chFile = "CrabHester.sed";
162  }
163  lgKeyword = true;
164  lgQuoteFound = true;
165  }
166  else if(p.nMatch("COOL"))
167  {
168  chFile = "cool.sed";
169  lgKeyword = true;
170  lgQuoteFound = true;
171  }
172  else if(p.nMatch("RUBI"))
173  {
174  chFile = "Rubin.sed";
175  lgKeyword = true;
176  lgQuoteFound = true;
177  }
178  else if(p.nMatch("TRAP"))
179  {
180  chFile = "Trapezium.sed";
181  lgKeyword = true;
182  lgQuoteFound = true;
183  }
184  else if(p.nMatch(" XDR"))
185  {
186  chFile = "XDR.sed";
187  lgKeyword = true;
188  lgQuoteFound = true;
189  }
190  }
191 
192  // we reserve space for 1000 table points, no worries if that is too small though...
193  ASSERT( rfield.tNu[rfield.nShape].empty() );
194  rfield.tNu[rfield.nShape].reserve( 1000 );
195  ASSERT( rfield.tslop[rfield.nShape].empty() );
196  rfield.tslop[rfield.nShape].reserve( 1000 );
197  ASSERT( rfield.tFluxLog[rfield.nShape].empty() );
198  rfield.tFluxLog[rfield.nShape].reserve( 1000 );
199 
200  /* set flag telling interpolate */
201  strcpy( rfield.chSpType[rfield.nShape], "INTER" );
202 
203  bool lgHM05 = false, lgHM12 = false;
204 
205  /* NB when adding more keys also change the comment at the end */
206  if( p.nMatch(" AGN") )
207  {
208  /* do Mathews and Ferland generic AGN continuum */
209  for( i=0; i < NAGN; i++ )
210  {
211  rfield.tNu[rfield.nShape].push_back( tnuagn[i] );
212  rfield.tslop[rfield.nShape].push_back( tslagn[i] );
213  }
215 
216  /* optional keyword break, to adjust IR cutoff */
217  if( p.nMatch("BREA") )
218  {
219  ConBreak = p.FFmtRead();
220 
221  if( p.lgEOL() )
222  {
223  /* no break, set to low energy limit of code */
224  if( p.nMatch(" NO ") )
225  {
226  ConBreak = rfield.emm()*1.01f;
227  }
228  else
229  {
230  fprintf( ioQQQ, " There must be a number for the break.\n Sorry.\n" );
232  }
233  }
234 
235  if( ConBreak == 0. )
236  {
237  fprintf( ioQQQ, " The break must be greater than 0.2 Ryd.\n Sorry.\n" );
239  }
240 
241  if( p.nMatch("MICR") )
242  {
243  /* optional keyword, ``microns'', convert to Rydbergs */
244  ConBreak = 0.0912/ConBreak;
245  }
246 
247  if( ConBreak < 0. )
248  {
249  /* option to enter break as LOG10 */
250  ConBreak = exp10(ConBreak);
251  }
252 
253  else if( ConBreak == 0. )
254  {
255  fprintf( ioQQQ, " An energy of 0 is not allowed.\n Sorry.\n" );
257  }
258 
259  if( ConBreak >= rfield.tNu[rfield.nShape][2].Ryd() )
260  {
261  fprintf( ioQQQ, " The energy of the break cannot be greater than%10.2e Ryd.\n Sorry.\n",
262  rfield.tNu[rfield.nShape][2].Ryd() );
264  }
265 
266  else if( ConBreak <= rfield.tNu[rfield.nShape][0].Ryd() )
267  {
268  fprintf( ioQQQ, " The energy of the break cannot be less than%10.2e Ryd.\n Sorry.\n",
269  rfield.tNu[rfield.nShape][0].Ryd() );
271  }
272 
273  rfield.tNu[rfield.nShape][1].set(ConBreak);
274 
275  rfield.tslop[rfield.nShape][1] =
276  (realnum)(rfield.tslop[rfield.nShape][2] +
277  log10(rfield.tNu[rfield.nShape][2].Ryd()/rfield.tNu[rfield.nShape][1].Ryd()));
278 
279  rfield.tslop[rfield.nShape][0] =
280  (realnum)(rfield.tslop[rfield.nShape][1] -
281  2.5*log10(rfield.tNu[rfield.nShape][1].Ryd()/rfield.tNu[rfield.nShape][0].Ryd()));
282  }
283  }
284 
285  else if( p.nMatchErase("HM96") )
286  {
287  /* this is the old Haardt & Madau continuum, one set of points
288  * with only the quasars
289  * this command does not include the CMB - do that separately with the CMB command */
290  /* set flag telling interpolate */
291  strcpy( rfield.chSpType[rfield.nShape], "INTER" );
292 
293  /* z=2 background,
294  * >>refer continuum background Haardt, Francesco, & Madau, Piero, 1996, ApJ, 461, 20 */
295  for( j=0; j < NHM96; j++ )
296  {
297  /* frequency was stored as log of ryd */
298  rfield.tNu[rfield.nShape].push_back( Energy(exp10( tnuHM96[j])) );
299  rfield.tslop[rfield.nShape].push_back( (realnum)fnuHM96[j] );
300  }
302 
303  /* optional scale factor to change default intensity from their value
304  * assumed to be log if negative, and linear otherwise */
305  scale = p.FFmtRead();
306  if( scale > 0. )
307  scale = log10(scale);
308 
309  /* this also sets continuum intensity*/
310  if( p.m_nqh >= LIMSPC )
311  {
312  fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
313  p.m_nqh );
315  }
316 
317  /* check that stack of shape and luminosity specifications
318  * is parallel, stop if not - this happens is background comes
319  * BETWEEN another set of shape and luminosity commands */
320  if( rfield.nShape != p.m_nqh )
321  {
322  fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
323  fprintf( ioQQQ, " Sorry.\n" );
325  }
326 
327  strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
328  strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
329 
330  /* this is an isotropic radiation field */
331  rfield.lgBeamed[p.m_nqh] = false;
333 
334  /* this will be flux density at some frequency on the table. the numbers
335  * are per Hz and sr so must multiply by 4 pi
336  * [2] is not special, could have been any within array*/
337  rfield.range[p.m_nqh][0] = exp10( tnuHM96[2] )*1.0001;
338 
339  /* convert intensity HM96 give to current units of code */
340  rfield.totpow[p.m_nqh] = (fnuHM96[2] + log10(PI4) + scale);
341 
342  ++p.m_nqh;
343  }
344 
345  else if( ( lgHM05 = p.nMatchErase("HM05") ) || ( lgHM12 = p.nMatchErase("HM12") ) )
346  {
347  bool lgQuasar = p.nMatch("QUAS");
348  if( lgHM12 && lgQuasar )
349  {
350  fprintf( ioQQQ, " The QUASAR option is not supported on the TABLE HM12 command.\n" );
352  }
353  // read requested redshift
354  double redshift = p.FFmtRead();
355  if( p.lgEOL() )
356  {
357  p.NoNumb("redshift");
358  }
359  // read optional scale factor, defaults to 1.
360  double scale = p.FFmtRead();
361  if( scale > 0. )
362  scale = log10(scale);
363  strcpy( rfield.chSpType[rfield.nShape], "VOLK " );
364  UNUSED double zlow, zhigh;
365  int version = lgHM05 ? 2005 : 2012;
366  // this does the hard work of interpolating the SEDs...
367  rfield.ncont[rfield.nShape] = HaardtMadauInterpolate( redshift, version, lgQuasar, &zlow, &zhigh );
368 
369  // choose some frequency for normalizing the spectrum, precise value is not important
370  long ip = rfield.ipointC(0.95);
371  rfield.range[p.m_nqh][0] = rfield.tNu[rfield.nShape][ip].Ryd();
372 
373  // the numbers given by Haardt&Madau are per Hz and sr so must multiply by 4 pi
374  rfield.totpow[p.m_nqh] = log10(rfield.tslop[rfield.nShape][ip]) + log10(PI4) + scale;
375 
376  /* this also sets continuum intensity*/
377  if( p.m_nqh >= LIMSPC )
378  {
379  fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n", p.m_nqh );
381  }
382 
383  /* check that stack of shape and luminosity specifications
384  * is parallel, stop if not - this happens is background comes
385  * BETWEEN another set of shape and luminosity commands */
386  if( rfield.nShape != p.m_nqh )
387  {
388  fprintf( ioQQQ, " This command has come between a previous ordered pair of"
389  " continuum shape and luminosity commands.\n Reorder the commands"
390  " to complete each continuum specification before starting another.\n" );
391  fprintf( ioQQQ, " Sorry.\n" );
393  }
394 
395  strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
396  strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
397 
398  /* this is an isotropic radiation field */
399  rfield.lgBeamed[p.m_nqh] = false;
401 
402  ++p.m_nqh;
403  }
404  else if( p.nMatchErase("KS18") )
405  {
406  // read requested redshift
407  double redshift = p.FFmtRead();
408  if( p.lgEOL() )
409  {
410  p.NoNumb("redshift");
411  }
412  // read optional scale factor, defaults to 1.
413  double scale = p.FFmtRead();
414  if( scale > 0. )
415  scale = log10(scale);
416  // read optional Q value
417  int Q = (int)p.FFmtRead();
418  if( p.lgEOL() )
419  Q = 18;
420  if( Q < 14 || Q > 20 )
421  {
422  fprintf( ioQQQ, " Invalid value Q=%d, should be 14 <= Q <= 20.\n", Q );
424  }
425  strcpy( rfield.chSpType[rfield.nShape], "VOLK " );
426  UNUSED double zlow, zhigh;
427  // this does the hard work of interpolating the SEDs...
428  rfield.ncont[rfield.nShape] = KhaireSrianandInterpolate( redshift, Q, &zlow, &zhigh );
429 
430  // choose some frequency for normalizing the spectrum, precise value is not important
431  long ip = rfield.ipointC(0.95);
432  rfield.range[p.m_nqh][0] = rfield.tNu[rfield.nShape][ip].Ryd();
433 
434  // the numbers given by Khaire&Srianand are per Hz and sr so must multiply by 4 pi
435  rfield.totpow[p.m_nqh] = log10(rfield.tslop[rfield.nShape][ip]) + log10(PI4) + scale;
436 
437  /* this also sets continuum intensity*/
438  if( p.m_nqh >= LIMSPC )
439  {
440  fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n", p.m_nqh );
442  }
443 
444  /* check that stack of shape and luminosity specifications
445  * is parallel, stop if not - this happens is background comes
446  * BETWEEN another set of shape and luminosity commands */
447  if( rfield.nShape != p.m_nqh )
448  {
449  fprintf( ioQQQ, " This command has come between a previous ordered pair of"
450  " continuum shape and luminosity commands.\n Reorder the commands"
451  " to complete each continuum specification before starting another.\n" );
452  fprintf( ioQQQ, " Sorry.\n" );
454  }
455 
456  strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
457  strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
458 
459  /* this is an isotropic radiation field */
460  rfield.lgBeamed[p.m_nqh] = false;
462 
463  ++p.m_nqh;
464  }
465  else if( p.nMatch(" ISM") )
466  {
467  /* local ISM radiation field from Black 1987, Interstellar Processes */
468  /* >>chng 04 mar 16, rm CMB from field so that it can be used at
469  * any redshift */
470  rfield.tNu[rfield.nShape].push_back( Energy( exp10(6.)/FR1RYD) );
471  rfield.tslop[rfield.nShape].push_back( (-21.21f - 6.f) );
472  for( i=6; i < NISM; i++ )
473  {
474  /* energies were stored as log Hz and intensity as log nu Fnu, as per John's plot.
475  * convert to Ryd and log Fnu */
476  rfield.tNu[rfield.nShape].push_back( Energy( exp10(tnuism[i])/FR1RYD ) );
477  rfield.tslop[rfield.nShape].push_back( (realnum)(fnuism[i] - tnuism[i]) );
478  }
479  rfield.ncont[rfield.nShape] = NISM -5;
480 
481  /* optional scale factor to change default luminosity
482  * from observed value
483  * want final number to be log
484  * assumed to be log if negative, and linear otherwise unless log option is present */
485  scale = p.FFmtRead();
486  if( scale > 0. && !p.nMatch(" LOG"))
487  scale = log10(scale);
488 
489  /* this also sets continuum intensity*/
490  if( p.m_nqh >= LIMSPC )
491  {
492  fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
493  p.m_nqh );
495  }
496 
497  /* check that stack of shape and luminosity specifications
498  * is parallel, stop if not - this happens is background comes
499  * BETWEEN another set of shape and luminosity commands */
500  if( rfield.nShape != p.m_nqh )
501  {
502  fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
503  fprintf( ioQQQ, " Sorry.\n" );
505  }
506 
507  strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
508  strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
509 
510  /* this is an isotropic radiation field */
511  rfield.lgBeamed[p.m_nqh] = false;
513 
514  /* this will be flux density at 1 Ryd
515  * >>chng 96 dec 18, from 1 Ryd to H mass Rydberg
516  * >>chng 97 jan 10, had HLevNIonRyd but not defined yet */
517  rfield.range[p.m_nqh][0] = HIONPOT;
518 
519  /* interpolated from Black 1987 */
520  rfield.totpow[p.m_nqh] = (-18.517 + scale);
521 
522  ++p.m_nqh;
523 
524  if( optimize.lgVarOn )
525  {
527  strcpy( optimize.chVarFmt[optimize.nparm], "TABLE ISM LOG %f");
528  /* pointer to where to write */
530  /* the scale factor */
531  optimize.vparm[0][optimize.nparm] = (realnum)scale;
532  optimize.vincr[optimize.nparm] = 0.2f;
533  ++optimize.nparm;
534  }
535  }
536  else if( p.nMatch("DRAI") )
537  {
538  rfield.lgMustBlockHIon = true;
539  /* local ISM radiation field from equation 23
540  *>>refer ISM continuum Draine & Bertoldi 1996 */
541  for( i=0; i < NDRAINE; i++ )
542  {
543  rfield.tNu[rfield.nShape].push_back( Energy(tnudrn[i]) );
544  rfield.tslop[rfield.nShape].push_back( (realnum)tsldrn[i] );
545  }
547 
548  /* optional scale factor to change default luminosity
549  * from observed value
550  * assumed to be log if negative, and linear otherwise unless log option is present */
551  scale = p.FFmtRead();
552  if( scale > 0. && !p.nMatch(" LOG") )
553  scale = log10(scale);
554 
555  /* this also sets continuum intensity*/
556  if( p.m_nqh >= LIMSPC )
557  {
558  fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
559  p.m_nqh );
561  }
562 
563  /* check that stack of shape and luminosity specifications
564  * is parallel, stop if not - this happens is background comes
565  * BETWEEN another set of shape and luminosity commands */
566  if( rfield.nShape != p.m_nqh )
567  {
568  fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
569  fprintf( ioQQQ, " Sorry.\n" );
571  }
572 
573  strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
574  strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
575 
576  /* this is an isotropic radiation field */
577  rfield.lgBeamed[p.m_nqh] = false;
579 
580  /* continuum normalization given by flux density at first point,
581  * must set energy a bit higher to make sure it is within energy bounds
582  * that results from float arithmetic */
583  rfield.range[p.m_nqh][0] = tnudrn[0]*1.01;
584 
585  /* this is f_nu at this first point */
586  rfield.totpow[p.m_nqh] = tsldrn[0] + scale;
587 
588  if( optimize.lgVarOn )
589  {
591  strcpy( optimize.chVarFmt[optimize.nparm], "TABLE DRAINE LOG %f");
592  /* pointer to where to write */
594  /* the scale factor */
595  optimize.vparm[0][optimize.nparm] = (realnum)scale;
596  optimize.vincr[optimize.nparm] = 0.2f;
597  ++optimize.nparm;
598  }
599 
600  ++p.m_nqh;
601  }
602 
603  // match LINES to avoid confusion with LINEAR
604  else if( p.nMatch("LINES") )
605  {
606  /* table lines command - way to check that lines within a data
607  * file are still valid */
608 
609  /* say that this is not a continuum command, so don't try to work with unallocated space */
610  /* this is not a continuum source - it is to read a table of lines */
611  lgNoContinuum = true;
612 
613  if( chLINE_LIST.size() > 0 )
614  {
615  fprintf(ioQQQ," sorry, only one table line per input stream\n");
617  }
618 
619  /* get file name within double quotes, if not present will use default
620  * return value of 1 indicates did not find double quotes on line */
621  if( lgQuoteFound && chFile.length() > 0 )
622  chLINE_LIST = chFile;
623  else
624  chLINE_LIST = "LineList_BLR.dat";
625 
626  // check if the file exists
627  FILE* ioData = open_data( chLINE_LIST.c_str(), "r", AS_LOCAL_DATA_TRY );
628  if( ioData == NULL )
629  {
630  /* did not find file, abort */
631  fprintf(ioQQQ,"\n DISASTER PROBLEM ParseTable could not find "
632  "line list file %s\n", chLINE_LIST.c_str() );
633  fprintf(ioQQQ," Please check the spelling of the file name and that it "
634  "is in either the local or data directory.\n\n");
636  }
637  else
638  {
639  fclose(ioData);
640  }
641  /* actually reading the data is done in lines_table() */
642  }
643 
644  else if( p.nMatch("POWE") )
645  {
646  /* simple power law continuum between 10 micron and 50 keV
647  * option to read in any slope for the intermediate continuum */
648  alpha = p.FFmtRead();
649 
650  /* default (no number on line) is f_nu proportional nu^-1 */
651  if( p.lgEOL() )
652  alpha = -1.;
653 
654  /* this is low energy for code */
655  rfield.tNu[rfield.nShape].push_back( Energy(rfield.emm()) );
656  /* and the value of the flux at this point (f_nu units)*/
657  rfield.tslop[rfield.nShape].push_back( -5.f );
658 
659  lgLogSet = false;
660 
661  /* option to adjust sub-millimeter break */
662  brakmm = p.FFmtRead();
663 
664  /* default is 10 microns */
665  if( p.lgEOL() )
666  {
667  lgLogSet = true;
668  brakmm = 9.115e-3;
669  }
670 
671  else if( brakmm == 0. )
672  {
673  /* if second number on line is zero then set lower limit to
674  * low-energy limit of the code. Also set linear mode,
675  * so that last number will also be linear. */
676  lgLogSet = false;
677  brakmm = rfield.tNu[rfield.nShape][0].Ryd()*1.001;
678  }
679 
680  else if( brakmm < 0. )
681  {
682  /* if number is negative then this and next are logs */
683  lgLogSet = true;
684  brakmm = exp10(brakmm);
685  }
686 
687  /* optional microns keyword - convert to Rydbergs */
688  if( p.nMatch("MICR") )
689  brakmm = RYDLAM / (1e4*brakmm);
690 
691  rfield.tNu[rfield.nShape].push_back( Energy(brakmm) );
692 
693  /* check whether this is a reasonable mm break */
694  if( brakmm > 1. )
695  fprintf(ioQQQ,
696  " Check the order of parameters on this table power law command - the low-energy break of %f Ryd seems high to me.\n",
697  brakmm );
698 
699  /* this is spectral slope, in F_nu units, between the low energy limit
700  * and the break that may have been adjusted above
701  * this is the slope appropriate for self-absorbed synchrotron, see eq 6.54, p.190
702  *>>refer continuum synchrotron Rybicki, G. B., & Lightman, A.P. 1979,
703  *>>refercon Radiative Processes in Astrophysics (New York: Wiley)*/
704  slopir = 5./2.;
705 
706  /* now extrapolate a flux at this energy using the flux entered for
707  * the first point, and this slope */
708  rfield.tslop[rfield.nShape].push_back(
709  (realnum)(rfield.tslop[rfield.nShape][0] +
710  slopir*log10(rfield.tNu[rfield.nShape][1].Ryd()/rfield.tNu[rfield.nShape][0].Ryd()))
711  );
712 
713  /* option to adjust hard X-ray break */
714  brakxr = p.FFmtRead();
715 
716  /* default is 50 keV */
717  if( p.lgEOL() )
718  {
719  brakxr = 3676.;
720  }
721 
722  else if( brakxr == 0. )
723  {
724  brakxr = rfield.egamry()/1.001;
725  }
726 
727  else if( lgLogSet )
728  {
729  /* first number was negative this is a logs */
730  brakxr = exp10(brakxr);
731  }
732 
733  /* note that this second cutoff does not have the micron keyword */
734  rfield.tNu[rfield.nShape].push_back( Energy(brakxr) );
735 
736  /* this is energy of the high-energy limit to code */
737  rfield.tNu[rfield.nShape].push_back( Energy(rfield.egamry()) );
738 
739  /* >>chng 03 jul 19, check that upper energy is greater than lower energy,
740  * quit if this is not the case */
741  if( brakmm >= brakxr )
742  {
743  fprintf( ioQQQ, " HELP!! The lower energy for the power law, %f, is greater than the upper energy, %f. This is not possible.\n Sorry.\n",
744  brakmm , brakxr );
746  }
747 
748  /* alpha was first option on line, is slope of mid-range */
749  rfield.tslop[rfield.nShape].push_back(
750  (realnum)(rfield.tslop[rfield.nShape][1] +
751  alpha*log10(rfield.tNu[rfield.nShape][2].Ryd()/rfield.tNu[rfield.nShape][1].Ryd()))
752  );
753 
754  /* high energy range is nu^-2 */
755  slopxr = -2.;
756 
757  rfield.tslop[rfield.nShape].push_back(
758  (realnum)(rfield.tslop[rfield.nShape][2] +
759  slopxr*log10(rfield.tNu[rfield.nShape][3].Ryd()/rfield.tNu[rfield.nShape][2].Ryd()))
760  );
761 
762  /* following is number of portions of continuum */
763  rfield.ncont[rfield.nShape] = 4;
764  }
765 
766  else if( p.nMatch("READ") )
767  {
768  /* set up eventual read of table of points previously punched by code
769  * get file name within double quotes, return as null terminated string
770  * also blank out original, chCard version of name so do not trigger */
771  if( !lgQuoteFound )
772  {
773  fprintf( ioQQQ, " Name of file must appear on TABLE READ.\n");
775  }
776 
777  /* set flag saying really just read in continuum exactly as punched */
778  strcpy( rfield.chSpType[rfield.nShape], "READ " );
779 
780  ReadTable(chFile.c_str());
781 
782  if ( p.nMatch("SCALE") )
783  {
784  /* optional scale factor to change default intensity from their value
785  * assumed to be log if negative, and linear otherwise
786  * increment i to move past the 96 in the keyword */
787  scale = p.FFmtRead();
788  if( scale > 0. && !p.nMatch(" LOG"))
789  scale = log10(scale);
790  /* this also sets continuum intensity*/
791  if( p.m_nqh >= LIMSPC )
792  {
793  fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
794  p.m_nqh );
796  }
797  /* check that stack of shape and luminosity specifications
798  * is parallel, stop if not - this happens is background comes
799  * BETWEEN another set of shape and luminosity commands */
800  if( rfield.nShape != p.m_nqh )
801  {
802  fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
803  fprintf( ioQQQ, " Sorry.\n" );
805  }
806 
807  strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
808  strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
809  // rfield.lgBeamed[p.m_nqh] = false;
810  // rfield.Illumination[p.m_nqh] = Illuminate::ISOTROPIC;
811 
812  rfield.range[p.m_nqh][0] = rfield.tNu[rfield.nShape][0].Ryd();
813  double fmax = -70.;
814  for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
815  {
816  if (rfield.tFluxLog[rfield.nShape][i] > fmax)
817  {
818  fmax = rfield.tFluxLog[rfield.nShape][i];
819  rfield.range[p.m_nqh][0] = rfield.tNu[rfield.nShape][i].Ryd();
820  }
821  }
822 
824 
825  // EN1RYD/FR1RYD == HPLANCK -- EN1RYD scaling already applied when saved
826  rfield.totpow[p.m_nqh] = scale+fmax-log10(FR1RYD);
827 
828  ++p.m_nqh;
829  }
830  /* number of spectra shapes that have been specified */
831  ++rfield.nShape;
832 
833  return;
834  }
835 
836  else if( p.nMatch("TLUSTY") && !p.nMatch("STAR") )
837  {
838  /* >>chng 04 nov 30, retired TABLE TLUSTY command, PvH */
839  fprintf( ioQQQ, " The TABLE TLUSTY command is no longer supported.\n" );
840  fprintf( ioQQQ, " Please use TABLE STAR TLUSTY instead. See Hazy for details.\n" );
842  }
843 
844  else if( p.nMatch("SED") || lgKeyword )
845  {
846  bool lgEOL;
847  char chLine[INPUT_LINE_LENGTH];
848 
849  if( !lgQuoteFound )
850  {
851  fprintf( ioQQQ, "PROBLEM in TABLE SED: No quotes were found.\n The SED table file must be designated in quotes.\n" );
853  }
854 
855  /* now open the data file -- first try local directory, then data/SED */
856  char chPath[FILENAME_PATH_LENGTH_2];
857  FILE *ioDATA;
858  if( (ioDATA = open_data( chFile.c_str(), "r", AS_LOCAL_ONLY_TRY ) ) == NULL )
859  {
860  strcpy( chPath, "SED" );
861  strcat( chPath, input.chDelimiter );
862  strcat( chPath, chFile.c_str() );
863  ioDATA = open_data( chPath, "r" );
864  }
865 
866  /* read data first time, counting size of required tables,
867  * and parse user optional keywords */
868  long int fileLength = 0;
869  const char *chEnergyUnits=NULL;
870  bool lgUnitsSet = false, lgNuFnu = false, lgExtrapolate = false;
871 
872  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
873  {
874  // field of stars ends data
875  if( chLine[0]=='*' )
876  break;
877 
878  /* skip comment line */
879  if( chLine[0]=='#' )
880  continue;
881 
882  if( chLine[0]=='\n' || chLine[0]=='\0' )
883  {
884  fprintf(ioQQQ, "PROBLEM in TABLE SED: Encountered unexpected empty line.\n");
886  }
887 
888  /* copy input to CAPS to parse keywords */
889  char chLineCAPS[INPUT_LINE_LENGTH];
890  strcpy( chLineCAPS , chLine );
891  caps( chLineCAPS );
892 
893  /* change photon units */
894  if(nMatch( "UNIT", chLineCAPS ) != 0)
895  {
896  chEnergyUnits = StandardEnergyUnit(chLineCAPS);
897  lgUnitsSet = true;
898  }
899 
900  /* default is F_nu - this makes it nu F_nu */
901  if(nMatch( "NUFN" , chLineCAPS ) != 0)
902  lgNuFnu = true;
903 
904  /* default is use exactly the specified continuum, extrapolate option
905  * will extrapolate lowest energy point to low-energy limit of the code */
906  if(nMatch( "EXTRAP" , chLineCAPS ) != 0)
907  lgExtrapolate = true;
908 
909  fileLength++;
910  }
911 
912  if( fileLength < 2 )
913  {
914  fprintf(ioQQQ, "PROBLEM in TABLE SED: less than two data pairs found.\n");
916  }
917 
918  vector<double> tnuInput(fileLength);
919  vector<double> tslopInput(fileLength);
920  rewind(ioDATA);
921 
922  /* read in and save data */
923  long int entryNum = 0;
924  Energy unitChange;
925  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
926  {
927  if( chLine[0]=='*' )
928  break;
929  /* skip comment */
930  if( chLine[0]=='#' )
931  continue;
932 
933  long i = 1;
934  tnuInput[entryNum] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
935  if( lgEOL)
936  {
937  fprintf(ioQQQ,"PROBLEM in TABLE SED: freq/wavl and flux pairs must be entered."
938  " The first number was not found on this line:\n%s\n", chLine);
940  }
941  tslopInput[entryNum] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
942  /* these were not checked during the first read through the data since
943  * no data were actually parsed
944  */
945  if( lgEOL)
946  {
947  fprintf(ioQQQ,"PROBLEM in TABLE SED: freq/wavl and flux pairs must be entered."
948  " The second number was not found on this line:\n%s\n", chLine);
950  }
951 
952  if( tslopInput[entryNum] <= 0. )
953  {
954  fprintf(ioQQQ,"PROBLEM flux in TABLE SED file entry %li is not positive\nphoton energy %e where the value is %e.\n",
955  entryNum , tnuInput[entryNum] , tslopInput[entryNum] );
956  fprintf(ioQQQ,"Line image follows:\n%s\n", chLine );
958  }
959 
960  if( tnuInput[entryNum] <= 0. )
961  {
962  fprintf(ioQQQ,"PROBLEM freq/wavl in TABLE SED file entry %li is not positive\n"
963  "freq/wavl is %e.\n",
964  entryNum , tnuInput[entryNum] );
965  fprintf(ioQQQ,"Line image follows:\n%s\n", chLine );
967  }
968 
969  if( entryNum > 1 && (tnuInput[entryNum-2]-tnuInput[entryNum-1])*
970  (tnuInput[entryNum-1]-tnuInput[entryNum]) <= 0. )
971  {
972  fprintf(ioQQQ,"PROBLEM freq/wavl values in TABLE SED file must be in increasing or decreasing order\n");
973  fprintf(ioQQQ, "The following freq/wavl sequence is not monotonically changing:\n");
974  for( long ij=entryNum-2; ij<=entryNum; ++ij )
975  fprintf(ioQQQ, "%li %e\n", ij, tnuInput[ij] );
977  }
978 
979  tslopInput[entryNum] = log10(tslopInput[entryNum]);
980  entryNum++;
981  }
982 
983  ASSERT( entryNum == fileLength );
984 
985  //Unit Conversion of frequency, we want linear Rydbergs
986  if( lgUnitsSet )
987  {
988  for( long i=0; i < entryNum; i++ )
989  {
990  unitChange.set(tnuInput[i], chEnergyUnits);
991  tnuInput[i] = unitChange.Ryd();
992  }
993  }
994 
995  /* we want to specify log F_nu - this converts linear nh F+nu to log F_nu */
996  if( lgNuFnu )
997  {
998  for( long i=0; i < entryNum; i++ )
999  tslopInput[i] -= log10(tnuInput[i]);
1000  }
1001 
1002  if( tnuInput[0] < tnuInput[1] )
1003  {
1004  for( long i=0; i < entryNum; i++ )
1005  {
1006  rfield.tNu[rfield.nShape].push_back( Energy(tnuInput[i]) );
1007  rfield.tslop[rfield.nShape].push_back( (realnum)tslopInput[i] );
1008  }
1009  }
1010  else
1011  {
1012  for( long i=0; i < entryNum; i++ )
1013  {
1014  rfield.tNu[rfield.nShape].push_back( Energy(tnuInput[entryNum-i-1]) );
1015  rfield.tslop[rfield.nShape].push_back( (realnum)tslopInput[entryNum-i-1] );
1016  }
1017  }
1018  rfield.ncont[rfield.nShape] = entryNum;
1019 
1020  /* option to extrapolate lowest energy point specified to low-energy limit of code */
1021  if( lgExtrapolate )
1023 
1024  fclose( ioDATA );
1025  }
1026 
1027  /* >>chng 06 jul 10, retired TABLE STARBURST command, PvH */
1028 
1029  else if( p.nMatch("STAR") )
1030  {
1031  char chMetalicity[5] = {'\0'}, chODFNew[8], chVaryFlag[7] = {'\0'};
1032  bool lgHCa = false, lgHNi = false;
1033  long nval, ndim=0;
1034  double Tlow = -1., Thigh = -1.;
1035  double val[MDIM];
1036 
1037  /* >>chng 06 jun 22, add support for 3 and 4-dimensional grids, PvH */
1038  if( p.nMatchErase("1-DI") )
1039  ndim = 1;
1040  else if( p.nMatchErase("2-DI") )
1041  ndim = 2;
1042  else if( p.nMatchErase("3-DI") )
1043  ndim = 3;
1044  else if( p.nMatchErase("4-DI") )
1045  ndim = 4;
1046 
1047  if( ndim != 0 )
1048  {
1049  /* remember keyword for possible use in optimization command */
1050  sprintf(chVaryFlag,"%1ld-DIM",ndim);
1051  }
1052 
1053  /* time option to vary only some continua with time */
1054  rfield.lgTimeVary[p.m_nqh] = p.nMatch( "TIME" );
1055 
1056  static const char table[][2][10] = {
1057  {"Z+1.0 ", "p10"},
1058  {"Z+0.75", "p075"},
1059  {"Z+0.5 ", "p05"},
1060  {"Z+0.4 ", "p04"},
1061  {"Z+0.3 ", "p03"},
1062  {"Z+0.25", "p025"},
1063  {"Z+0.2 ", "p02"},
1064  {"Z+0.1 ", "p01"},
1065  {"Z+0.0 ", "p00"},
1066  {"Z-0.1 ", "m01"},
1067  {"Z-0.2 ", "m02"},
1068  {"Z-0.25", "m025"},
1069  {"Z-0.3 ", "m03"},
1070  {"Z-0.4 ", "m04"},
1071  {"Z-0.5 ", "m05"},
1072  {"Z-0.7 ", "m07"},
1073  {"Z-0.75", "m075"},
1074  {"Z-1.0 ", "m10"},
1075  {"Z-1.3 ", "m13"},
1076  {"Z-1.5 ", "m15"},
1077  {"Z-1.7 ", "m17"},
1078  {"Z-2.0 ", "m20"},
1079  {"Z-2.5 ", "m25"},
1080  {"Z-3.0 ", "m30"},
1081  {"Z-3.5 ", "m35"},
1082  {"Z-4.0 ", "m40"},
1083  {"Z-4.5 ", "m45"},
1084  {"Z-5.0 ", "m50"},
1085  {"Z-INF ", "m99"}
1086  };
1087 
1088  strncpy( chMetalicity, "p00", sizeof(chMetalicity) ); // default
1089  for (i=0; i < (long)(sizeof(table)/sizeof(*table)); ++i)
1090  {
1091  if (p.nMatchErase(table[i][0]))
1092  {
1093  strncpy( chVaryFlag, table[i][0], sizeof(chVaryFlag)-1 );
1094  strncpy( chMetalicity, table[i][1], sizeof(chMetalicity)-1 );
1095  break;
1096  }
1097  }
1098 
1099 
1100  /* there are two abundance sets (solar and halo) for CoStar and Rauch H-Ca/H-Ni models.
1101  * If halo keyword appears use halo, else use solar unless other metalicity was requested */
1102  bool lgHalo = p.nMatch("HALO");
1103  bool lgSolar = ( !lgHalo && strcmp( chMetalicity, "p00" ) == 0 );
1104 
1105  /* >>chng 05 oct 27, added support for PG1159 Rauch models, PvH */
1106  bool lgPG1159 = p.nMatchErase("PG1159");
1107 
1108  bool lgList = p.nMatch("LIST");
1109 
1110  if( p.nMatch("AVAI") )
1111  {
1112  AtmospheresAvail();
1114  }
1115 
1116  /* now that all the keywords are out of the way, scan numbers from line image */
1117  for( nval=0; nval < MDIM; nval++ )
1118  {
1119  val[nval] = p.FFmtRead();
1120  if( p.lgEOL() )
1121  break;
1122  }
1123 
1124  if( nval == 0 && !lgList )
1125  {
1126  fprintf( ioQQQ, " The stellar temperature MUST be entered.\n" );
1128  }
1129 
1130  /* match on "LOG " rather than " LOG" to avoid confusion with "LOG(G)" */
1131  if( p.nMatch("LOG ") )
1132  {
1133  if( val[0] < log10(BIGDOUBLE) )
1134  val[0] = exp10(val[0]);
1135  else
1136  {
1137  fprintf(ioQQQ," DISASTER the log of the parameter was specified, "
1138  "but the numerical value is too large.\n Sorry.\n\n");
1139  cdEXIT(EXIT_FAILURE );
1140  }
1141  }
1142 
1143  if( lgQuoteFound )
1144  {
1145  nstar = GridInterpolate(val,&nval,&ndim,chFile.c_str(),lgList,&Tlow,&Thigh);
1146  }
1147 
1148  else if( p.nMatch("ATLA") )
1149  {
1150  /* this sub-branch added by Kevin Volk, to read in large
1151  * grid of Kurucz atmospheres */
1152  /* >>chng 05 nov 19, option TRACE is no longer supported, PvH */
1153 
1154  if( p.nMatch("ODFN") )
1155  strncpy( chODFNew, "_odfnew", sizeof(chODFNew) );
1156  else
1157  strncpy( chODFNew, "", sizeof(chODFNew) );
1158 
1159  /* >>chng 05 nov 19, add support for non-solar metalicities, PvH */
1160  nstar = AtlasInterpolate(val,&nval,&ndim,chMetalicity,chODFNew,lgList,&Tlow,&Thigh);
1161  }
1162 
1163  else if( p.nMatch("COST") )
1164  {
1165  /* >>chng 99 apr 30, added CoStar stellar atmospheres */
1166  /* this subbranch reads in CoStar atmospheres, no log(g),
1167  * but second parameter is age sequence, a number between 1 and 7,
1168  * default is 1 */
1169  if( p.nMatch("INDE") )
1170  {
1171  imode = IM_COSTAR_TEFF_MODID;
1172  if( nval == 1 )
1173  {
1174  val[1] = 1.;
1175  nval++;
1176  }
1177 
1178  /* now make sure that second parameter is within allowed range -
1179  * we do not have enough information at this time to verify temperature */
1180  if( val[1] < 1. || val[1] > 7. )
1181  {
1182  fprintf( ioQQQ, " The second number must be the id sequence number, 1 to 7.\n" );
1183  fprintf( ioQQQ, " reset to 1.\n" );
1184  val[1] = 1.;
1185  }
1186  }
1187  else if( p.nMatch("ZAMS") )
1188  {
1189  imode = IM_COSTAR_MZAMS_AGE;
1190  if( nval == 1 )
1191  {
1192  fprintf( ioQQQ, " A second number, the age of the star, must be present.\n" );
1194  }
1195  }
1196  else if( p.nMatch(" AGE") )
1197  {
1198  imode = IM_COSTAR_AGE_MZAMS;
1199  if( nval == 1 )
1200  {
1201  fprintf( ioQQQ, " A second number, the ZAMS mass of the star, must be present.\n" );
1203  }
1204  }
1205  else
1206  {
1207  if( nval == 1 )
1208  {
1209  imode = IM_COSTAR_TEFF_MODID;
1210  /* default is to use ZAMS models, i.e. use index 1 */
1211  val[1] = 1.;
1212  nval++;
1213  }
1214  else
1215  {
1216  /* Caution: the code in CoStarInterpolate implicitly
1217  * assumes that the dependence of log(g) with age is
1218  * strictly monotonic. As of June 1999 this is the case. */
1219  imode = IM_COSTAR_TEFF_LOGG;
1220  }
1221  }
1222 
1223  if( ! ( lgSolar || lgHalo ) )
1224  {
1225  fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" );
1227  }
1228 
1229  nstar = CoStarInterpolate(val,&nval,&ndim,imode,lgHalo,lgList,&Tlow,&Thigh);
1230  }
1231 
1232  else if( p.nMatch("KURU") )
1233  {
1234  nstar = Kurucz79Interpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1235  }
1236 
1237  else if( p.nMatch("MIHA") )
1238  {
1239  nstar = MihalasInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1240  }
1241 
1242  else if( p.nMatch("RAUC") )
1243  {
1244  if( ! ( lgSolar || lgHalo ) )
1245  {
1246  fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" );
1248  }
1249 
1250  /* >>chng 97 aug 23, K Volk added Rauch stellar atmospheres */
1251  /* >>chng 05 oct 27, added support for PG1159 Rauch models, PvH */
1252  /* >>chng 06 jun 26, added support for H, He, H+He Rauch models, PvH */
1253  if( p.nMatch("H-CA") || p.nMatch(" OLD") )
1254  {
1255  lgHCa = true;
1256  nstar = RauchInterpolateHCa(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh);
1257  }
1258  else if( p.nMatch("HYDR") )
1259  {
1260  nstar = RauchInterpolateHydr(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1261  }
1262  else if( p.nMatch("HELI") )
1263  {
1264  nstar = RauchInterpolateHelium(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1265  }
1266  else if( p.nMatch("H+HE") )
1267  {
1268  nstar = RauchInterpolateHpHe(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1269  }
1270  else if( lgPG1159 ) /* the keyword PG1159 was already matched and erased above */
1271  {
1272  nstar = RauchInterpolatePG1159(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1273  }
1274  else if( p.nMatch("CO W") )
1275  {
1276  nstar = RauchInterpolateCOWD(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1277  }
1278  else
1279  {
1280  lgHNi = true;
1281  nstar = RauchInterpolateHNi(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh);
1282  }
1283  }
1284 
1285  else if( p.nMatch("TLUS") )
1286  {
1287  if( p.nMatch("OBST") )
1288  {
1289  /* >>chng 09 dec 24, this sub-branch added to read in
1290  * merged Tlusty O-star & B-star atmospheres, PvH */
1291  nstar = TlustyInterpolate(val,&nval,&ndim,TL_OBSTAR,chMetalicity,lgList,&Tlow,&Thigh);
1292  }
1293  else if( p.nMatch("BSTA") )
1294  {
1295  /* >>chng 06 oct 19, this sub-branch added to read in
1296  * large 2006 grid of Tlusty B-star atmospheres, PvH */
1297  nstar = TlustyInterpolate(val,&nval,&ndim,TL_BSTAR,chMetalicity,lgList,&Tlow,&Thigh);
1298  }
1299  else if( p.nMatch("OSTA") )
1300  {
1301  /* >>chng 05 nov 21, this sub-branch added to read in
1302  * large 2002 grid of Tlusty O-star atmospheres, PvH */
1303  nstar = TlustyInterpolate(val,&nval,&ndim,TL_OSTAR,chMetalicity,lgList,&Tlow,&Thigh);
1304  }
1305  else
1306  {
1307  fprintf( ioQQQ, " There must be a third key on TABLE STAR TLUSTY;" );
1308  fprintf( ioQQQ, " the options are OBSTAR, BSTAR, OSTAR.\n" );
1310  }
1311  }
1312 
1313  else if( p.nMatch("WERN") )
1314  {
1315  /* this subbranch added by Kevin Volk, to read in large
1316  * grid of hot PN atmospheres computed by Klaus Werner */
1317  nstar = WernerInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1318  }
1319 
1320  else if( p.nMatch("WMBA") )
1321  {
1322  /* >>chng 06 jun 27, this subbranch added to read in
1323  * grid of hot atmospheres computed by Pauldrach */
1324  nstar = WMBASICInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
1325  }
1326 
1327  else
1328  {
1329  fprintf( ioQQQ, " There must be a second key on TABLE STAR;" );
1330  fprintf( ioQQQ, " the options are ATLAS, KURUCZ, MIHALAS, RAUCH, WERNER, and WMBASIC.\n" );
1332  }
1333 
1334  /* set flag saying really just read in continuum exactly as punched */
1335  strcpy( rfield.chSpType[rfield.nShape], "VOLK " );
1336 
1337  rfield.ncont[rfield.nShape] = nstar;
1338 
1339  /* vary option */
1340  if( optimize.lgVarOn )
1341  {
1343 
1344  if( lgQuoteFound )
1345  {
1346  /* this is number of parameters to feed onto the input line */
1347  optimize.nvarxt[optimize.nparm] = nval;
1348  sprintf( optimize.chVarFmt[optimize.nparm], "TABLE STAR \"%s\"", chFile.c_str() );
1349  strcat( optimize.chVarFmt[optimize.nparm], " %f LOG" );
1350  for( i=1; i < nval; i++ )
1351  strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1352  }
1353 
1354  if( p.nMatch("ATLA") )
1355  {
1356  /* this is number of parameters to feed onto the input line */
1357  optimize.nvarxt[optimize.nparm] = ndim;
1358  strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR ATLAS " );
1359  strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag );
1360  if( p.nMatch("ODFN") )
1361  strcat( optimize.chVarFmt[optimize.nparm], " ODFNEW" );
1362 
1363  strcat( optimize.chVarFmt[optimize.nparm], " %f LOG %f" );
1364 
1365  if( ndim == 3 )
1366  strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1367 
1368  }
1369 
1370  else if( p.nMatch("COST") )
1371  {
1372  /* this is number of parameters to feed onto the input line */
1374 
1375  strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR COSTAR " );
1376  if( lgHalo )
1377  strcat( optimize.chVarFmt[optimize.nparm], "HALO " );
1378 
1379  if( imode == IM_COSTAR_TEFF_MODID )
1380  {
1381  strcat( optimize.chVarFmt[optimize.nparm], "%f LOG , INDEX %f" );
1382  }
1383  else if( imode == IM_COSTAR_TEFF_LOGG )
1384  {
1385  strcat( optimize.chVarFmt[optimize.nparm], "%f LOG %f" );
1386  }
1387  else if( imode == IM_COSTAR_MZAMS_AGE )
1388  {
1389  strcat( optimize.chVarFmt[optimize.nparm], "ZAMS %f LOG %f" );
1390  }
1391  else if( imode == IM_COSTAR_AGE_MZAMS )
1392  {
1393  strcat( optimize.chVarFmt[optimize.nparm], "AGE %f LOG %f" );
1395  }
1396  }
1397 
1398  else if( p.nMatch("KURU") )
1399  {
1400  /* this is number of parameters to feed onto the input line */
1402  strcpy( optimize.chVarFmt[optimize.nparm],
1403  "TABLE STAR KURUCZ %f LOG" );
1404  }
1405 
1406  else if( p.nMatch("MIHA") )
1407  {
1408  /* this is number of parameters to feed onto the input line */
1410  strcpy( optimize.chVarFmt[optimize.nparm],
1411  "TABLE STAR MIHALAS %f LOG" );
1412  }
1413 
1414  else if( p.nMatch("RAUC") )
1415  {
1416  /* this is number of parameters to feed onto the input line */
1417  optimize.nvarxt[optimize.nparm] = ndim;
1418 
1419  strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR RAUCH " );
1420 
1421  if( p.nMatch("HYDR") )
1422  strcat( optimize.chVarFmt[optimize.nparm], "HYDR " );
1423  else if( p.nMatch("HELI") )
1424  strcat( optimize.chVarFmt[optimize.nparm], "HELIUM " );
1425  else if( p.nMatch("H+HE") )
1426  strcat( optimize.chVarFmt[optimize.nparm], "H+HE " );
1427  else if( lgPG1159 )
1428  strcat( optimize.chVarFmt[optimize.nparm], "PG1159 " );
1429  else if( p.nMatch("CO W") )
1430  strcat( optimize.chVarFmt[optimize.nparm], "CO WD " );
1431  else if( lgHCa )
1432  strcat( optimize.chVarFmt[optimize.nparm], "H-CA " );
1433 
1434  if( ( lgHCa || lgHNi ) && ndim == 2 )
1435  {
1436  if( lgHalo )
1437  strcat( optimize.chVarFmt[optimize.nparm], "HALO " );
1438  else
1439  strcat( optimize.chVarFmt[optimize.nparm], "SOLAR " );
1440  }
1441 
1442  strcat( optimize.chVarFmt[optimize.nparm], "%f LOG %f" );
1443 
1444  if( ndim == 3 )
1445  {
1446  if( p.nMatch("H+HE") )
1447  strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1448  else
1449  strcat( optimize.chVarFmt[optimize.nparm], " %f 3-DIM" );
1450  }
1451  }
1452 
1453  if( p.nMatch("TLUS") )
1454  {
1455  /* this is number of parameters to feed onto the input line */
1456  optimize.nvarxt[optimize.nparm] = ndim;
1457  strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR TLUSTY " );
1458  if( p.nMatch("OBST") )
1459  strcat( optimize.chVarFmt[optimize.nparm], "OBSTAR " );
1460  else if( p.nMatch("BSTA") )
1461  strcat( optimize.chVarFmt[optimize.nparm], "BSTAR " );
1462  else if( p.nMatch("OSTA") )
1463  strcat( optimize.chVarFmt[optimize.nparm], "OSTAR " );
1464  else
1465  TotalInsanity();
1466  strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag );
1467  strcat( optimize.chVarFmt[optimize.nparm], " %f LOG %f" );
1468 
1469  if( ndim == 3 )
1470  strcat( optimize.chVarFmt[optimize.nparm], " %f" );
1471  }
1472 
1473  else if( p.nMatch("WERN") )
1474  {
1475  /* this is number of parameters to feed onto the input line */
1477  strcpy( optimize.chVarFmt[optimize.nparm],
1478  "TABLE STAR WERNER %f LOG %f" );
1479  }
1480 
1481  else if( p.nMatch("WMBA") )
1482  {
1483  /* this is number of parameters to feed onto the input line */
1485  strcpy( optimize.chVarFmt[optimize.nparm],
1486  "TABLE STAR WMBASIC %f LOG %f %f" );
1487  }
1488 
1489  if( rfield.lgTimeVary[p.m_nqh] )
1490  strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
1491 
1492  /* pointer to where to write */
1494 
1495  ASSERT( nval <= LIMEXT );
1496 
1497  /* log of temp will be pointer */
1498  optimize.vparm[0][optimize.nparm] = (realnum)log10(val[0]);
1499  for( i=1; i < nval; i++ )
1500  optimize.vparm[i][optimize.nparm] = (realnum)val[i];
1501 
1502  /* following are upper and lower limits to temperature range */
1503  optimize.varang[optimize.nparm][0] = (realnum)log10(Tlow);
1504  optimize.varang[optimize.nparm][1] = (realnum)log10(Thigh);
1505 
1506  /* finally increment this counter */
1507  ++optimize.nparm;
1508  }
1509  }
1510 
1511  else
1512  {
1513  fprintf( ioQQQ, " There MUST be a keyword on this line. The keys are:"
1514  "AGN, DRAINE, HM96, HM05, HM12, _ISM, LINES, POWERlaw, "
1515  "READ, SED, and STAR.\n Sorry.\n" );
1517  }
1518 
1519  /* table star and table hm05 / hm12 are special cases
1520  * they are not really tables at all,
1521  * so return if chSpType is "VOLK " */
1522  if( strcmp(rfield.chSpType[rfield.nShape],"VOLK ") == 0 )
1523  {
1524  ++rfield.nShape;
1525  return;
1526  }
1527 
1528  /* this flag set true if we did not parse a continuum source, thus creating
1529  * the arrays that are needed - return at this point */
1530  if( lgNoContinuum )
1531  return;
1532 
1533  // must have at least 2 points to interpolate
1534  ASSERT( rfield.ncont[rfield.nShape] > 1 );
1535 
1536  ASSERT( rfield.tNu[rfield.nShape].size() == (size_t)rfield.ncont[rfield.nShape] );
1537  ASSERT( rfield.tslop[rfield.nShape].size() == (size_t)rfield.ncont[rfield.nShape] );
1538 
1539 
1540  ASSERT( rfield.tNu[rfield.nShape][0].Ryd() > 0. );
1541 
1542  /* tFluxLog will be log(fnu) at that spot, read into tslop */
1543  for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
1544  {
1545  rfield.tFluxLog[rfield.nShape].push_back( rfield.tslop[rfield.nShape][i] );
1546  }
1547 
1548  ASSERT( rfield.tFluxLog[rfield.nShape].size() == (size_t)rfield.ncont[rfield.nShape] );
1549 
1550  /* at this point tslop is the log of the intensity */
1551  for( i=0; i < rfield.ncont[rfield.nShape]-1; i++ )
1552  {
1553  rfield.tslop[rfield.nShape][i] =
1554  (realnum)((rfield.tslop[rfield.nShape][i+1] -
1555  rfield.tslop[rfield.nShape][i])/
1556  log10(rfield.tNu[rfield.nShape][i+1].Ryd()/
1557  rfield.tNu[rfield.nShape][i].Ryd()));
1558  }
1560 
1561  if( trace.lgConBug && trace.lgTrace )
1562  {
1563  fprintf( ioQQQ, " Table for this continuum; TNU, TFAC, TSLOP\n" );
1564  for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
1565  {
1566  fprintf( ioQQQ, "%12.4e %12.4e %12.4e\n", rfield.tNu[rfield.nShape][i].Ryd(),
1568  }
1569  }
1570 
1571  /* renormalize continuum so that quantity passes through unity at 1 Ryd
1572  * lgHit will be set false when we get a hit in following loop,
1573  * then test made to make sure it happened*/
1574  lgHit = false;
1575  /*following will be reset when hit occurs*/
1576  fac = -DBL_MAX;
1577  /* >>chng 04 mar 16, chng loop so breaks when hit, previously had cycled
1578  * until last point reached, so last point always used */
1579  for( i=0; i < rfield.ncont[rfield.nShape] && !lgHit; i++ )
1580  {
1581  if( rfield.tNu[rfield.nShape][i].Ryd() > 1. )
1582  {
1583  fac = rfield.tFluxLog[rfield.nShape][i];
1584  lgHit = true;
1585  }
1586  }
1587 
1588  if( lgHit )
1589  {
1590  /* do the renormalization */
1591  for( i=0; i < rfield.ncont[rfield.nShape] ; i++ )
1592  rfield.tFluxLog[rfield.nShape][i] -= (realnum)fac;
1593  }
1594 
1595  ++rfield.nShape;
1596 
1597  return;
1598 }
1599 
1600 
1601 
1602 /*ZeroContin store sets of built in continua */
1603 STATIC void ZeroContin(void)
1604 {
1605 
1606  DEBUG_ENTRY( "ZeroContin()" );
1607 
1608  /* Draine 1978 ISM continuum */
1609  /* freq in ryd */
1610  tnudrn[0] = 0.3676;
1611  tnudrn[1] = 0.4144;
1612  tnudrn[2] = 0.4558;
1613  tnudrn[3] = 0.5064;
1614  tnudrn[4] = 0.5698;
1615  tnudrn[5] = 0.6511;
1616  tnudrn[6] = 0.7012;
1617  tnudrn[7] = 0.7597;
1618  tnudrn[8] = 0.8220;
1619  tnudrn[9] = 0.9116;
1620  tnudrn[10] = 0.9120;
1621  tnudrn[11] = 0.9306;
1622  tnudrn[12] = 0.9600;
1623  tnudrn[13] = 0.9806;
1624  /* >>chng 05 aug 03, this energy is so close to 1 ryd that it spills over
1625  * into the H-ionizing continuum - move down to 0.99 */
1626  /* >>chng 05 aug 08, this destabilized pdr_leiden_hack_v4 (!) so put back to
1627  * original value and include extinguish command */
1628  tnudrn[14] = 0.9999;
1629 
1630  /* f_nu from equation 23 of
1631  * >>refer ism field Draine, B.T. & Bertoldi, F. 1996, ApJ, 468, 269 */
1632  tsldrn[0] = -17.8063;
1633  tsldrn[1] = -17.7575;
1634  tsldrn[2] = -17.7268;
1635  tsldrn[3] = -17.7036;
1636  tsldrn[4] = -17.6953;
1637  tsldrn[5] = -17.7182;
1638  tsldrn[6] = -17.7524;
1639  tsldrn[7] = -17.8154;
1640  tsldrn[8] = -17.9176;
1641  tsldrn[9] = -18.1675;
1642  tsldrn[10] = -18.1690;
1643  tsldrn[11] = -18.2477;
1644  tsldrn[12] = -18.4075;
1645  tsldrn[13] = -18.5624;
1646  tsldrn[14] = -18.7722;
1647 
1648  /* generic AGN continuum taken from
1649  * >>refer AGN cont Mathews and Ferland ApJ Dec15 '87
1650  * except self absorption at 10 microns, nu**-2.5 below that */
1651  tnuagn[0] = Energy(1e-5);
1652  tnuagn[1] = Energy(9.12e-3);
1653  tnuagn[2] = Energy(.206);
1654  tnuagn[3] = Energy(1.743);
1655  tnuagn[4] = Energy(4.13);
1656  tnuagn[5] = Energy(26.84);
1657  tnuagn[6] = Energy(7353.);
1658  tnuagn[7] = Energy(7.4e6);
1659 
1660  tslagn[0] = -3.388;
1661  tslagn[1] = 4.0115;
1662  tslagn[2] = 2.6576;
1663  tslagn[3] = 2.194;
1664  tslagn[4] = 1.819;
1665  tslagn[5] = -.6192;
1666  tslagn[6] = -2.326;
1667  tslagn[7] = -7.34;
1668  resetBltin( tnuagn , tslagn , true );
1669 
1670  /* interstellar radiation field from Black 1987, "Interstellar Processes"
1671  * table of log nu, log nu*fnu taken from his figure 2 */
1672  /* >>chng 99 jun 14 energy range lowered to 1e-8 ryd */
1673  tnuism[0] = 6.00;
1674  /*tnuism[0] = 9.00;*/
1675  tnuism[1] = 10.72;
1676  tnuism[2] = 11.00;
1677  tnuism[3] = 11.23;
1678  tnuism[4] = 11.47;
1679  tnuism[5] = 11.55;
1680  tnuism[6] = 11.85;
1681  tnuism[7] = 12.26;
1682  tnuism[8] = 12.54;
1683  tnuism[9] = 12.71;
1684  tnuism[10] = 13.10;
1685  tnuism[11] = 13.64;
1686  tnuism[12] = 14.14;
1687  tnuism[13] = 14.38;
1688  tnuism[14] = 14.63;
1689  tnuism[15] = 14.93;
1690  tnuism[16] = 15.08;
1691  tnuism[17] = 15.36;
1692  tnuism[18] = 15.43;
1693  tnuism[19] = 16.25;
1694  tnuism[20] = 17.09;
1695  tnuism[21] = 18.00;
1696  tnuism[22] = 23.00;
1697  /* these are log nu*Fnu */
1698  fnuism[0] = -16.708;
1699  /*fnuism[0] = -7.97;*/
1700  fnuism[1] = -2.96;
1701  fnuism[2] = -2.47;
1702  fnuism[3] = -2.09;
1703  fnuism[4] = -2.11;
1704  fnuism[5] = -2.34;
1705  fnuism[6] = -3.66;
1706  fnuism[7] = -2.72;
1707  fnuism[8] = -2.45;
1708  fnuism[9] = -2.57;
1709  fnuism[10] = -3.85;
1710  fnuism[11] = -3.34;
1711  fnuism[12] = -2.30;
1712  fnuism[13] = -1.79;
1713  fnuism[14] = -1.79;
1714  fnuism[15] = -2.34;
1715  fnuism[16] = -2.72;
1716  fnuism[17] = -2.55;
1717  fnuism[18] = -2.62;
1718  fnuism[19] = -5.68;
1719  fnuism[20] = -6.45;
1720  fnuism[21] = -6.30;
1721  fnuism[22] = -11.3;
1722 
1723  return;
1724 }
1725 
1726 /*lines_table invoked by table lines command, check if we can find all lines in a given list */
1728 {
1729  DEBUG_ENTRY( "lines_table()" );
1730 
1731  if( chLINE_LIST.empty() )
1732  return 0;
1733 
1734  vector<string> chLabel;
1735  vector<realnum> wl;
1736 
1737  long nLINE_TABLE = cdGetLineList( chLINE_LIST.c_str(), chLabel, wl );
1738 
1739  // the check if the file exists has already been done by the parser
1740  if( nLINE_TABLE == 0 )
1741  return 0;
1742 
1743  fprintf( ioQQQ , "lines_table checking lines within data table %s\n", chLINE_LIST.c_str() );
1744  long miss = 0;
1745 
1746  /* \todo 2 DOS carriage return on linux screws this up. Can we overlook the CR? Skip in cdgetlinelist? */
1747  for( long n=0; n < nLINE_TABLE; ++n )
1748  {
1749  if( LineSave.findline( chLabel[n].c_str(), wl[n] ) <= 0 )
1750  {
1751  ++miss;
1752  fprintf(ioQQQ,"lines_table in parse_table.cpp did not find line %s ",chLabel[n].c_str());
1753  prt_wl(ioQQQ,wl[n]);
1754  fprintf(ioQQQ,"\n");
1755  }
1756  }
1757  if( miss )
1758  {
1759  /* this is so that we pick up problem automatically */
1760  fprintf( ioQQQ , " BOTCHED MONITORS!!! Botched Monitors!!! lines_table could not find a total of %li lines\n\n", miss );
1761  }
1762  else
1763  {
1764  fprintf( ioQQQ , "lines_table found all lines\n\n" );
1765  }
1766 
1767  // deallocate the memory allocated in cdGetLineList()
1768  chLabel.clear();
1769 
1770  return miss;
1771 }
1772 
1773 /*ReadTable called by TABLE READ to read in continuum from SAVE TRANSMITTED CONTINUUM */
1774 STATIC void ReadTable(const char *fnam)
1775 {
1776  char chLine[INPUT_LINE_LENGTH];
1777  size_t i;
1778  FILE *io;
1779 
1780  DEBUG_ENTRY( "ReadTable()" );
1781 
1782  io = open_data( fnam, "r", AS_LOCAL_ONLY );
1783 
1784  /* skip comment lines at the start of the file */
1785  do
1786  {
1787  if( read_whole_line( chLine, (int)sizeof(chLine), io ) == NULL )
1788  {
1789  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1790  goto error;
1791  }
1792  }
1793  while( chLine[0] == '#' );
1794 
1795  /* line 3: the version number */
1796  long version;
1797  sscanf( chLine, "%ld", &version );
1798  if( version != VERSION_TRNCON )
1799  {
1800  fprintf( ioQQQ,
1801  " the input continuum file has the wrong version number, I expected %li and found %li.\n",
1802  VERSION_TRNCON, version);
1803  goto error;
1804  }
1805 
1806  char md5sum[NMD5];
1807  long nflux;
1808  double mesh_lo, mesh_hi;
1809  union {
1810  double x;
1811  uint32 i[2];
1812  } u;
1813 
1814  /* line 4: the MD5 checksum */
1815  if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
1816  {
1817  for( unsigned int i=0; i < NMD5; ++i )
1818  md5sum[i] = chLine[i];
1819  }
1820  else
1821  {
1822  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1823  goto error;
1824  }
1825 
1826  /* line 5: the lower limit of the energy grid */
1827  if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
1828  {
1829  if( cpu.i().big_endian() )
1830  sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
1831  else
1832  sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
1833  mesh_lo = u.x;
1834  }
1835  else
1836  {
1837  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1838  goto error;
1839  }
1840 
1841  /* line 6: the upper limit of the energy grid */
1842  if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
1843  {
1844  if( cpu.i().big_endian() )
1845  sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
1846  else
1847  sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
1848  mesh_hi = u.x;
1849  }
1850  else
1851  {
1852  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1853  goto error;
1854  }
1855 
1856  if( strncmp( md5sum, rfield.mesh_md5sum().c_str(), NMD5 ) != 0 ||
1857  !fp_equal_tol( mesh_lo, rfield.emm(), 1.e-11*rfield.emm() ) ||
1858  !fp_equal_tol( mesh_hi, rfield.egamry(), 1.e-7*rfield.egamry() ) )
1859  {
1860  fprintf( ioQQQ, " the input continuum file has an incompatible energy grid.\n" );
1861  goto error;
1862  }
1863 
1864  /* line 7: the energy grid resolution scale factor */
1865  if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
1866  {
1867  if( cpu.i().big_endian() )
1868  sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
1869  else
1870  sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
1871  rfield.RSFCheck[rfield.nShape] = u.x;
1872  }
1873  else
1874  {
1875  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1876  goto error;
1877  }
1878 
1879  /* line 8: the outer radius of the saved model */
1880  if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
1881  {
1882  sscanf( chLine, "%le", &rfield.TableRadius[rfield.nShape] );
1883  }
1884  else
1885  {
1886  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1887  goto error;
1888  }
1889 
1890  /* line 9: the number of frequency grid points contained in the file */
1891  if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
1892  {
1893  sscanf( chLine, "%ld", &nflux );
1894  }
1895  else
1896  {
1897  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1898  goto error;
1899  }
1900 
1901  /* line 10: empty comment line */
1902  if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
1903  {
1904  fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
1905  goto error;
1906  }
1907 
1908  /* now read in the file of numbers */
1909  i = 0;
1910  /* keep reading until we hit eol or run out of room in the array */
1911  while( read_whole_line(chLine, (int)sizeof(chLine),io) != NULL )
1912  {
1913  double help[2];
1914  sscanf( chLine, "%lf%lf ", &help[0], &help[1] );
1915  // check that frequency grid matches, frequencies are printed with 6 significant digits
1916  if( !fp_equal_tol(help[0],rfield.anu(i),1e-5*rfield.anu(i)) )
1917  {
1918  fprintf(ioQQQ,"\n\nPROBLEM in TABLE READ continuum frequencies: point %li should "
1919  "have %e and has %e\n", (unsigned long)i, rfield.anu(i), help[0] );
1920  fprintf(ioQQQ,"Please recreate this file.\n" );
1922  }
1923  rfield.tNu[rfield.nShape].push_back( Energy(help[0]) );
1924  // Convert to standard FluxLog units
1925  if (help[1] > 0.0)
1926  {
1927  rfield.tFluxLog[rfield.nShape].push_back( (realnum)log10(help[1]/help[0]) );
1928  }
1929  else
1930  {
1931  rfield.tFluxLog[rfield.nShape].push_back( realnum(-70.) );
1932  }
1933  ++i;
1934  }
1935  /* put pointer at last good value */
1936  rfield.ncont[rfield.nShape] = i;
1937 
1938  ASSERT( rfield.tNu[rfield.nShape].size() == (size_t)rfield.ncont[rfield.nShape] );
1939  ASSERT( rfield.tFluxLog[rfield.nShape].size() == (size_t)rfield.ncont[rfield.nShape] );
1940 
1941  /* check that sane number of values entered */
1942  if( rfield.ncont[rfield.nShape] != nflux )
1943  {
1944  fprintf( ioQQQ, " the input continuum file has the wrong number of points: %ld\n",
1946  goto error;
1947  }
1948 
1949  fclose(io);
1950  return;
1951 
1952 error:
1953  fprintf( ioQQQ, " please recreate this file using the SAVE TRANSMITTED CONTINUUM command.\n" );
1955 }
double emm() const
Definition: mesh.h:93
bool nMatch(const char *chKey) const
Definition: parser.h:150
static bool lgCalled
Definition: cddrive.cpp:424
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
bool lgBeamed[LIMSPC]
Definition: rfield.h:294
long RauchInterpolateHydr(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1295
static double tsldrn[NDRAINE]
Definition: parse_table.cpp:51
double FFmtRead(void)
Definition: parser.cpp:472
static const int NDRAINE
Definition: parse_table.cpp:50
double exp10(double x)
Definition: cddefines.h:1368
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:296
T * get_ptr(T *v)
Definition: cddefines.h:1109
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
t_input input
Definition: input.cpp:12
const double BIGDOUBLE
Definition: cpu.h:249
long int nvfpnt[LIMPAR]
Definition: optimize.h:198
bool lgSphericalDilution[LIMSPC]
Definition: rfield.h:321
STATIC void ReadTable(const char *fnam)
long findline(const char *chLabel, realnum wavelength)
Definition: lines.cpp:294
long int m_nqh
Definition: parser.h:54
t_cpu_i & i()
Definition: cpu.h:419
double totpow[LIMSPC]
Definition: rfield.h:284
string mesh_md5sum() const
Definition: mesh.h:112
long CoStarInterpolate(double val[], long *nval, long *ndim, IntMode imode, bool lgHalo, bool lgList, double *val0_lo, double *val0_hi)
Definition: stars.cpp:653
static const long VERSION_TRNCON
Definition: save.h:16
char chRSpec[LIMSPC][5]
Definition: rfield.h:335
long TlustyInterpolate(double val[], long *nval, long *ndim, tl_grid tlg, const char *chMetalicity, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1590
void set(double energy)
Definition: energy.h:26
static const int NAGN
Definition: parse_table.cpp:45
void ParseTable(Parser &p)
Definition: parse_table.cpp:94
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:525
long HaardtMadauInterpolate(double val, int version, bool lgQuasar, double *zlow, double *zhigh)
Definition: stars.cpp:799
static const double tnuHM96[NHM96]
Definition: parse_table.cpp:35
int GetQuote(string &chLabel)
Definition: parser.cpp:213
realnum varang[LIMPAR][2]
Definition: optimize.h:201
long int nRead
Definition: input.h:105
static double tnuism[NISM]
Definition: parse_table.cpp:27
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:904
t_LineSave LineSave
Definition: lines.cpp:10
double TableRadius[LIMSPC]
Definition: rfield.h:319
long RauchInterpolateHCa(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1175
double RSFCheck[LIMSPC]
Definition: rfield.h:327
bool nMatchErase(const char *chKey)
Definition: parser.h:173
vector< Energy > tNu[LIMSPC]
Definition: rfield.h:314
vector< realnum > tFluxLog[LIMSPC]
Definition: rfield.h:317
static double tnudrn[NDRAINE]
Definition: parse_table.cpp:51
long int cdGetLineList(const char chFile[], vector< string > &chLabels, vector< realnum > &wl)
bool big_endian() const
Definition: cpu.h:354
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
long WMBASICInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1763
FILE * ioQQQ
Definition: cddefines.cpp:7
static const int NHM96
Definition: parse_table.cpp:33
long ncont[LIMSPC]
Definition: rfield.h:323
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:192
vector< realnum > tslop[LIMSPC]
Definition: rfield.h:315
static realnum tslagn[NAGN]
Definition: parse_table.cpp:47
Definition: parser.h:43
double anu(size_t i) const
Definition: mesh.h:120
bool lgVarOn
Definition: optimize.h:207
bool lgTimeVary[LIMSPC]
Definition: rfield.h:290
static string chLINE_LIST
Definition: parse_table.cpp:23
double range[LIMSPC][2]
Definition: rfield.h:331
IntMode
Definition: stars.h:16
const int LIMSPC
Definition: rfield.h:21
t_trace trace
Definition: trace.cpp:5
long KhaireSrianandInterpolate(double val, int Q, double *zlow, double *zhigh)
Definition: stars.cpp:872
size_t ipointC(double anu) const
Definition: mesh.h:161
static const int NISM
Definition: parse_table.cpp:26
bool lgConBug
Definition: trace.h:97
long RauchInterpolateHpHe(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1351
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
long MihalasInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:973
long int nparm
Definition: optimize.h:204
static Energy tnuagn[NAGN]
Definition: parse_table.cpp:46
t_rfield rfield
Definition: rfield.cpp:9
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
#define MDIM
Definition: stars.h:9
Illuminate::IlluminationType Illumination[LIMSPC]
Definition: rfield.h:300
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
#define cdEXIT(FAIL)
Definition: cddefines.h:482
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:345
static const double fnuHM96[NHM96]
Definition: parse_table.cpp:41
bool lgMustBlockHIon
Definition: rfield.h:92
char * chLabel
Definition: species.h:38
t_optimize optimize
Definition: optimize.cpp:6
STATIC void ZeroContin(void)
char chSpNorm[LIMSPC][5]
Definition: rfield.h:335
realnum vincr[LIMPAR]
Definition: optimize.h:195
long RauchInterpolateCOWD(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1267
#define ASSERT(exp)
Definition: cddefines.h:613
const long LIMEXT
Definition: optimize.h:60
const char * StandardEnergyUnit(const char *chCard)
Definition: energy.cpp:44
long WernerInterpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1690
double Ryd() const
Definition: energy.h:33
Definition: energy.h:9
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
char chDelimiter[3]
Definition: input.h:83
double egamry() const
Definition: mesh.h:97
static const unsigned int NMD5
Definition: thirdparty.h:451
long RauchInterpolateHelium(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1323
bool lgEOL(void) const
Definition: parser.h:113
long GridInterpolate(double val[], long *nval, long *ndim, const char *FileName, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:760
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
Definition: stars.h:22
void AtmospheresAvail()
Definition: stars.cpp:254
long int nShape
Definition: rfield.h:306
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
STATIC void resetBltin(Energy *tnu, realnum *fluxlog, bool lgLog)
Definition: parse_table.cpp:61
void caps(char *chCard)
Definition: service.cpp:295
long Kurucz79Interpolate(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:921
static t_cpu cpu
Definition: cpu.h:427
long RauchInterpolatePG1159(double val[], long *nval, long *ndim, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1239
long RauchInterpolateHNi(double val[], long *nval, long *ndim, bool lgHalo, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:1207
long int nvarxt[LIMPAR]
Definition: optimize.h:198
char chSpType[LIMSPC][6]
Definition: rfield.h:335
int lines_table()
Definition: stars.h:22
long AtlasInterpolate(double val[], long *nval, long *ndim, const char *chMetalicity, const char *chODFNew, bool lgList, double *Tlow, double *Thigh)
Definition: stars.cpp:553
static double fnuism[NISM]
Definition: parse_table.cpp:27
#define EXIT_SUCCESS
Definition: cddefines.h:166
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
#define UNUSED
Definition: cpu.h:14