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