cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
monitor_results.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 /*InitMonitorResults, this must be called first, done at startup of ParseCommands*/
4 /*ParseMonitorResults - parse input stream */
5 /*lgCheckMonitors, checks asserts, last thing cloudy calls, returns true if all are ok, false if problems */
6 #include "cddefines.h"
7 #include "input.h"
8 #include "conv.h"
9 #include "optimize.h"
10 #include "iso.h"
11 #include "called.h"
12 #include "atmdat.h"
13 #include "hcmap.h"
14 #include "thermal.h"
15 #include "pressure.h"
16 #include "struc.h"
17 #include "wind.h"
18 #include "h2.h"
19 #include "colden.h"
20 #include "dense.h"
21 #include "lines.h"
22 #include "secondaries.h"
23 #include "radius.h"
24 #include "version.h"
25 #include "hmi.h"
26 #include "prt.h"
27 #include "grainvar.h"
28 #include "cddrive.h"
29 #include "elementnames.h"
30 #include "monitor_results.h"
31 #include "parser.h"
32 #include "mole.h"
33 #include "rfield.h"
34 #include "lines_service.h"
35 #include "service.h"
36 #include "generic_state.h"
37 
38 bool lgMonitorsOK , lgBigBotch, lgPrtSciNot;
39 t_monitorresults MonitorResults;
40 
41 /* flag to remember that InitMonitorResults was called */
42 static bool lgInitDone=false ,
43  /* will be set true when space for asserts is allocated */
44  lgSpaceAllocated=false;
45 
46 /* number of asserts we can handle, used in dim of space */
47 static const int NASSERTS = 450;
48 
49 /* default relative error for monitored physical quantities */
50 static realnum ErrorDefault;
51 
52 /* default relative error for monitored performance metrics */
53 static realnum ErrorDefaultPerformance;
54 
55 static vector<string> chAssertLineLabel;
56 
57 /* type of line integral to report:
58  * 0: intrinsic, 1: emergent, 2: cumul intrinsic, 3: cumul emergent */
59 static vector<int> iLineType;
60 
61 static vector<std::vector<TransitionProxy>* > assertBlends;
62 
63 /* this will be = for equal, < or > for limit */
64 static vector<char> chAssertLimit;
65 
66 /* this will be a two character label identifying which type of monitor */
67 static char **chAssertType;
68 
69 /* the values and error in the asserted quantity */
70 static vector<double> AssertQuantity, AssertQuantity2 ,AssertError;
71 static multi_arr<double,2> Param;
72 
73 /* this flag says where we print linear or log quantity */
74 static vector<bool> lgQuantityLog;
75 static long nAsserts=0;
76 static vector<realnum> wavelength;
77 
78 static vector<string> strAssertSpecies;
79 
80 void PrtOneMonitor( FILE *ioMONITOR, const char* chAssertType, const string chAssertLineLabel,
81  const realnum wavelength, const int iLineType, const double PredQuan, const char chAssertLimit, const double AssertQuantity,
82  const double RelError, const double AssertError, const bool lg1OK, const bool lgQuantityLog, const bool lgFound );
83 
84 inline double ForcePass( char chAssertLimit1 )
85 {
86  // force monitors to pass by returning a safe value
87  if( chAssertLimit1 == '=' )
88  return 0.;
89  else if( chAssertLimit1 == '<' )
90  return 1.;
91  else if( chAssertLimit1 == '>' )
92  return -1.;
93  else
94  TotalInsanity();
95 }
96 
97 /*======================================================================*/
98 /*InitMonitorResults, this must be called first, done at startup of ParseCommands*/
99 void InitMonitorResults(void)
100 {
101  /* set flag that init was called, and set number of asserts to zero.
102  * this is done by ParseComments for every model, even when no asserts will
103  * be done, so do not allocate space at this time */
104  lgInitDone = true;
105  nAsserts = 0;
106 
107  // following occur in hazy1 default for monitor commands
108  // default error, changed with MONITOR SET ERROR
109  ErrorDefault = 0.05f;
110  // default performance monitor, changed with MONITOR SET PERFORMANCE ERROR
111  ErrorDefaultPerformance = 0.2f;
112 }
113 
114 /*======================================================================*/
115 /*ParseMonitorResults parse the monitor command */
116 void ParseMonitorResults(Parser &p)
117 {
118  long nelem,
119  n2;
120 
121  DEBUG_ENTRY( "ParseMonitorResults()" );
122 
123  if( !lgInitDone )
124  {
125  fprintf( ioQQQ, " ParseMonitorResults called before InitAsserResults\n" );
126  cdEXIT(EXIT_FAILURE);
127  }
128 
129  /* has space been allocated yet? */
130  if( !lgSpaceAllocated )
131  {
132  /* - no, we must allocate it */
133  /* remember that space has been allocated */
134  lgSpaceAllocated = true;
135 
136  /* create space for the array of labels*/
137  chAssertLineLabel.resize(NASSERTS);
138 
139  assertBlends.resize(NASSERTS);
140 
141  /* the 2-character string saying what type of monitor */
142  chAssertType = ((char **)MALLOC(NASSERTS*sizeof(char *)));
143 
144  /* these are a pair of optional parameters */
145  Param.reserve( NASSERTS );
146 
147  /* now fill out the 2D arrays */
148  for( long i=0; i<NASSERTS; ++i )
149  {
150  chAssertLineLabel[i] = "unkn" ;
151 
152  assertBlends[i] = NULL;
153 
154  /* two char plus eol */
155  chAssertType[i] = ((char *)MALLOC(3*sizeof(char )));
156  strcpy(chAssertType[i],"un" );
157 
158  /* these are a pair of optional parameters */
159  Param.reserve( i, 5 );
160  }
161 
162  Param.alloc();
163 
164  /* now make space for the asserted quantities */
165  AssertQuantity.resize(NASSERTS);
166 
167  AssertQuantity2.resize(NASSERTS);
168 
169  /* now the errors */
170  AssertError.resize(NASSERTS);
171 
172  /* now the line wavelengths */
173  wavelength.resize(NASSERTS);
174 
175  /* label for a species */
176  strAssertSpecies.resize(NASSERTS);
177 
178  /* now the flag saying whether should be log */
179  lgQuantityLog.resize(NASSERTS);
180 
181  /* the flag for upper, lower limit, or equal */
182  chAssertLimit.resize(NASSERTS);
183 
184  iLineType.resize(NASSERTS, -1);
185  }
186  /* end space allocation - we are ready to roll */
187 
188  /* read asserted values from file if GRID keyword is present */
189  vector<double> AssertVector;
190  if( p.nMatch( "GRID" ) )
191  {
192  ASSERT( optimize.nOptimiz >= 0 );
193  AssertVector.resize( optimize.nOptimiz+1 );
194 
195  /* this file should contain all the values that are asserted */
196  string chLabel;
197  if (p.GetQuote( chLabel ))
198  p.StringError();
199 
200  bool lgEOF;
201  input_readvector( chLabel.c_str(), get_ptr(AssertVector), optimize.nOptimiz+1, &lgEOF );
202  if( lgEOF )
203  fprintf(ioQQQ,"PROBLEM the file %s does not have enough values. sorry\n", chLabel.c_str() );
204  }
205 
206  /* first say we do not know what job to do */
207  strcpy( chAssertType[nAsserts] , "un" );
208 
209  /* false means print linear quantity - will be set true if entered
210  * quantity comes in as a log */
211  lgQuantityLog[nAsserts] = false;
212 
213  /* all asserts have option for quantity to be a limit, or the quantity itself */
214  if( p.nMatch("<" ) )
215  {
216  chAssertLimit[nAsserts] = '<';
217  }
218  else if( p.nMatch(">" ) )
219  {
220  chAssertLimit[nAsserts] = '>';
221  }
222  else
223  {
224  chAssertLimit[nAsserts] = '=';
225  }
226 
227  /* which quantity will we check?, first is */
228 
229  if( p.nMatch(" SET" ) )
230  {
231  /* set an option for the monitor command, not an actual asserted
232  * quantity */
233 
234  /* decrement number of asserts since will be incremented below,
235  * this is not an actual asserted quantity */
236  if( nAsserts >0 )
237  fprintf(ioQQQ," The default monitor error is being changed after"
238  " some asserts were entered. \n This will only affect asserts"
239  " that come after this command.\n");
240  --nAsserts;
241 
242  if( p.nMatch("ERRO" ) )
243  {
244  if( p.nMatch("PERF" ) )
245  {
246  // set performance monitor error
247  ErrorDefaultPerformance =
248  (realnum)p.FFmtRead();
249  if( p.lgEOL() )
250  p.NoNumb("error");
251  }
252  /* set default error */
253  ErrorDefault =
254  (realnum)p.FFmtRead();
255  if( p.lgEOL() )
256  p.NoNumb("error");
257  }
258  else
259  {
260  /* problem - no recognized quantity */
261  fprintf( ioQQQ,
262  " I could not identify an option on this ASSERT SET XXX command.\n");
263  fprintf( ioQQQ, " Sorry.\n" );
264  cdEXIT(EXIT_FAILURE);
265  }
266  }
267 
268  /* monitor mean ionization */
269  else if( p.nMatch("IONI" ) )
270  {
271 
272  /* say that this will be mean ionization fraction */
273 
274  /* f will indicate average over radius, F over volume -
275  * check whether keyword radius or volume occurs,
276  * default will be radius */
277  if( p.nMatch("VOLU" ) )
278  {
279  strcpy( chAssertType[nAsserts] , "F " );
280  }
281  else
282  {
283  /* this is default case, Fraction over radius */
284  strcpy( chAssertType[nAsserts] , "f " );
285  }
286 
287  /* first get element label and make null terminated string*/
288  if( (nelem = p.GetElem()) < 0 )
289  {
290  fprintf( ioQQQ,
291  " I could not identify an element name on this line.\n");
292  fprintf( ioQQQ, " Sorry.\n" );
293  cdEXIT(EXIT_FAILURE);
294  }
295  ASSERT( nelem>= 0);
296  ASSERT( nAsserts>= 0);
297  /* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem])
298  * into array that will be used to get ionization after calculation */
299  chAssertLineLabel[nAsserts] = elementnames.chElementNameShort[nelem];
300 
301  /* now get ionization stage, which will be saved into wavelength */
302  wavelength[nAsserts] = (realnum)p.FFmtRead();
303  if( p.lgEOL() )
304  {
305  p.NoNumb("ionization stage");
306  }
307  /* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */
308  if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 )
309  {
310  fprintf( ioQQQ,
311  " The ionization stage is inappropriate for this element.\n");
312  fprintf( ioQQQ, " Sorry.\n" );
313  cdEXIT(EXIT_FAILURE);
314  }
315 
316  if( p.nMatch( "GRID" ) )
317  {
318  AssertQuantity[nAsserts] = AssertVector[optimize.nOptimiz];
319  }
320  else
321  {
322  /* now get ionization fraction, log if number is negative or == 0,
323  * linear if positive but less than or equal to 1.*/
324  AssertQuantity[nAsserts] = p.FFmtRead();
325  if( p.lgEOL() )
326  p.NoNumb("ionization fraction");
327  }
328 
329  /* optional error, default available */
330  AssertError[nAsserts] = p.FFmtRead();
331  if( p.lgEOL() )
332  AssertError[nAsserts] = ErrorDefault;
333 
334  /* now make sure we end up with positive linear ionization fraction that
335  * is greater than 0 but less than or equal to 1. */
336  if( AssertQuantity[nAsserts] <= 0. )
337  {
338  /* log since negative or 0 */
339  AssertQuantity[nAsserts] =
340  exp10(AssertQuantity[nAsserts] );
341  /* entered as a log, so print as a log too */
342  lgQuantityLog[nAsserts] = true;
343  }
344  else if( AssertQuantity[nAsserts] > 1. )
345  {
346  fprintf( ioQQQ,
347  " The ionization fraction must be less than one.\n");
348  fprintf( ioQQQ, " Sorry.\n" );
349  cdEXIT(EXIT_FAILURE);
350  }
351 
352  /* result cannot be zero */
353  if( fabs(AssertQuantity[nAsserts]) <= SMALLDOUBLE )
354  {
355  fprintf( ioQQQ,
356  " The ionization ionization fraction is too small, or zero. Check input\n");
357  fprintf( ioQQQ, " Sorry.\n" );
358  cdEXIT(EXIT_FAILURE);
359  }
360  }
361 
362  /* molecular fraction averaged over model */
363  else if( p.nMatch("MOLE" )&&p.nMatch("FRAC" ) )
364  {
365 
366  /* say that this will be mean molecular fraction */
367 
368  /* mf will indicate average over radius, MF over vol -
369  * check whether keyword radius or volume occurs,
370  * default will be radius */
372  if( p.nMatch("VOLU" ) )
373  {
374  strcpy( chAssertType[nAsserts] , "MF" );
375  }
376  else
377  {
378  /* this is default case, Fraction over radius */
379  strcpy( chAssertType[nAsserts] , "mf" );
380  }
381 
382  if( p.nMatchErase(" H2 " ) )
383  {
384  chAssertLineLabel[nAsserts] = "H2" ;
385  /* increment to get past the label */
386  }
387  else if( p.nMatch(" CO " ) )
388  {
389  chAssertLineLabel[nAsserts] = "CO" ;
390  }
391  else
392  {
393  fprintf( ioQQQ,
394  " I could not identify CO or H2 on this line.\n");
395  fprintf( ioQQQ, " Sorry.\n" );
396  cdEXIT(EXIT_FAILURE);
397  }
398 
399  /* not meaningful */
400  wavelength[nAsserts] = 0;
401 
402  /* now get log of molecular fraction */
403  AssertQuantity[nAsserts] = p.FFmtRead();
404  if( p.lgEOL() )
405  {
406  p.NoNumb("molecular fraction");
407  }
408  if( AssertQuantity[nAsserts] <= 0. )
409  {
410  /* if negative then entered as log, but we will compare with linear */
411  AssertQuantity[nAsserts] =
412  exp10(AssertQuantity[nAsserts] );
413  }
414 
415  /* optional error, default available (cannot do before loop since we
416  * do not know how many numbers are on line */
417  AssertError[nAsserts] = p.FFmtRead();
418  if( p.lgEOL() )
419  {
420  /* default error was set above */
421  AssertError[nAsserts] = ErrorDefault;
422  }
423  /* print results as logs */
424  lgQuantityLog[nAsserts] = true;
425  }
426 
427  /* monitor line "LINE" -- key is ine_ since linear option appears on some commands */
428  else if( p.nMatch(" LINE " ) )
429  {
430  if( p.nMatch("LUMI") || p.nMatch("INTE"))
431  {
432  /* say that this is a line luminosity or intensity*/
433  strcpy( chAssertType[nAsserts] , "Ll" );
434 
435  /* entered as a log, so print as a log too */
436  lgQuantityLog[nAsserts] = true;
437  }
438  else if( p.nMatch( "CASE" ) && p.nMatch( " B ") )
439  {
440  /* say that this is line relative to norm line - this is the default */
441  strcpy( chAssertType[nAsserts] , "Lb" );
442 
443  /* entered linear quantity, so print as linear too */
444  lgQuantityLog[nAsserts] = false;
445  }
446  else
447  {
448  /* say that this is line relative to norm line - this is the default */
449  strcpy( chAssertType[nAsserts] , "Lr" );
450 
451  /* entered linear quantity, so print as linear too */
452  lgQuantityLog[nAsserts] = false;
453  }
454 
455  iLineType[nAsserts] = 0;
456 
457  if( p.nMatch( "EMER" ) )
458  {
459  iLineType[nAsserts] = 1;
460  }
461 
462  if( p.nMatch( "CUMU" ) )
463  {
464  iLineType[nAsserts] += 2;
465  }
466 
467  /* this will check a line intensity, first get labels within quotes,
468  * will be null terminated */
469  string chLabel;
470  if (p.GetQuote( chLabel ))
471  p.StringError();
472  trimTrailingWhiteSpace( chLabel );
473 
474  blend_iterator b=blends.find(chLabel);
475  if(0 && b != blends.end())
476  {
477  assertBlends[nAsserts] = &(b->second);
478  }
479  else if( chLabel.length() > NCHLAB-1 )
480  {
481  fprintf( ioQQQ, " The label must be no more than %d char long, between double quotes.\n",
482  NCHLAB-1 );
483  fprintf( ioQQQ, " Sorry.\n" );
484  cdEXIT(EXIT_FAILURE);
485  }
486 
487  /* copy string into array */
488  chAssertLineLabel[nAsserts] = chLabel;
489 
490  /* now get line wavelength */
491  wavelength[nAsserts] = (realnum)p.getWave();
492 
493  /* now get intensity or luminosity -
494  * rel intensity is linear and intensity or luminosity are log */
495  AssertQuantity[nAsserts] = p.FFmtRead();
496  if( p.lgEOL() )
497  {
498  p.NoNumb("intensity/luminosity");
499  }
500  /* luminosity was entered as a log */
501  if( lgQuantityLog[nAsserts] )
502  {
503  if( AssertQuantity[nAsserts] > DBL_MAX_10_EXP ||
504  AssertQuantity[nAsserts] < -DBL_MAX_10_EXP )
505  {
506  fprintf( ioQQQ,
507  " The asserted quantity is a log, but is too large or small, value is %e.\n",
508  AssertQuantity[nAsserts] );
509  fprintf( ioQQQ, " I would crash if I tried to use it.\n Sorry.\n" );
510  cdEXIT(EXIT_FAILURE);
511  }
512  AssertQuantity[nAsserts] =
513  exp10(AssertQuantity[nAsserts] );
514  }
515  if( AssertQuantity[nAsserts]< 0. )
516  {
517  fprintf( ioQQQ,
518  " The relative intensity must be non-negative, and was %e.\n",AssertQuantity[nAsserts] );
519  fprintf( ioQQQ, " Sorry.\n" );
520  cdEXIT(EXIT_FAILURE);
521  }
522 
523  /* optional error, default available */
524  AssertError[nAsserts] = p.FFmtRead();
525  if( p.lgEOL() )
526  {
527  /* default error was set above */
528  AssertError[nAsserts] = ErrorDefault;
529  }
530  }
531 
532  /* monitor line predictions relative to case B */
533  else if( p.nMatch("CASE" ) )
534  {
535  /* this is Case B for some element */
536  strcpy( chAssertType[nAsserts] , "CB" );
537  /* this is relative error */
538  AssertError[nAsserts] = p.FFmtRead();
539  if( p.lgEOL() )
540  /* default error was set above */
541  AssertError[nAsserts] = ErrorDefault;
542  AssertQuantity[nAsserts] = 0;
543  wavelength[nAsserts] = 0.;
544 
545  /* faint option - do not test line if relative intensity is less
546  * than entered value */
547  if( p.GetParam("FAINT",&Param[nAsserts][4]) )
548  {
549  if( p.lgEOL() )
550  {
551  /* did not get 2 numbers */
552  fprintf(ioQQQ," The monitor Case B faint option must have a number,"
553  " the relative intensity of the fainest line to monitor.\n");
554  cdEXIT(EXIT_FAILURE);
555  }
556  /* number is log if <= 0 */
557  if( Param[nAsserts][4]<=0. )
558  Param[nAsserts][4] = exp10( Param[nAsserts][4] );
559  }
560  else
561  {
562  /* use default - include everything*/
563  Param[nAsserts][4] = SMALLFLOAT;
564  }
565 
566  /* range option - to limit check on a certain wavelength range */
567  if( p.GetRange("RANG",&Param[nAsserts][2],&Param[nAsserts][3]) )
568  {
569  if( p.lgEOL() )
570  {
571  /* did not get 2 numbers */
572  fprintf(ioQQQ," The monitor Case B range option must have two numbers,"
573  " the lower and upper limit to the wavelengths in Angstroms.\n");
574  fprintf(ioQQQ," There must be a total of three numbers on the line,"
575  " the relative error followed by the lower and upper limits to the "
576  "wavelength in Angstroms.\n");
577  cdEXIT(EXIT_FAILURE);
578  }
579  if( Param[nAsserts][2]>Param[nAsserts][3])
580  {
581  /* make sure in increasing order */
582  double sav = Param[nAsserts][3];
583  Param[nAsserts][3] = Param[nAsserts][2];
584  Param[nAsserts][2] = sav;
585  }
586  }
587  else
588  {
589  /* use default - include everything*/
590  Param[nAsserts][2] = 0.;
591  Param[nAsserts][3] = 1e30;
592  }
593  /* monitor case b for H - O checking against Hummer & Storey tables */
594  if( p.nMatch("H-LI" ) )
595  {
596  /* H-like - now get an element */
597  if( (nelem = p.GetElem()) < 0 )
598  {
599  /* no name found */
600  fprintf(ioQQQ, "monitor H-like case B did not find an element on this line, sorry\n");
601  p.PrintLine(ioQQQ);
602  cdEXIT(EXIT_FAILURE);
603  }
604  if( nelem>7 )
605  {
606  /* beyond reach of tables */
607  fprintf(ioQQQ, "monitor H-like cannot do elements heavier than O, sorry\n");
608  p.PrintLine(ioQQQ);
609  cdEXIT(EXIT_FAILURE);
610  }
611  Param[nAsserts][0] = ipH_LIKE;
612  Param[nAsserts][1] = nelem;
613  /* generate string to find simple prediction, as in "Ca B" */
614  chAssertLineLabel[nAsserts] = "Ca ";
615  if( p.nMatch("CASE A " ) )
616  chAssertLineLabel[nAsserts] += "A";
617  else
618  chAssertLineLabel[nAsserts] += "B";
619  }
620  else if( p.nMatch("HE-L") )
621  {
622  /* He-like - only helium itself */
623  Param[nAsserts][0] = ipHE_LIKE;
624  Param[nAsserts][1] = ipHELIUM;
625 
626  chAssertLineLabel[nAsserts] = "Ca B";
627  }
628  else
629  {
630  /*no option found */
631  fprintf( ioQQQ,
632  " I could not identify an iso-sequence on this Case A/B command.\n");
633  fprintf( ioQQQ, " Sorry.\n" );
634  cdEXIT(EXIT_FAILURE);
635  }
636  }
637 
638  /* monitor departure coefficients */
639  else if( p.nMatch("DEPA") )
640  {
641  string chLabel;
642  if ( p.GetQuote(chLabel) == 0 )
643  {
644  strAssertSpecies[nAsserts] = chLabel;
645  // printed label must still be less than NCHLAB in length
646  ASSERT(chLabel.length() <= NCHLAB-1);
647  // Pad with spaces
648  chAssertLineLabel[nAsserts] = chLabel;
649  for (long i=chLabel.length();i<NCHLAB-1;++i)
650  chAssertLineLabel[nAsserts] += ' ';
651  }
652 
653  /* get expected average departure coefficient, almost certainly 1 */
654  AssertQuantity[nAsserts] = p.FFmtRead();
655  if( p.lgEOL() )
656  p.NoNumb("average departure coefficient");
657 
658  /* this is relative error, max departure from unity of any level or std */
659  AssertError[nAsserts] = p.FFmtRead();
660  if( p.lgEOL() )
661  /* default error was set above */
662  AssertError[nAsserts] = ErrorDefault;
663 
664  if( !strAssertSpecies[nAsserts].empty() )
665  {
666  // remember this is departure coefficient for some species
667  strcpy( chAssertType[nAsserts] , "DC" );
668  }
669  /* H-like key means do one of the hydrogenic ions */
670  else if( p.nMatch("H-LI" ) )
671  {
672  Param[nAsserts][0] = ipH_LIKE;
673  chAssertLineLabel[nAsserts] = "dHyd" ;
674  /* remember this is departure coefficient for some element */
675  strcpy( chAssertType[nAsserts] , "DI" );
676  /* now get element number for h ion from element name on card */
677  if( (wavelength[nAsserts] = (realnum)p.GetElem()) < 0 )
678  {
679  fprintf( ioQQQ,
680  " I could not identify an element name on this line.\n");
681  fprintf( ioQQQ, " Sorry.\n" );
682  cdEXIT(EXIT_FAILURE);
683  }
684  if( p.nMatch("ZEROOK") )
685  Param[nAsserts][1] = 1.;
686  else
687  Param[nAsserts][1] = 0.;
688  }
689 
690  /* He-like key means do one of the helike ions */
691  else if( p.nMatch("HE-L" ) )
692  {
693  Param[nAsserts][0] = ipHE_LIKE;
694  chAssertLineLabel[nAsserts] = "dHel" ;
695  /* remember this is departure coefficient for some element */
696  strcpy( chAssertType[nAsserts] , "DI" );
697  /* now get element number for h ion from element name on card */
698  if( (wavelength[nAsserts] = (realnum)p.GetElem()) < 0 )
699  {
700  fprintf( ioQQQ,
701  " I could not identify an element name on this line.\n");
702  fprintf( ioQQQ, " Sorry.\n" );
703  cdEXIT(EXIT_FAILURE);
704  }
705  if( p.nMatch("ZEROOK") )
706  Param[nAsserts][1] = 1.;
707  else
708  Param[nAsserts][1] = 0.;
709  }
710 
711  /* this is H- h minus */
712  else if( p.nMatch("HMIN" ) )
713  {
714  /* label */
715  chAssertLineLabel[nAsserts] = "d H-" ;
716  /* remember this is departure coefficient for H- */
717  strcpy( chAssertType[nAsserts] , "d-" );
718  /* the wavelength is meaningless */
719  wavelength[nAsserts] = -1;
720  }
721  else
722  {
723  fprintf( ioQQQ,
724  " There must be a second key: H-LIke, HMINus, or HE-Like followed by element.\n");
725  fprintf( ioQQQ, " Sorry.\n" );
726  cdEXIT(EXIT_FAILURE);
727  }
728 
729  /* last check for key "excited" - which means to skip the ground state */
730  if( p.nMatch("EXCI" ) )
731  {
732  /* this is lowest level - do not do 0 */
733  AssertQuantity2[nAsserts] = 1.;
734  }
735  else
736  {
737  /* do the ground state */
738  AssertQuantity2[nAsserts] = 0.;
739  }
740  }
741 
742  /* monitor some results from map */
743  else if( p.nMatch(" MAP" ) )
744  {
745 
746  /* must have heating or cooling, since will check one or the other */
747  /* check heating cooling results from map at some temperature */
748  if( p.nMatch("HEAT" ) )
749  {
750  strcpy( chAssertType[nAsserts] , "mh" );
751  }
752  else if( p.nMatch("COOL" ) )
753  {
754  strcpy( chAssertType[nAsserts] , "mc" );
755  }
756  else
757  {
758  fprintf( ioQQQ,
759  " There must be a second key, HEATing or COOLing.\n");
760  fprintf( ioQQQ, " Sorry.\n" );
761  cdEXIT(EXIT_FAILURE);
762  }
763 
764  /* now get temperature for AssertQuantity2 array*/
765  AssertQuantity2[nAsserts] = p.FFmtRead();
766  if( p.lgEOL() )
767  {
768  p.NoNumb("temperature");
769  }
770 
771  if( AssertQuantity2[nAsserts] <= 10. )
772  {
773  /* entered as log, but we will compare with linear */
774  AssertQuantity2[nAsserts] =
775  exp10(AssertQuantity2[nAsserts] );
776  }
777 
778  /* print the temperature in the wavelength column */
779  wavelength[nAsserts] = (realnum)AssertQuantity2[nAsserts];
780 
781  /* heating or cooling, both log, put into error */
782  AssertQuantity[nAsserts] = p.FFmtRead();
783  if( p.lgEOL() )
784  {
785  p.NoNumb("heating/cooling");
786  }
787 
788  /* AssertQuantity array will have heating or cooling */
789  AssertQuantity[nAsserts] = exp10( AssertQuantity[nAsserts]);
790 
791  /* optional error, default available (cannot do before loop since we
792  * do not know how many numbers are on line */
793  AssertError[nAsserts] = p.FFmtRead();
794  if( p.lgEOL() )
795  {
796  /* default error was set above */
797  AssertError[nAsserts] = ErrorDefault;
798  }
799 
800  /* entered as a log, so print as a log too */
801  lgQuantityLog[nAsserts] = true;
802  }
803 
804  /* monitor column density of something */
805  else if( p.nMatch("COLU" ) )
806  {
807  /* this is column density */
808  strcpy( chAssertType[nAsserts] , "cd" );
809 
810  /* this says to look for molecular column density, also could be ion stage */
811  wavelength[nAsserts] = 0;
812 
813  string chLabel;
814 
815  if ( p.GetQuote(chLabel) == 0 )
816  {
817  // cdColm not yet able to cope with longer species
818  ASSERT(chLabel.length() <= NCHLAB-1);
819  trimTrailingWhiteSpace( chLabel );
820  chAssertLineLabel[nAsserts] = chLabel;
821  for (long i=chLabel.length();i<NCHLAB-1;++i)
822  chAssertLineLabel[nAsserts] += ' ';
823  }
824  else if( p.nMatchErase(" H2 ") )
825  {
826  chAssertLineLabel[nAsserts] = "H2" ;
827  }
828  else if( p.nMatchErase("H3+ "))
829  {
830  chAssertLineLabel[nAsserts] = "H3+" ;
831  }
832  else if( p.nMatchErase("H2+ "))
833  {
834  chAssertLineLabel[nAsserts] = "H2+" ;
835  }
836  else if( p.nMatchErase(" H- "))
837  {
838  chAssertLineLabel[nAsserts] = "H-" ;
839  }
840  else if( p.nMatchErase("H2G "))
841  {
842  chAssertLineLabel[nAsserts] = "H2g" ;
843  }
844  else if( p.nMatchErase("H2* "))
845  {
846  chAssertLineLabel[nAsserts] = "H2*" ;
847  }
848  else if( p.nMatchErase("HEH+"))
849  {
850  chAssertLineLabel[nAsserts] = "HeH+" ;
851  }
852  else if( p.nMatchErase(" O2 "))
853  {
854  chAssertLineLabel[nAsserts] = "O2" ;
855  }
856  else if( p.nMatchErase("H2O "))
857  {
858  chAssertLineLabel[nAsserts] = "H2O" ;
859  }
860  else if( p.nMatchErase(" C2 "))
861  {
862  chAssertLineLabel[nAsserts] = "C2" ;
863  }
864  else if( p.nMatchErase(" C3 "))
865  {
866  chAssertLineLabel[nAsserts] = "C3" ;
867  }
868  else if( p.nMatch(" CO "))
869  {
870  chAssertLineLabel[nAsserts] = "CO" ;
871  }
872  else if( p.nMatch("SIO "))
873  {
874  chAssertLineLabel[nAsserts] = "SiO" ;
875  }
876  else if( p.nMatch(" OH "))
877  {
878  chAssertLineLabel[nAsserts] = "OH" ;
879  }
880  else if( p.nMatch(" CN ") )
881  {
882  chAssertLineLabel[nAsserts] = "CN";
883  }
884  else if( p.nMatch(" CH ") )
885  {
886  chAssertLineLabel[nAsserts] = "CH" ;
887  }
888  else if( p.nMatch(" CH+") )
889  {
890  chAssertLineLabel[nAsserts] = "CH+" ;
891  }
892  else
893  {
894  fprintf( ioQQQ,
895  " I could not identify H2, H3+, H2+, H2g, H2*, H2H+, CO, O2, SiO, OH, C2 or C3 or on this line.\n");
896  fprintf( ioQQQ, " Sorry.\n" );
897  cdEXIT(EXIT_FAILURE);
898  }
899 
900  if( p.nMatch( "LEVE" ) )
901  {
902  if (chAssertLineLabel[nAsserts] != "H2")
903  {
904  fprintf( ioQQQ, " LEVEL option only implemented for H2.\n");
905  fprintf( ioQQQ, " Sorry.\n" );
906  cdEXIT(EXIT_FAILURE);
907  }
908  /* this is option for level-specific column density,
909  * next two numbers must be v then J */
910  Param[nAsserts][0] = p.FFmtRead();
911  if( p.lgEOL() )
912  p.NoNumb("level v" );
913  Param[nAsserts][1] = p.FFmtRead();
914  if( p.lgEOL() )
915  p.NoNumb("level J" );
916  /* wavelength will be 10. * vib + rot */
917  wavelength[nAsserts] = (realnum)(100.*Param[nAsserts][0] + Param[nAsserts][1]);
918  }
919  else
920  {
921  /* these are flags saying not to do state specific column densities */
922  Param[nAsserts][0] = -1.;
923  Param[nAsserts][1] = -1.;
924  }
925 
926  /* i was set above for start of scan */
927  /* now get log of column density */
928  AssertQuantity[nAsserts] = p.FFmtRead();
929  if( p.lgEOL() )
930  {
931  p.NoNumb("column density");
932  }
933  /* entered as log, but we will compare with linear */
934  AssertQuantity[nAsserts] =
935  exp10(AssertQuantity[nAsserts] );
936 
937  /* optional error, default available (cannot do before loop since we
938  * do not know how many numbers are on line */
939  AssertError[nAsserts] = p.FFmtRead();
940  if( p.lgEOL() )
941  {
942  /* default error was set above */
943  AssertError[nAsserts] = ErrorDefault;
944  }
945  /* the keyword log is special for this case, since H2 and CO column densities can
946  * be so very unstable. look for work log, in which case the error is log not linear.
947  * main column is always a log */
948  if( p.nMatch( " LOG" ) )
949  {
950  AssertError[nAsserts] = exp10( AssertError[nAsserts] );
951  }
952 
953  /* entered as a log, so print as a log too although asserted quantity is now linear */
954  lgQuantityLog[nAsserts] = true;
955  }
956 
957  /* monitor rate H2 forms on grain surfaces */
958  else if( p.nMatch("GRAI") && p.nMatch(" H2 ") )
959  {
960  /* this flag will mean h2 form on grains */
961  strcpy( chAssertType[nAsserts] , "g2" );
962  /* a label */
963  chAssertLineLabel[nAsserts] = "R H2" ;
964  /* now get the first number on the line, which must be the 2 in H2 */
965  nelem = (long int)p.FFmtRead();
966  if( nelem!=2 )
967  {
968  fprintf( ioQQQ,
969  " I did not find a 2, as in H2, as the first number on this line.\n");
970  fprintf( ioQQQ,
971  " The rate should be the second number.\n");
972  fprintf( ioQQQ, " Sorry.\n" );
973  cdEXIT(EXIT_FAILURE);
974  }
975  /* now past the 2 in h2, get the real number */
976  AssertQuantity[nAsserts] = p.FFmtRead();
977  if( p.lgEOL() )
978  {
979  p.NoNumb("grain H2 formation");
980  }
981  /* if negative (almost certainly the case) then the log of the rate coefficient */
982  if( AssertQuantity[nAsserts] <=0. )
983  AssertQuantity[nAsserts] = exp10(AssertQuantity[nAsserts] );
984  /* will not use this */
985  wavelength[nAsserts] = 0;
986 
987  /* optional error, default available (cannot do before loop since we
988  * do not know how many numbers are on line */
989  AssertError[nAsserts] = p.FFmtRead();
990  if( p.lgEOL() )
991  {
992  /* default error was set above */
993  AssertError[nAsserts] = ErrorDefault;
994  /* want to print as a log since so small */
995  lgQuantityLog[nAsserts] = true;
996  }
997  }
998 
999  /* monitor grain potential */
1000  else if( p.nMatch( "GRAI" ) && p.nMatch( "POTE") )
1001  {
1002  /* this flag will mean grain potential */
1003  strcpy( chAssertType[nAsserts] , "gp" );
1004  /* a label */
1005  chAssertLineLabel[nAsserts] = "GPot" ;
1006  /* now get the first number on the line */
1007  /* grain bin number */
1008  wavelength[nAsserts] = (realnum)p.FFmtRead();
1009  /* the potential itself, in volt, always linear */
1010  AssertQuantity[nAsserts] = p.FFmtRead();
1011 
1012  if( p.lgEOL() )
1013  {
1014  p.NoNumb("grain potential");
1015  }
1016 
1017  /* optional error, default available (cannot do before loop since we
1018  * do not know how many numbers are on line */
1019  AssertError[nAsserts] = p.FFmtRead();
1020 
1021  if( p.lgEOL() )
1022  {
1023  /* default error was set above */
1024  AssertError[nAsserts] = ErrorDefault;
1025  }
1026  }
1027 
1028  /* monitor mean temperature, monitor temperature hydrogen 2 8000 */
1029  else if( p.nMatch("TEMP") )
1030  {
1031  /* say that this will be mean temperature, electron or grain */
1032 
1033  /* t will indicate temperature average over radius, T over volume -
1034  * check whether keyword radius or volume occurs,
1035  * default will be radius */
1036  if( p.nMatch("VOLU") )
1037  {
1038  strcpy( chAssertType[nAsserts] , "T ");
1039  }
1040  else
1041  {
1042  /* this is default case, Fraction over radius */
1043  strcpy( chAssertType[nAsserts] , "t ");
1044  }
1045 
1046  /* first look for keyword Grains, since label silicate may be on it,
1047  * and this would trigger the element search */
1048  if( p.nMatch("GRAI") )
1049  {
1050  chAssertLineLabel[nAsserts] = "GTem" ;
1051  /* this is to make sure that pointer to grain type is valid, we check
1052  * that it is less than this below */
1053  nelem = LONG_MAX-2;
1054  /* the first numerical arg on the temper grain command is the grain
1055  * type as defined in Hazy, counting from 1. When it is used to
1056  * to get mean temps below we will subtract 1 from this to get onto
1057  * the c scale. but must leave on physical scale here to pass sanity
1058  * checks that occur below */
1059  }
1060 
1061  /* face is temperature at illuminated face of cloud */
1062  else if( p.nMatch("FACE") )
1063  {
1064  chAssertLineLabel[nAsserts] = "face" ;
1065  /* not used so zero out */
1066  nelem = 0;
1067  }
1068 
1069  /* mean H2 temperature */
1070  else if( p.nMatch(" H2 ") )
1071  {
1072  chAssertLineLabel[nAsserts] = "H2" ;
1073  /* not used so zero out */
1074  nelem = 0;
1075  }
1076 
1077  else if( p.nMatch( "21CM" ) )
1078  {
1079  /* not used so zero out */
1080  nelem = 0;
1081 
1082  if( p.nMatch( "MEAN" ) )
1083  {
1084  chAssertLineLabel[nAsserts] = "21CM" ;
1085  }
1086  else if( p.nMatch( "SPIN" ) )
1087  {
1088  chAssertLineLabel[nAsserts] = "SPIN" ;
1089  }
1090  else if( p.nMatch( "OPTI" ) )
1091  {
1092  chAssertLineLabel[nAsserts] = "OPTI" ;
1093  }
1094  else
1095  {
1096  fprintf( ioQQQ,
1097  "No option given."
1098  " One of MEAN, SPIN, and OPTICAL may be used.\n" );
1099  fprintf( ioQQQ, " Sorry.\n" );
1100  cdEXIT(EXIT_FAILURE);
1101  }
1102  }
1103 
1104  /* look for element label within quotes,
1105  * will be null terminated */
1106  else if( (nelem = p.GetElem()) >= 0 )
1107  {
1108  /* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem])
1109  * into array that will be used to get ionization after calculation */
1110  chAssertLineLabel[nAsserts] = elementnames.chElementNameShort[nelem];
1111  trimTrailingWhiteSpace( chAssertLineLabel[nAsserts] );
1112  }
1113  else
1114  {
1115  fprintf( ioQQQ,
1116  " I could not identify an element name on this line.\n");
1117  fprintf( ioQQQ, " Sorry.\n" );
1118  cdEXIT(EXIT_FAILURE);
1119  }
1120 
1121  /* now get ionization stage (or grain type), which will be saved into wavelength */
1122  if( chAssertLineLabel[nAsserts] == "face" ||
1123  chAssertLineLabel[nAsserts] == "21CM" ||
1124  chAssertLineLabel[nAsserts] == "SPIN" ||
1125  chAssertLineLabel[nAsserts] == "OPTI" )
1126  {
1127  /* this option is temperature at illuminated face so no element */
1128  wavelength[nAsserts] = 0;
1129  }
1130  else if( chAssertLineLabel[nAsserts] == "H2" )
1131  {
1132  int j;
1133  /* this option is temperature of H2 so element is zero */
1134  wavelength[nAsserts] = 0;
1135  /* first find the 2 in H2 */
1136  j = (int)p.FFmtRead();
1137  if( j != 2 )
1138  TotalInsanity();
1139  }
1140  else
1141  {
1142  wavelength[nAsserts] = (realnum)p.FFmtRead();
1143  if( p.lgEOL() )
1144  {
1145  p.NoNumb("ionization stage");
1146  }
1147 
1148  /* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */
1149  if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 )
1150  {
1151  fprintf( ioQQQ,
1152  " The ionization stage is inappropriate for this element.\n");
1153  fprintf( ioQQQ, " Sorry.\n" );
1154  cdEXIT(EXIT_FAILURE);
1155  }
1156  }
1157 
1158  if( p.nMatch( "GRID" ) )
1159  {
1160  AssertQuantity[nAsserts] = AssertVector[optimize.nOptimiz];
1161  }
1162  else
1163  {
1164  /* now get temperature, log if number is <= 10, else linear */
1165  AssertQuantity[nAsserts] = p.FFmtRead();
1166  if( p.lgEOL() )
1167  p.NoNumb("temperature");
1168 
1169  /* 21 is picked up by previous read; repeat */
1170  if( chAssertLineLabel[nAsserts] == "21CM" )
1171  {
1172  AssertQuantity[nAsserts] = p.FFmtRead();
1173  if( p.lgEOL() )
1174  p.NoNumb("temperature");
1175  }
1176  }
1177 
1178  /* optional error, default available */
1179  AssertError[nAsserts] = p.FFmtRead();
1180  if( p.lgEOL() )
1181  AssertError[nAsserts] = ErrorDefault;
1182 
1183  /* now make sure we end up with positive linear temperature
1184  * number is log if <=10 unless linear keyword appears */
1185  if( AssertQuantity[nAsserts] <= 10. && !p.nMatch( "LINE" ) )
1186  {
1187  /* log since negative or 0 */
1188  AssertQuantity[nAsserts] =
1189  exp10(AssertQuantity[nAsserts] );
1190  /* entered as a log, so print as a log too */
1191  lgQuantityLog[nAsserts] = true;
1192  }
1193  }
1194 
1195  /* monitor log of helium hydrogen ionization correction factor */
1196  else if( p.nMatch("HHEI") )
1197  {
1198  /* this flag will mean H-He icf */
1199  strcpy( chAssertType[nAsserts] , "hh" );
1200  /* say that this is zone numbering */
1201  chAssertLineLabel[nAsserts] = "HHei" ;
1202 
1203  /* now get the ionization correction factor, it is always the linear
1204  * quantity itself, since can be positive or negative*/
1205  AssertQuantity[nAsserts] = p.FFmtRead();
1206  if( p.lgEOL() )
1207  {
1208  p.NoNumb("ionization correction factor");
1209  }
1210  /* will not use this */
1211  wavelength[nAsserts] = 0;
1212 
1213  /* optional error, default available (cannot do before loop since we
1214  * do not know how many numbers are on line */
1215  AssertError[nAsserts] = p.FFmtRead();
1216  if( p.lgEOL() )
1217  {
1218  /* default error was set above */
1219  AssertError[nAsserts] = ErrorDefault;
1220  }
1221  }
1222 
1223  /* this large H2 molecule */
1224  else if( p.nMatch(" H2 ") )
1225  {
1226  /* ortho to para ratio for last computed zone */
1227  if( p.nMatch("ORTH") )
1228  {
1229  /* this flag will mean ortho to para ratio */
1230  strcpy( chAssertType[nAsserts] , "or" );
1231  /* say that this is ortho to para density ratio */
1232  chAssertLineLabel[nAsserts] = "orth" ;
1233  }
1234  else
1235  {
1236  fprintf( ioQQQ,
1237  " I could not identify a second keyword on this line.\n");
1238  fprintf( ioQQQ, " Sorry.\n" );
1239  cdEXIT(EXIT_FAILURE);
1240  }
1241 
1242  /* now get the first number, which better be the 2 in H2 */
1243  n2 = (long)p.FFmtRead();
1244  if( p.lgEOL() || n2 != 2 )
1245  {
1246  p.NoNumb("the 2 in H2 ?!");
1247  }
1248  AssertQuantity[nAsserts] = p.FFmtRead();
1249  /* will not use this */
1250  wavelength[nAsserts] = 0;
1251 
1252  /* optional error, default available (cannot do before loop since we
1253  * do not know how many numbers are on line */
1254  AssertError[nAsserts] = p.FFmtRead();
1255  if( p.lgEOL() )
1256  {
1257  /* default error was set above */
1258  AssertError[nAsserts] = ErrorDefault;
1259  }
1260  }
1261 
1262  /* monitor we are running in MPI mode */
1263  else if( p.nMatch(" MPI") )
1264  {
1265  /* this flag will mean number of zones */
1266  strcpy( chAssertType[nAsserts] , "mp" );
1267  /* say that this is zone numbering */
1268  chAssertLineLabel[nAsserts] = "mpi" ;
1269 
1270  wavelength[nAsserts] = 0.;
1271  AssertQuantity[nAsserts] = 1.;
1272  AssertError[nAsserts] = ErrorDefault;
1273  }
1274 
1275  /* monitor number of zones */
1276  else if( p.nMatch("NZON") )
1277  {
1278  /* this flag will mean number of zones */
1279  strcpy( chAssertType[nAsserts] , "z " );
1280  /* say that this is zone numbering */
1281  chAssertLineLabel[nAsserts] = "zone" ;
1282 
1283  /* now get number of zones */
1284  wavelength[nAsserts] = 0.;
1285  AssertQuantity[nAsserts] = p.FFmtRead();
1286  if( p.lgEOL() )
1287  p.NoNumb("zone number");
1288 
1289  /* optional error */
1290  AssertError[nAsserts] = p.FFmtRead();
1291  if( p.lgEOL() )
1292  AssertError[nAsserts] = ErrorDefault;
1293  }
1294 
1295  /* monitor (probably upper limit to) error in pressure across model */
1296  else if( p.nMatch("PRES") && p.nMatch("ERRO") )
1297  {
1298  /* this flag indicates ratio of standard deviation to the mean pressure */
1299  strcpy( chAssertType[nAsserts] , "pr" );
1300  /* say that this is error in pressure */
1301  chAssertLineLabel[nAsserts] = "pres" ;
1302 
1303  /* now get the pressure error, which will be saved into wavelength
1304  * in nearly all cases this is limit to error */
1305  wavelength[nAsserts] = 0;
1306  AssertQuantity[nAsserts] = (double)p.FFmtRead();
1307  if( p.lgEOL() )
1308  {
1309  p.NoNumb("pressure error");
1310  }
1311  else if( AssertQuantity[nAsserts] <= 0.)
1312  {
1313  /* number <= 0 is log of error */
1314  AssertQuantity[nAsserts] = exp10(AssertQuantity[nAsserts]);
1315  }
1316 
1317  /* optional error, default available (cannot do before loop since we
1318  * do not know how many numbers are on line */
1319  AssertError[nAsserts] = p.FFmtRead();
1320  if( p.lgEOL() )
1321  {
1322  /* default error was set above */
1323  AssertError[nAsserts] = ErrorDefault;
1324  }
1325  }
1326 
1327  else if( p.nMatch("PRADMAX") )
1328  {
1329  /* monitor pradmax - max ratio of rad to gas pressure */
1330  /* this flag indicates ratio of rad to gas pressure */
1331  strcpy( chAssertType[nAsserts] , "RM" );
1332  /* say that this is error in pressure */
1333  chAssertLineLabel[nAsserts] = "Prad" ;
1334 
1335  /* now get the pressure error, which will be saved into wavelength
1336  * in nearly all cases this is limit to error */
1337  wavelength[nAsserts] = 0;
1338  AssertQuantity[nAsserts] = (double)p.FFmtRead();
1339  if( p.lgEOL() )
1340  {
1341  p.NoNumb("PRADMAX");
1342  }
1343  else if( AssertQuantity[nAsserts] <= 0.)
1344  {
1345  /* number <= 0 is log of error */
1346  AssertQuantity[nAsserts] = exp10(AssertQuantity[nAsserts]);
1347  }
1348 
1349  /* optional error, default available (cannot do before loop since we
1350  * do not know how many numbers are on line */
1351  AssertError[nAsserts] = p.FFmtRead();
1352  if( p.lgEOL() )
1353  {
1354  /* default error was set above */
1355  AssertError[nAsserts] = ErrorDefault;
1356  }
1357  }
1358 
1359  /* monitor secondary ionization rate, csupra */
1360  else if( p.nMatch("CSUP") )
1361  {
1362  /* this flag will mean secondary ionization, entered as log */
1363  strcpy( chAssertType[nAsserts] , "sc" );
1364  /* say that this is sec ioniz */
1365  chAssertLineLabel[nAsserts] = "sion" ;
1366 
1367  /* now get rate, saved into monitor quantity */
1368  AssertQuantity[nAsserts] = p.FFmtRead();
1369  if( p.lgEOL() )
1370  {
1371  p.NoNumb("secondary ionization rate");
1372  }
1373  /* entered as log, make linear */
1374  AssertQuantity[nAsserts] = exp10( AssertQuantity[nAsserts] );
1375 
1376  /* no wavelength */
1377  wavelength[nAsserts] = 0;
1378 
1379  /* optional error, default available (cannot do before loop since we
1380  * do not know how many numbers are on line */
1381  AssertError[nAsserts] = p.FFmtRead();
1382  if( p.lgEOL() )
1383  {
1384  /* default error was set above */
1385  AssertError[nAsserts] = ErrorDefault;
1386  }
1387 
1388  /* we want to print the log of eden, not linear value */
1389  lgQuantityLog[nAsserts] = true;
1390  }
1391 
1392  /* monitor heating rate, erg/cm3/s, htot */
1393  else if( p.nMatch("HTOT") )
1394  {
1395  /* this flag will mean heating, entered as log */
1396  strcpy( chAssertType[nAsserts] , "ht" );
1397 
1398  /* say that this is heating rate */
1399  chAssertLineLabel[nAsserts] = "htot" ;
1400 
1401  /* now get rate, saved into monitor quantity */
1402  AssertQuantity[nAsserts] = p.FFmtRead();
1403  if( p.lgEOL() )
1404  {
1405  p.NoNumb("heating rate");
1406  }
1407  /* entered as log, make linear */
1408  AssertQuantity[nAsserts] = exp10( AssertQuantity[nAsserts] );
1409 
1410  /* no wavelength */
1411  wavelength[nAsserts] = 0;
1412 
1413  /* optional error, default available (cannot do before loop since we
1414  * do not know how many numbers are on line */
1415  AssertError[nAsserts] = p.FFmtRead();
1416  if( p.lgEOL() )
1417  {
1418  /* default error was set above */
1419  AssertError[nAsserts] = ErrorDefault;
1420  }
1421 
1422  /* we want to print the log of the heating, not linear value */
1423  lgQuantityLog[nAsserts] = true;
1424  }
1425 
1426  /* monitor cooling rate, erg/cm3/s, ctot */
1427  else if( p.nMatch("CTOT") )
1428  {
1429  /* this flag will mean cooling */
1430  strcpy( chAssertType[nAsserts] , "ct" );
1431 
1432  /* say that this is cooling rate */
1433  chAssertLineLabel[nAsserts] = "ctot";
1434 
1435  /* Look for GRID command */
1436  if( p.nMatch( "GRID" ) )
1437  {
1438  AssertQuantity[nAsserts] = AssertVector[optimize.nOptimiz];
1439  }
1440  else
1441  {
1442  AssertQuantity[nAsserts] = p.FFmtRead();
1443  /* now get rate, saved into monitor quantity */
1444  if( p.lgEOL() )
1445  {
1446  p.NoNumb("cooling rate");
1447  }
1448  }
1449  /* entered as log, make linear */
1450  AssertQuantity[nAsserts] = exp10( AssertQuantity[nAsserts] );
1451 
1452  /* no wavelength */
1453  wavelength[nAsserts] = 0;
1454 
1455  /* optional error, default available (cannot do before loop since we
1456  * do not know how many numbers are on line */
1457  AssertError[nAsserts] = p.FFmtRead();
1458  if( p.lgEOL() )
1459  {
1460  /* default error was set above */
1461  AssertError[nAsserts] = ErrorDefault;
1462  }
1463 
1464  /* we want to print the log of the heating, not linear value */
1465  lgQuantityLog[nAsserts] = true;
1466  }
1467 
1468  /* monitor number of iterations per zone, a test of convergence */
1469  else if( p.nMatch("ITRZ") )
1470  {
1471  /* this flag will mean number of iterations per zone */
1472  strcpy( chAssertType[nAsserts] , "iz" );
1473  /* say that this is iterations per zone */
1474  chAssertLineLabel[nAsserts] = "itrz" ;
1475 
1476  /* now get quantity */
1477  AssertQuantity[nAsserts] = p.FFmtRead();
1478  if( p.lgEOL() )
1479  {
1480  p.NoNumb("iterations per zone");
1481  }
1482  /* wavelength is meaningless */
1483  wavelength[nAsserts] = 0;
1484 
1485  /* optional error, default available */
1486  AssertError[nAsserts] = p.FFmtRead();
1487  if( p.lgEOL() )
1488  AssertError[nAsserts] = ErrorDefaultPerformance;
1489  }
1490 
1491  /* monitor electron density of the last zone */
1492  else if( p.nMatch("EDEN") )
1493  {
1494  /* this flag will mean electron density of the last zone */
1495  strcpy( chAssertType[nAsserts] , "e " );
1496  /* say that this is electron density */
1497  chAssertLineLabel[nAsserts] = "eden" ;
1498 
1499  /* now get electron density, which is a log */
1500  AssertQuantity[nAsserts] =
1501  exp10( p.FFmtRead() );
1502  if( p.lgEOL() )
1503  {
1504  p.NoNumb(" electron density of the last zone");
1505  }
1506 
1507  /* optional error, default available (cannot do before loop since we
1508  * do not know how many numbers are on line */
1509  AssertError[nAsserts] = p.FFmtRead();
1510  if( p.lgEOL() )
1511  {
1512  /* default error was set above */
1513  AssertError[nAsserts] = ErrorDefault;
1514  }
1515  wavelength[nAsserts] = 0;
1516 
1517  /* we want to print the log of eden, not linear value */
1518  lgQuantityLog[nAsserts] = true;
1519  }
1520 
1521  /* monitor energy density - Tu - of last zone */
1522  else if( p.nMatch(" TU ") )
1523  {
1524  /* this flag will mean energy density of the last zone */
1525  strcpy( chAssertType[nAsserts] , "Tu" );
1526  chAssertLineLabel[nAsserts] = "Tu" ;
1527 
1528  /* get temperature */
1529  AssertQuantity[nAsserts] = p.FFmtRead();
1530  if( p.lgEOL() )
1531  p.NoNumb("energy density of last zone");
1532 
1533  /* now make sure we end up with positive linear temperature
1534  * number is log if <=10 unless linear keyword appears */
1535  lgQuantityLog[nAsserts] = false;
1536  if( AssertQuantity[nAsserts] <= 10. && !p.nMatch( "LINE" ) )
1537  {
1538  /* log since negative or 0 */
1539  AssertQuantity[nAsserts] =
1540  exp10(AssertQuantity[nAsserts] );
1541  /* entered as a log, so print as a log too */
1542  }
1543 
1544  /* optional error, default available (cannot do before loop since we
1545  * do not know how many numbers are on line */
1546  AssertError[nAsserts] = p.FFmtRead();
1547  if( p.lgEOL() )
1548  /* default error was set above */
1549  AssertError[nAsserts] = ErrorDefault;
1550  wavelength[nAsserts] = 0;
1551  }
1552 
1553  /* monitor thickness or depth of model */
1554  else if( p.nMatch("THIC") || p.nMatch("DEPT") )
1555  {
1556  /* this flag will mean thickness or depth */
1557  strcpy( chAssertType[nAsserts] , "th" );
1558  /* say that this is thickness */
1559  chAssertLineLabel[nAsserts] = "thic" ;
1560 
1561  /* now get thickness, which is a log */
1562  AssertQuantity[nAsserts] = exp10( p.FFmtRead() );
1563  if( p.lgEOL() )
1564  {
1565  p.NoNumb("thickness or depth of model");
1566  }
1567 
1568  /* optional error, default available (cannot do before loop since we
1569  * do not know how many numbers are on line */
1570  AssertError[nAsserts] = p.FFmtRead();
1571  if( p.lgEOL() )
1572  {
1573  /* default error was set above */
1574  AssertError[nAsserts] = ErrorDefault;
1575  }
1576  wavelength[nAsserts] = 0;
1577 
1578  /* we want to print the log of eden, not linear value */
1579  lgQuantityLog[nAsserts] = true;
1580  }
1581 
1582  /* monitor outer radius of model */
1583  else if( p.nMatch("RADI") )
1584  {
1585  /* this flag will mean radius */
1586  strcpy( chAssertType[nAsserts] , "ra" );
1587  /* say that this is radius */
1588  chAssertLineLabel[nAsserts] = "radi" ;
1589 
1590  /* now get radius, which is a log */
1591  AssertQuantity[nAsserts] = exp10( p.FFmtRead() );
1592  if( p.lgEOL() )
1593  {
1594  p.NoNumb("outer radius");
1595  }
1596 
1597  /* optional error, default available (cannot do before loop since we
1598  * do not know how many numbers are on line */
1599  AssertError[nAsserts] = p.FFmtRead();
1600  if( p.lgEOL() )
1601  {
1602  /* default error was set above */
1603  AssertError[nAsserts] = ErrorDefault;
1604  }
1605  wavelength[nAsserts] = 0;
1606 
1607  /* we want to print the log of radius, not linear value */
1608  lgQuantityLog[nAsserts] = true;
1609  }
1610 
1611  /* monitor number of iterations */
1612  else if( p.nMatch("NITE") )
1613  {
1614  /* this flag will mean number of iterations */
1615  strcpy( chAssertType[nAsserts] , "Z " );
1616  /* say that this is iteration numbering */
1617  chAssertLineLabel[nAsserts] = "iter" ;
1618 
1619  wavelength[nAsserts] = 0.f;
1620 
1621  /* now get number of iterations, which will be saved into wavelength */
1622  AssertQuantity[nAsserts] = p.FFmtRead();
1623  if( p.lgEOL() )
1624  {
1625  p.NoNumb("number of iterations");
1626  }
1627 
1628  /* optional error, default available (cannot do before loop since we
1629  * do not know how many numbers are on line */
1630  AssertError[nAsserts] = p.FFmtRead();
1631  if( p.lgEOL() )
1632  {
1633  /* default error was set above */
1634  AssertError[nAsserts] = ErrorDefault;
1635  }
1636  }
1637 
1638  /* monitor terminal velocity, at end of calculation
1639  * input in km/s and saved that way, even though code uses cm/s */
1640  else if( p.nMatch("VELO") )
1641  {
1642  /* this flag will mean velocity */
1643  strcpy( chAssertType[nAsserts] , "Ve" );
1644  /* say that this is velocity */
1645  chAssertLineLabel[nAsserts] = "vel" ;
1646  wavelength[nAsserts] = 0;
1647 
1648  /* now get velocity */
1649  AssertQuantity[nAsserts] = p.FFmtRead();
1650  if( p.lgEOL() )
1651  p.NoNumb("terminal velocity");
1652 
1653  /* optional error, default available (cannot do before loop since we
1654  * do not know how many numbers are on line */
1655  AssertError[nAsserts] = p.FFmtRead();
1656  if( p.lgEOL() )
1657  {
1658  /* default error was set above */
1659  AssertError[nAsserts] = ErrorDefault;
1660  }
1661  }
1662  /* monitor nothing - a pacifier */
1663  else if( p.nMatch("NOTH") )
1664  {
1665  strcpy( chAssertType[nAsserts] , "NO" );
1666  chAssertLineLabel[nAsserts] = "noth" ;
1667  wavelength[nAsserts] = 0;
1668  AssertQuantity[nAsserts] = 0.;
1669  AssertError[nAsserts] = ErrorDefault;
1670  }
1671  else
1672  {
1673  /* did not recognize a command */
1674  fprintf( ioQQQ,
1675  " Unrecognized command. The line image was\n");
1676  p.PrintLine(ioQQQ);
1677  fprintf( ioQQQ,
1678  " The options I know about are: ionization, line, departure coefficient, map, column, "
1679  "temperature, nzone, csupre, htot, itrz, eden, thickness, niter, \n");
1680  fprintf( ioQQQ, " Sorry.\n" );
1681  cdEXIT(EXIT_FAILURE);
1682  }
1683 
1684  /* increment number of asserts and confirm that limit not exceeded */
1685  ++nAsserts;
1686  if( nAsserts >= NASSERTS )
1687  {
1688  fprintf(ioQQQ,
1689  " ParseMonitorResults: too many asserts, limit is NASSERT=%d\n",
1690  NASSERTS );
1691  cdEXIT(EXIT_FAILURE);
1692  }
1693  return;
1694 }
1695 
1696 inline double get_error_ratio ( double pred, double assert )
1697 {
1698  return 1. - safe_div(pred,assert,1.);
1699 }
1700 
1701 /*============================================================================*/
1702 /*lgCheckMonitors, checks asserts, last thing cloudy calls, returns true if all are ok, false if problems */
1703 bool lgCheckMonitors(
1704  /* this is the file we will write this to, usually standard output,
1705  * but can be save */
1706  FILE * ioMONITOR )
1707 {
1708  double PredQuan[NASSERTS] , RelError[NASSERTS];
1709  /* this will be true if the quantity was found, and false if not. Used to prevent
1710  * big botch flag when quantity not found (as in removed old he atom) */
1711  bool lgFound[NASSERTS];
1712  double relint , absint;
1713  bool lg1OK[NASSERTS];
1714  long i,j;
1715  /* This structure is for reporting another close match for asserts of line
1716  * intensities only. The zeroth, first, and second elements for each monitor are,
1717  * respectively, the first, second, and third matches the code finds, if any.
1718  * A negative number means no match. A positive number indicates the pointer
1719  * in the line stack of that match. */
1720  /* use multi_arr here to prevent bogus array bounds violations being reported by pgCC */
1721  multi_arr<long,2> ipDisambiguate(NASSERTS,3);
1722  long lgDisambiguate = false;
1723  char chCaps[NCHLAB], chFind[NCHLAB];
1724  realnum errorwave;
1725 
1726  DEBUG_ENTRY( "lgCheckMonitors()" );
1727 
1728  /* this is a global variable in monitor_results.h, and can be checked by
1729  * other routines to see if asserts are ok - (most runs will not use asserts) */
1730  lgMonitorsOK = true;
1731 
1732  /* will be used if there were big botched monitors */
1733  lgBigBotch = false;
1734 
1735  /* the optimize*.in and stars_oppim*.in tests in the test suite include
1736  * asserts while optimizing. We do not want to test the asserts during
1737  * the optimization process, since we will find blown asserts and report
1738  * major problems. We do want to test asserts on the final model however.
1739  * Printout will usually be turned off in all except the final model,
1740  * so do not to the monitor tests if we are optimizing but not printing */
1741  if( !called.lgTalk && optimize.lgOptimize )
1742  {
1743  /* just return */
1744  return true;
1745  }
1746 
1747  /*fprintf(ioQQQ , "DEBUG grid %li\n", optimize.nOptimiz );*/
1748 
1749  /* this will usually just return, but with table lines will check
1750  * existence of some lines */
1751  if( lines_table() )
1752  {
1753  lgBigBotch = true;
1754  lgMonitorsOK = false;
1755  }
1756 
1757  /* first check all asserts, there probably are none */
1758  for(i=0; i<nAsserts; ++i )
1759  {
1760  lg1OK[i] = true;
1761  PredQuan[i] = 0.;
1762  RelError[i] = 0.;
1763  for(j=0; j<3; ++j )
1764  ipDisambiguate[i][j] = -1;
1765 
1766  /* this flag is set false if we don't find the requested quantity */
1767  lgFound[i] = true;
1768 
1769  /* which type of monitor? */
1770  /* is it intensity? */
1771  if( strcmp( chAssertType[i] , "Lr")==0 )
1772  {
1773  /* this is an intensity, get the line, returns false if could not find it */
1774  ipDisambiguate[i][0] = cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint, iLineType[i] );
1775  if( ipDisambiguate[i][0] <= 0 )
1776  {
1777  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1778  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1779 
1780  fprintf( ioMONITOR,
1781  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1782  /* go to next line */
1783  lg1OK[i] = false;
1784  RelError[i] = 100000.;
1785  PredQuan[i] = 0;
1786  lg1OK[i] = false;
1787  lgMonitorsOK = false;
1788  lgFound[i] = false;
1789  continue;
1790  }
1791  else
1792  {
1793  /********* LINE DISAMBIGUATION *************/
1794  /* Here we look for lines with same label and small wavelength
1795  * differences so that we can disambiguate below */
1796 
1797  /* change chLabel to all caps */
1798  strcpy(chFind, chAssertLineLabel[i].c_str());
1799  caps(chFind);
1800 
1801  /* get the error associated with specified significant figures */
1802  errorwave = WavlenErrorGet( wavelength[i], LineSave.sig_figs );
1803 
1804  /* go through rest of line stack to look for close matches */
1805  for( j=1; j < LineSave.nsum; j++ )
1806  {
1807  /* don't bother with this one, we've already identified it. */
1808  if( j==ipDisambiguate[i][0] )
1809  continue;
1810 
1811  /* change chLabel to all caps to be like input chALab */
1812  cap4(chCaps, LineSave.lines[j].chALab());
1813 
1814  /* look for wavelengths within 3 error bars.
1815  * For example, for a line entered in Angstroms with
1816  * four significant figures, the error bar is 0.5 Ang.
1817  * So here we will find any lines within 1.5 Angstroms
1818  * of the
1819  * asserted wavelength. check wavelength and chLabel for a match */
1820  if( fabs(LineSave.lines[j].wavelength()-wavelength[i]) < 3.f*errorwave )
1821  {
1822  /* now see if labels agree */
1823  if( strcmp(chCaps,chFind) == 0 )
1824  {
1825  double relint1, absint1, current_error;
1826 
1827  cdLine_ip( j, &relint1, &absint1, iLineType[i] );
1828 
1829  if (AssertError[i] > 0.)
1830  current_error = fabs(1. - relint1/AssertQuantity[i]);
1831  else
1832  current_error = relint1 - AssertQuantity[i];
1833 
1834  if( current_error < 2.*fabs(AssertError[i]) ||
1835  current_error < 2.*fabs(RelError[i]) )
1836  {
1837  lgDisambiguate = true;
1838  /* if second match (element 1) is already set,
1839  * this is third match (element 2). Set and break out. */
1840  if( ipDisambiguate[i][1] > 0 )
1841  {
1842  ipDisambiguate[i][2] = j;
1843  break;
1844  }
1845  else
1846  {
1847  ipDisambiguate[i][1] = j;
1848  }
1849  }
1850  }
1851  }
1852  }
1853  }
1854 
1855  PredQuan[i] = relint;
1856  if (AssertError[i] > 0.0)
1857  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1858  else
1859  RelError[i] = PredQuan[i]-AssertQuantity[i] ;
1860  }
1861 
1862  else if( strcmp( chAssertType[i], "Lb" ) == 0 )
1863  {
1864  if( cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint, iLineType[i] ) <= 0 )
1865  {
1866  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1867  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1868 
1869  fprintf( ioMONITOR,
1870  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1871  /* go to next line */
1872  RelError[i] = 10000000.;
1873  PredQuan[i] = 0;
1874  lg1OK[i] = false;
1875  lgFound[i] = false;
1876  lgMonitorsOK = false;
1877  continue;
1878  }
1879 
1880  double relint_cb = 0.,
1881  absint_cb = 0.;
1882  if( cdLine( "Ca B", wavelength[i], &relint_cb, &absint_cb, iLineType[i] ) <= 0 )
1883  {
1884  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1885  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1886 
1887  fprintf( ioMONITOR,
1888  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1889  /* go to next line */
1890  RelError[i] = 10000000.;
1891  PredQuan[i] = 0;
1892  lg1OK[i] = false;
1893  lgFound[i] = false;
1894  lgMonitorsOK = false;
1895  continue;
1896  }
1897 
1898  PredQuan[i] = absint / absint_cb;
1899  if (AssertError[i] > 0.0)
1900  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1901  else
1902  RelError[i] = PredQuan[i] - AssertQuantity[i];
1903  // printf("i=%ld\t %s\t %g\t pred = %g\t asse = %g\t err = %g\n",
1904  }
1905 
1906  /*this is line luminosity */
1907  else if( strcmp( chAssertType[i] , "Ll")==0 )
1908  {
1909 // long indice=cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint );
1910 // fprintf(ioQQQ,"index line %ld %s %g \n", indice, LineSave.lines[indice].chALab(), LineSave.lines[indice].wavelength());
1911  /* this is a luminosity, get the line, returns false if could not find it */
1912  if( cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint, iLineType[i] ) <= 0 )
1913  //if (indice <= 0 )
1914  {
1915  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1916  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1917 
1918  fprintf( ioMONITOR,
1919  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1920  /* go to next line */
1921  RelError[i] = 10000000.;
1922  PredQuan[i] = 0;
1923  lg1OK[i] = false;
1924  lgFound[i] = false;
1925  lgMonitorsOK = false;
1926  continue;
1927  }
1928  PredQuan[i] = absint;
1929  if (AssertError[i] > 0.0)
1930  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1931  else
1932  RelError[i] = PredQuan[i] - AssertQuantity[i];
1933  // printf("i=%ld\t %s\t %g\t pred = %g\t asse = %g\t err = %g\n",
1934  // i, chAssertLineLabel[i].c_str(), wavelength[i], PredQuan[i], AssertQuantity[i], RelError[i] );
1935  }
1936  else if( strcmp( chAssertType[i] , "hh" ) == 0)
1937  {
1938  double hfrac , hefrac;
1939  /* get H ionization fraction, returns false if could not find it */
1940  if( cdIonFrac(
1941  /* four char string, null terminated, giving the element name */
1942  "hydr",
1943  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
1944  1,
1945  /* will be fractional ionization */
1946  &hfrac,
1947  /* how to weight the average, must be "VOLUME" or "RADIUS" */
1948  "VOLUME" ,
1949  /* do not want extra factor of density */
1950  false) )
1951  {
1952  fprintf( ioMONITOR,
1953  " monitor error: lgCheckMonitors could not find h ionization fraction \n");
1954  lg1OK[i] = false;
1955  RelError[i] = 0;
1956  PredQuan[i] = 0;
1957  lgFound[i] = false;
1958  continue;
1959  }
1960  if( cdIonFrac(
1961  /* four char string, null terminated, giving the element name */
1962  "heli",
1963  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
1964  1,
1965  /* will be fractional ionization */
1966  &hefrac,
1967  /* how to weight the average, must be "VOLUME" or "RADIUS" */
1968  "VOLUME" ,
1969  /* do not want extra factor of density */
1970  false) )
1971  {
1972  fprintf( ioMONITOR,
1973  " monitor error: lgCheckMonitors could not find h ionization fraction \n");
1974  lg1OK[i] = false;
1975  RelError[i] = 0;
1976  PredQuan[i] = 0;
1977  lgFound[i] = false;
1978  continue;
1979  }
1980  /* the helium hydrogen ionization correction factor */
1981  PredQuan[i] = hefrac-hfrac;
1982  /* two icf's in difference, no need to div by mean since already on scale with unity */
1983  RelError[i] = fabs(AssertQuantity[i] - (hefrac-hfrac) );
1984  }
1985 
1986  else if( strcmp( chAssertType[i] , "mp" ) == 0)
1987  {
1988  PredQuan[i] = cpu.i().lgMPI() ? 1. : 0.;
1989  /* use absolute error */
1990  RelError[i] = AssertQuantity[i] - PredQuan[i];
1991  }
1992 
1993  else if( strcmp( chAssertType[i] , "z " ) == 0)
1994  {
1995  /* this is the number of zones */
1996  PredQuan[i] = (double)nzone;
1997  if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
1998  RelError[i] = ForcePass(chAssertLimit[i]);
1999  else
2000  {
2001  /* two integers in difference */
2002  if (AssertError[i] > 0.)
2003  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2004  else
2005  RelError[i] = PredQuan[i] - AssertQuantity[i];
2006  }
2007  }
2008 
2009  else if( strcmp( chAssertType[i] , "or" ) == 0)
2010  {
2011  /* ortho to para ratio for large H2 molecule in last zone */
2012  PredQuan[i] = h2.ortho_density / SDIV( h2.para_density );
2013 
2014  /* this is relative error */
2015  if (AssertError[i] > 0.)
2016  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2017  else
2018  RelError[i] = PredQuan[i] - AssertQuantity[i];
2019  }
2020 
2021  else if( strcmp( chAssertType[i] , "g2" ) == 0)
2022  {
2023  /* check Jura rate, rate per vol that H2 forms on grain surfaces */
2024  PredQuan[i] = gv.rate_h2_form_grains_used_total;
2025  /* this is relative error */
2026  if (AssertError[i] > 0.)
2027  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2028  else
2029  RelError[i] = PredQuan[i] - AssertQuantity[i];
2030  }
2031 
2032  else if( strcmp( chAssertType[i] , "RM" ) == 0)
2033  {
2034  /* check Jura rate, rate per vol that H2 forms on grain surfaces */
2035  PredQuan[i] = pressure.RadBetaMax;
2036  /* this is relative error */
2037  if (AssertError[i] > 0.)
2038  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2039  else
2040  RelError[i] = PredQuan[i] - AssertQuantity[i];
2041  }
2042 
2043  else if( strcmp( chAssertType[i] , "pr" ) == 0)
2044  {
2045  /* standard deviation of the pressure */
2046  double sumx=0., sumx2=0., average;
2047  long int n;
2048  /* do sums to form standard deviation */
2049  for( n=0; n<nzone; n++ )
2050  {
2051  sumx += struc.pressure[n];
2052  sumx2 += POW2(struc.pressure[n]);
2053  }
2054  if( nzone>1 )
2055  {
2056  /* this is average */
2057  average = sumx/nzone;
2058  /* this is abs std */
2059  sumx = sqrt( (sumx2-POW2(sumx)/nzone)/(nzone-1) );
2060  /* save the relative std */
2061  PredQuan[i] = sumx / average;
2062  }
2063  else
2064  {
2065  PredQuan[i] = 0.;
2066  }
2067 
2068  // this is already relative error, do not need 1-ratio
2069  RelError[i] = PredQuan[i];
2070  }
2071 
2072  else if( strcmp( chAssertType[i] , "iz") == 0 )
2073  {
2074  /* this is number of iterations per zone, a test of convergence properties */
2075  if( nzone > 0 )
2076  PredQuan[i] = (double)(conv.nTotalIoniz-conv.nTotalIoniz_start)/(double)(nzone);
2077  else
2078  /* something big so monitor will botch. */
2079  PredQuan[i] = 1e10;
2080 
2081  if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
2082  RelError[i] = ForcePass(chAssertLimit[i]);
2083  else
2084  {
2085  /* this is relative error */
2086  if (AssertError[i] > 0.)
2087  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2088  else
2089  RelError[i] = PredQuan[i]- AssertQuantity[i];
2090  }
2091  }
2092 
2093  else if( strcmp( chAssertType[i] , "e ") == 0 )
2094  {
2095  /* this is electron density of the last zone */
2096  PredQuan[i] = dense.eden;
2097  /* this is relative error */
2098  if (AssertError[i] > 0.)
2099  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2100  else
2101  RelError[i] = PredQuan[i]- AssertQuantity[i];
2102  }
2103 
2104  else if( strcmp( chAssertType[i] , "Tu") == 0 )
2105  {
2106  /* this is radiation energy density of the last zone */
2107  PredQuan[i] = powpq((rfield.EnergyIncidCont+rfield.EnergyDiffCont)/
2108  (4.*STEFAN_BOLTZ),1,4);
2109  /* this is relative error */
2110  if (AssertError[i] > 0.)
2111  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2112  else
2113  RelError[i] = PredQuan[i]- AssertQuantity[i];
2114  }
2115 
2116  else if( strcmp( chAssertType[i] , "th") == 0 )
2117  {
2118  /* this is thickness */
2119  PredQuan[i] = radius.depth;
2120  /* this is relative error */
2121  if (AssertError[i] > 0.)
2122  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2123  else
2124  RelError[i] = PredQuan[i]- AssertQuantity[i];
2125  }
2126 
2127  else if( strcmp( chAssertType[i] , "ra") == 0 )
2128  {
2129  /* this is thickness */
2130  PredQuan[i] = radius.Radius;
2131  /* this is relative error */
2132  if (AssertError[i] > 0.)
2133  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2134  else
2135  RelError[i] = PredQuan[i]- AssertQuantity[i];
2136  }
2137 
2138  else if( strcmp( chAssertType[i] , "Ve") == 0 )
2139  {
2140  /* this is final velocity of wind in km/s (code uses cm/s) */
2141  PredQuan[i] = wind.windv/1e5;
2142  /* this is relative error */
2143  if (AssertError[i] > 0.)
2144  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2145  else
2146  RelError[i] = PredQuan[i]- AssertQuantity[i];
2147  }
2148 
2149  else if( strcmp( chAssertType[i] , "NO") == 0 )
2150  {
2151  /* this is nothing */
2152  PredQuan[i] = 0;
2153  /* this is relative error */
2154  RelError[i] = 0.;
2155  }
2156 
2157  else if( strcmp( chAssertType[i] , "sc") == 0 )
2158  {
2159  /* this is secondary ionization rate */
2160  PredQuan[i] = secondaries.csupra[ipHYDROGEN][0];
2161  /* this is relative error */
2162  if (AssertError[i] > 0.)
2163  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2164  else
2165  RelError[i] = PredQuan[i]- AssertQuantity[i];
2166  }
2167 
2168  else if( strcmp( chAssertType[i] , "ht") == 0 )
2169  {
2170  /* this is heating rate */
2171  PredQuan[i] = thermal.htot;
2172  /* this is relative error */
2173  if (AssertError[i] > 0.)
2174  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2175  else
2176  RelError[i] = PredQuan[i]- AssertQuantity[i];
2177  }
2178 
2179  else if( strcmp( chAssertType[i] , "ct") == 0 )
2180  {
2181  /* this is cooling rate */
2182  PredQuan[i] = thermal.ctot;
2183  /* this is relative error */
2184  if (AssertError[i] > 0.)
2185  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2186  else
2187  RelError[i] = PredQuan[i]- AssertQuantity[i];
2188  }
2189 
2190  else if( strcmp( chAssertType[i] , "Z ") == 0 )
2191  {
2192  /* this is the number of iterations */
2193  PredQuan[i] = (double)iteration;
2194  if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
2195  RelError[i] = ForcePass(chAssertLimit[i]);
2196  else
2197  {
2198  /* two integers in difference */
2199  if (AssertError[i] > 0.)
2200  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2201  else
2202  RelError[i] = PredQuan[i]- AssertQuantity[i];
2203  }
2204  }
2205 
2206  else if( strcmp( chAssertType[i] , "CB") == 0 )
2207  {
2208  long int nISOCaseB = (long)Param[i][0];
2209  long int nelemCaseB = (long)Param[i][1];
2210  string chElemLabelCaseB = chIonLbl( nelemCaseB+1, nelemCaseB+1-nISOCaseB );
2211 
2212  /* sets lowest quantum number index */
2213  int iCase;
2214  if( "Ca A" == chAssertLineLabel[i])
2215  iCase = 0;
2216  else if( "Ca B" == chAssertLineLabel[i] )
2217  iCase = 1;
2218  else
2219  TotalInsanity();
2220 
2221  iLineType[i] = 0;
2222 
2223  RelError[i] = 0.;
2224  long nHighestPrinted = iso_sp[nISOCaseB][nelemCaseB].n_HighestResolved_max;
2225  if( nISOCaseB == ipH_LIKE )
2226  {
2227  fprintf(ioMONITOR," Species nHi nLo Wl Computed Asserted error\n");
2228  /* limit of 10 is because that is all we printed and saved in prt_lines_hydro
2229  * wavelengths will come out of atmdat.WaveLengthCaseB - first index is
2230  * nelem on C scale, H is 0, second two are configurations of line on
2231  * physics scale, so Ha is 3-2 */
2232  for( long int ipLo=1+iCase; ipLo< MIN2(10,nHighestPrinted-1); ++ipLo )
2233  {
2234  for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi )
2235  {
2236  /* monitor the line */
2237  realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo];
2238  /* range option to restrict wavelength coverage */
2239  if( wl < Param[i][2] || wl > Param[i][3] )
2240  continue;
2241 
2242  double relint, absint, CBrelint, CBabsint;
2243  /* find the predicted line intensity */
2244  cdLine( chAssertLineLabel[i].c_str() , wl , &CBrelint , &CBabsint, iLineType[i] );
2245  if( CBrelint < Param[i][4] )
2246  continue;
2247  double error;
2248  /* now find the Case B intensity - may not all be present */
2249  if( cdLine( chElemLabelCaseB.c_str(), wl, &relint, &absint, iLineType[i] ) > 0 )
2250  {
2251  if (AssertError[i] > 0.)
2252  error = (CBabsint - absint)/MAX2(CBabsint , absint);
2253  else
2254  error = (CBabsint - absint);
2255  double RelativeError = fabs(error / AssertError[i]);
2256  /* start of line, flag problems */
2257  if( RelativeError < 1. )
2258  {
2259  if( RelativeError < 0.25 )
2260  {
2261  fprintf( ioMONITOR, " ChkMonitor ");
2262  }
2263  else if( RelativeError < 0.50 )
2264  {
2265  fprintf( ioMONITOR, " ChkMonitor - ");
2266  }
2267  else if( RelativeError < 0.75 )
2268  {
2269  fprintf( ioMONITOR, " ChkMonitor -- ");
2270  }
2271  else if( RelativeError < 0.90 )
2272  {
2273  fprintf( ioMONITOR, " ChkMonitor --- ");
2274  }
2275  else if( RelativeError < 0.95 )
2276  {
2277  fprintf( ioMONITOR, " ChkMonitor ---- ");
2278  }
2279  else if( RelativeError < 0.98 )
2280  {
2281  fprintf( ioMONITOR, " ChkMonitor ----- ");
2282  }
2283  else
2284  {
2285  fprintf( ioMONITOR, " ChkMonitor ------ ");
2286  }
2287 
2288  }
2289  else
2290  {
2291  fprintf( ioMONITOR, " ChkMonitor botch>>");
2292  }
2293  fprintf(ioMONITOR," %s %3li %3li ",
2294  chElemLabelCaseB.c_str() , ipHi , ipLo );
2295  prt_wl(ioMONITOR, wl );
2296  fprintf(ioMONITOR," %.2e %.2e %10.3f",
2297  log10(absint) , log10(CBabsint) , error );
2298  }
2299  else
2300  TotalInsanity();
2301  if( fabs(error) > fabs(AssertError[i]) )
2302  fprintf(ioMONITOR , " botch \n");
2303  else
2304  fprintf(ioMONITOR , "\n");
2305 
2306  PredQuan[i] = 0;
2307  AssertQuantity[i] = 0.;
2308  RelError[i] = MAX2( RelError[i] , fabs(error) );
2309 
2310  /* save sum which we will report later */
2311  MonitorResults.SumErrorCaseMonitor += RelError[i];
2312  ++MonitorResults.nSumErrorCaseMonitor;
2313 
2314  }
2315  }
2316  fprintf(ioMONITOR,"\n");
2317  }
2318  else if( nISOCaseB == ipHE_LIKE )
2319  {
2320  if( !dense.lgElmtOn[ipHELIUM] )
2321  {
2322  fprintf(ioQQQ,"PROBLEM monitor case B for a He is requested but He is not "
2323  "included.\n");
2324  fprintf(ioQQQ,"Do not turn off He if you want to monitor its spectrum.\n");
2325  cdEXIT(EXIT_FAILURE);
2326  }
2327 
2328  /* do He I as special case */
2329  fprintf(ioMONITOR," Wl Computed Asserted error\n");
2330  for( unsigned int ipLine=0; ipLine< atmdat.CaseBWlHeI.size(); ++ipLine )
2331  {
2332  /* monitor the line */
2333  realnum wl = atmdat.CaseBWlHeI[ipLine];
2334  /* range option to restrict wavelength coverage */
2335  if( wl < Param[i][2] || wl > Param[i][3] )
2336  continue;
2337  double relint , absint,CBrelint , CBabsint;
2338  cdLine( chAssertLineLabel[i].c_str(), wl, &CBrelint, &CBabsint, iLineType[i] );
2339  if( CBrelint < Param[i][4] )
2340  continue;
2341  double error;
2342  if( cdLine( chElemLabelCaseB.c_str(), wl, &relint, &absint, iLineType[i] ) > 0)
2343  {
2344  if (AssertError[i] > 0.0)
2345  error = (CBabsint - absint)/MAX2(CBabsint , absint);
2346  else
2347  error = (CBabsint - absint);
2348  double RelativeError = fabs(error / AssertError[i]);
2349  /* start of line, flag problems */
2350  if( RelativeError < 1. )
2351  {
2352  if( RelativeError < 0.25 )
2353  {
2354  fprintf( ioMONITOR, " ChkMonitor ");
2355  }
2356  else if( RelativeError < 0.50 )
2357  {
2358  fprintf( ioMONITOR, " ChkMonitor - ");
2359  }
2360  else if( RelativeError < 0.75 )
2361  {
2362  fprintf( ioMONITOR, " ChkMonitor -- ");
2363  }
2364  else if( RelativeError < 0.90 )
2365  {
2366  fprintf( ioMONITOR, " ChkMonitor --- ");
2367  }
2368  else if( RelativeError < 0.95 )
2369  {
2370  fprintf( ioMONITOR, " ChkMonitor ---- ");
2371  }
2372  else if( RelativeError < 0.98 )
2373  {
2374  fprintf( ioMONITOR, " ChkMonitor ----- ");
2375  }
2376  else
2377  {
2378  fprintf( ioMONITOR, " ChkMonitor ------ ");
2379  }
2380 
2381  }
2382  else
2383  {
2384  fprintf( ioMONITOR, " ChkMonitor botch>>");
2385  }
2386  prt_wl(ioMONITOR, wl );
2387  fprintf(ioMONITOR," %.2e %.2e %10.3f",
2388  absint , CBabsint , error );
2389  }
2390  else
2391  TotalInsanity();
2392  if( fabs(error) > fabs(AssertError[i]) )
2393  fprintf(ioMONITOR , " botch \n");
2394  else
2395  fprintf(ioMONITOR , "\n");
2396 
2397  PredQuan[i] = 0;
2398  AssertQuantity[i] = 0.;
2399  RelError[i] = MAX2( RelError[i] , fabs(error) );
2400 
2401  /* save sum which we will report later */
2402  MonitorResults.SumErrorCaseMonitor += RelError[i];
2403  ++MonitorResults.nSumErrorCaseMonitor;
2404  }
2405  fprintf(ioMONITOR,"\n");
2406  }
2407  else
2408  TotalInsanity();
2409  }
2410 
2411  // departure coefficients for a species given in quotes
2412  else if( strcmp( chAssertType[i] , "DC") == 0 )
2413  {
2414  string this_species = strAssertSpecies[i];
2415  if( this_species.find( '[' ) == string::npos )
2416  {
2417  // keyword EXCITED appeared
2418  if( AssertQuantity2[i] == 1 )
2419  this_species += "[2:]";
2420  else
2421  this_species += "[:]";
2422  }
2423 
2424  vector<long> speciesLevels;
2425  const molezone *sp = getLevelsGeneric( this_species.c_str(), false, speciesLevels );
2426  if( sp == null_molezone )
2427  {
2428  fprintf( ioQQQ, "PROBLEM Could not find species between quotes: \"%s\".\n",
2429  strAssertSpecies[i].c_str() );
2430  cdEXIT(EXIT_FAILURE);
2431  }
2432 
2433  if( sp->levels == NULL )
2434  {
2435  fprintf( ioQQQ, "WARNING Species '%s' has no internal structure."
2436  " Cannot compute departure coefficient\n",
2437  strAssertSpecies[i].c_str() );
2438  // Just report unity?
2439  PredQuan[i] = 1.;
2440  RelError[i] = 0.;
2441  }
2442  else if( speciesLevels.size() > 0 )
2443  {
2444  ASSERT( sp->levels->size() > 0 );
2445  RelError[i] = 0.;
2446  PredQuan[i] = 0.;
2447 
2448  long numPrintLevels = 0;
2449  for( vector<long>::const_iterator ilvl = speciesLevels.begin(); ilvl != speciesLevels.end(); ilvl++ )
2450  {
2451  const qStateConstProxy& st = (*sp->levels)[ *ilvl ];
2452  if( st.status() == LEVEL_INACTIVE )
2453  continue;
2454  ++numPrintLevels;
2455  PredQuan[i] += st.DepartCoef();
2456  }
2457  ASSERT( numPrintLevels > 0 );
2458  PredQuan[i] /= (double)(numPrintLevels);
2459  RelError[i] = AssertQuantity[i] - PredQuan[i];
2460 
2461 #if 0
2462  if( 0 && fp_equal( Param[i][1], 1. ) && PredQuan[i]==0. )
2463  {
2464  // this should only happen if either the present stage or the parent has zero density.
2465  ASSERT( dense.xIonDense[nelem][nelem-ipISO]==0. || dense.xIonDense[nelem][nelem+1-ipISO]==0. );
2466  PredQuan[i] = AssertQuantity[i];
2467  RelError[i] = 0.;
2468  }
2469 #endif
2470  }
2471  else
2472  {
2473  fprintf( ioQQQ, "WARNING Requested level(s) in '%s' do not exist\n",
2474  strAssertSpecies[i].c_str() );
2475  PredQuan[i] = 0.;
2476  RelError[i] = AssertQuantity[i];
2477  }
2478  }
2479 
2480  /* departure coefficients for something in isoelectronic sequences */
2481  else if( strcmp( chAssertType[i] , "DI") == 0 )
2482  {
2483  /* this is departure coefficient for XX-like ion given by wavelength */
2484  /* stored number was element number on C scale */
2485  long ipISO = (long)Param[i][0];
2486  long nelem = (long)wavelength[i];
2487  if( !dense.lgElmtOn[nelem] )
2488  {
2489  fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",nelem);
2490  PredQuan[i] = 0.;
2491  RelError[i] = 0.;
2492  }
2493  else
2494  {
2495  RelError[i] = 0.;
2496  PredQuan[i] = 0.;
2497  long numPrintLevels = iso_sp[ipISO][nelem].numLevels_local - (long)AssertQuantity2[i];
2498  for( long n=(long)AssertQuantity2[i]; n<numPrintLevels+(long)AssertQuantity2[i]; ++n )
2499  {
2500  PredQuan[i] += iso_sp[ipISO][nelem].st[n].DepartCoef();
2501  }
2502  ASSERT( numPrintLevels > 0 );
2503  PredQuan[i] /= (double)(numPrintLevels);
2504  RelError[i] = AssertQuantity[i] - PredQuan[i];
2505 
2506  if( fp_equal( Param[i][1], 1. ) && PredQuan[i]==0. )
2507  {
2508  // this should only happen if either the present stage or the parent has zero density.
2509  ASSERT( dense.xIonDense[nelem][nelem-ipISO]==0. || dense.xIonDense[nelem][nelem+1-ipISO]==0. );
2510  PredQuan[i] = AssertQuantity[i];
2511  RelError[i] = 0.;
2512  }
2513  }
2514  }
2515 
2516  /* this is H- departure coefficient */
2517  else if( strcmp( chAssertType[i] , "d-") == 0 )
2518  {
2519  PredQuan[i] = hmi.hmidep;
2520  RelError[i] = AssertQuantity[i] - hmi.hmidep;
2521  }
2522 
2523  /* this would be ionization fraction */
2524  else if( (strcmp( chAssertType[i] , "f ") == 0) ||
2525  (strcmp( chAssertType[i] , "F ") == 0) )
2526  {
2527  char chWeight[7];
2528  if( strcmp( chAssertType[i] , "F ") == 0 )
2529  {
2530  strcpy( chWeight , "VOLUME" );
2531  }
2532  else
2533  {
2534  /* this is default case, Fraction over radius */
2535  strcpy( chWeight , "RADIUS" );
2536  }
2537  /* get ionization fraction, returns false if could not find it */
2538  if( cdIonFrac(
2539  /* four char string, null terminated, giving the element name */
2540  chAssertLineLabel[i].c_str(),
2541  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
2542  (long)wavelength[i],
2543  /* will be fractional ionization */
2544  &relint,
2545  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2546  chWeight ,
2547  /* do not want extra factor of density */
2548  false) )
2549  {
2550  fprintf( ioMONITOR,
2551  " monitor error: lgCheckMonitors could not find a line with label %s %f \n",
2552  chAssertLineLabel[i].c_str() , wavelength[i] );
2553  /* go to next line */
2554  lg1OK[i] = false;
2555  RelError[i] = 0;
2556  PredQuan[i] = 0;
2557  lgFound[i] = false;
2558  continue;
2559  }
2560  /* this is ionization fraction */
2561  PredQuan[i] = relint;
2562  if (AssertError[i] > 0.)
2563  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2564  else
2565  RelError[i] = PredQuan[i]- AssertQuantity[i];
2566  }
2567 
2568  /* this would be column density of several molecules */
2569  else if( strcmp( chAssertType[i] , "cd") == 0)
2570  {
2571  /* for H2 column density - total or for a state? */
2572  if( ( chAssertLineLabel[i] == "H2" ) && (Param[i][0] >= 0.) )
2573  {
2574  /* this branch get state specific column density */
2575  /* get total H2 column density */
2576  if( (relint = cdH2_colden( (long)Param[i][0] , (long)Param[i][1] ) ) < 0. )
2577  {
2578  fprintf(ioQQQ," PROBLEM lgCheckMonitors did not find v=%li, J=%li for H2 column density.\n",
2579  (long)Param[i][0] , (long)Param[i][1] );
2580  lg1OK[i] = false;
2581  RelError[i] = 0;
2582  PredQuan[i] = 0;
2583  lgFound[i] = false;
2584  continue;
2585  }
2586  }
2587  else
2588  {
2589  /* get ionization fraction, returns 0 if all ok */
2590  if( cdColm(
2591  /* four char string, null terminated, giving the element name */
2592  chAssertLineLabel[i].c_str(),
2593  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2594  * zero for molecule*/
2595  (long)wavelength[i],
2596  /* will be fractional ionization */
2597  &relint) )
2598  {
2599  fprintf( ioMONITOR,
2600  " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
2601  chAssertLineLabel[i].c_str() , wavelength[i] );
2602  /* go to next line */
2603  lg1OK[i] = false;
2604  RelError[i] = 0;
2605  PredQuan[i] = 0;
2606  lgFound[i] = false;
2607  continue;
2608  }
2609  }
2610  /* this is ionization fraction */
2611  PredQuan[i] = relint;
2612  if (AssertError[i] > 0.)
2613  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2614  else
2615  RelError[i] = PredQuan[i]- AssertQuantity[i];
2616  }
2617 
2618  /* this would be molecular fraction of CO or H2 */
2619  else if( (strcmp( chAssertType[i] , "MF") == 0) || (strcmp( chAssertType[i] , "mf") == 0) )
2620  {
2621  /* get molecular fraction, returns 0 if all ok */
2622  relint = 0.;
2623  if( chAssertLineLabel[i] == "H2" )
2624  {
2625  /* get total H2 column density */
2626  if( cdColm("H2 " , 0,
2627  /* will be fractional ionization */
2628  &relint) )
2629  TotalInsanity();
2630 
2631  relint = relint / colden.colden[ipCOL_HTOT];
2632  }
2633  else
2634  {
2635  fprintf( ioMONITOR,
2636  " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
2637  chAssertLineLabel[i].c_str() , wavelength[i] );
2638  /* go to next line */
2639  lg1OK[i] = false;
2640  RelError[i] = 0;
2641  PredQuan[i] = 0;
2642  continue;
2643  }
2644  /* this is ionization fraction */
2645  PredQuan[i] = relint;
2646  if (AssertError[i] > 0.)
2647  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2648  else
2649  RelError[i] = PredQuan[i]- AssertQuantity[i];
2650  }
2651 
2652  /* check heating/cooling at some temperature in a thermal map */
2653  else if( strcmp( chAssertType[i] , "mh") == 0 || strcmp( chAssertType[i] , "mc") == 0)
2654  {
2655  /* check heating or cooling (stored in error array) at temperature in monitor results */
2656  /* check that map was done, and arrays have nmap elements */
2657  if( hcmap.nMapAlloc == 0 )
2658  {
2659  /* this happens if map not done and space for h/c not allocated */
2660  fprintf( ioMONITOR,
2661  " monitor error: lgCheckMonitors cannot check map since map not done.\n");
2662  /* go to next line */
2663  lg1OK[i] = false;
2664  RelError[i] = 0;
2665  PredQuan[i] = 0;
2666  continue;
2667  }
2668  /* now check that requested temperature is within the range of the map we computed */
2669  if( AssertQuantity2[i]< hcmap.temap[0] || AssertQuantity2[i]> hcmap.temap[hcmap.nmap-1] )
2670  {
2671  fprintf( ioMONITOR,
2672  " monitor error: lgCheckMonitors cannot check map since temperature not within range.\n");
2673  /* go to next line */
2674  lg1OK[i] = false;
2675  RelError[i] = 0;
2676  PredQuan[i] = 0;
2677  continue;
2678  }
2679 
2680  /* we should have valid data - find closest temperature >- requested temperature */
2681  j = 0;
2682  while( AssertQuantity2[i]>hcmap.temap[j]*1.001 && j < hcmap.nmap )
2683  {
2684  ++j;
2685  }
2686 
2687  /* j points to correct cell in heating cooling array */
2688  /* we will not interpolate, just use this value, and clobber te to prove it*/
2689  if( strcmp( chAssertType[i] , "mh") == 0 )
2690  {
2691  /* heating */
2692  PredQuan[i] = hcmap.hmap[j];
2693  chAssertLineLabel[i] = "MapH" ;
2694  }
2695  else if( strcmp( chAssertType[i] , "mc") == 0)
2696  {
2697  /* cooling */
2698  PredQuan[i] = hcmap.cmap[j];
2699  chAssertLineLabel[i] = "MapC" ;
2700  }
2701  if (AssertError[i] > 0.)
2702  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2703  else
2704  RelError[i] = PredQuan[i]- AssertQuantity[i];
2705  }
2706 
2707  /* this will be an average temperature */
2708  else if( (strcmp( chAssertType[i] , "t ") == 0) ||
2709  (strcmp( chAssertType[i] , "T ") == 0) )
2710  {
2711  char chWeight[7];
2712  if( strcmp( chAssertType[i] , "T ") == 0 )
2713  {
2714  strcpy( chWeight , "VOLUME" );
2715  }
2716  else
2717  {
2718  /* this is default case, Fraction over radius */
2719  strcpy( chWeight , "RADIUS" );
2720  }
2721 
2722  /* options are average Te for ion, temp at ill face, or temp for grain */
2723  if( chAssertLineLabel[i] == "GTem" )
2724  {
2725  size_t nd;
2726  /* the minus one is because the grain types are counted from one,
2727  * but stuffed into the c array, that counts from zero */
2728  nd = (size_t)wavelength[i]-1;
2729  if( nd >= gv.bin.size() )
2730  {
2731  fprintf( ioQQQ, "Illegal grain number found: %f\n" , wavelength[i] );
2732  fprintf( ioQQQ, "Use 1 for first grain that is turned on, " );
2733  fprintf( ioQQQ, "2 for second, etc....\n" );
2734  fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" );
2735  cdEXIT(EXIT_FAILURE);
2736  }
2737  relint = gv.bin[nd]->avdust/radius.depth_x_fillfac;
2738  }
2739  else if( chAssertLineLabel[i] == "face" )
2740  {
2741  /* this is the temperature at the illuminated face */
2742  relint = struc.testr[0];
2743  }
2744  else
2745  {
2746  long ion_stage = long( wavelength[i] );
2747  if( chAssertLineLabel[i] == "21CM" ||
2748  chAssertLineLabel[i] == "SPIN" ||
2749  chAssertLineLabel[i] == "OPTI" )
2750  ion_stage = 0;
2751 
2752  /* get temperature, returns false if could not find it */
2753  if( cdTemp(
2754  /* four char string, null terminated, giving the element name */
2755  chAssertLineLabel[i].c_str(),
2756  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
2757  ion_stage,
2758  /* will be mean temperatue */
2759  &relint,
2760  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2761  chWeight ) )
2762  {
2763  fprintf( ioMONITOR,
2764  " monitor error: lgCheckMonitors could not find an ion with label %s ion %li \n",
2765  chAssertLineLabel[i].c_str() , (long)wavelength[i] );
2766  /* go to next line */
2767  lg1OK[i] = false;
2768  RelError[i] = 0;
2769  PredQuan[i] = 0;
2770  lgFound[i] = false;
2771  continue;
2772  }
2773  }
2774  /* this is the temperature */
2775  PredQuan[i] = relint;
2776  if (AssertError[i] > 0.)
2777  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2778  else
2779  RelError[i] = PredQuan[i]- AssertQuantity[i];
2780  }
2781 
2782  /* this would be grain potential in volt */
2783  else if( strcmp( chAssertType[i], "gp" ) == 0 )
2784  {
2785  /* the minus one is because the grain types are counted from one,
2786  * but stuffed into the c array, that counts from zero */
2787  size_t nd = (size_t)wavelength[i]-1;
2788  if( nd >= gv.bin.size() )
2789  {
2790  fprintf( ioQQQ, "Illegal grain number found: %g\n" , wavelength[i] );
2791  fprintf( ioQQQ, "Use 1 for first grain that is turned on, " );
2792  fprintf( ioQQQ, "2 for second, etc....\n" );
2793  fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" );
2794  cdEXIT(EXIT_FAILURE);
2795  }
2796 
2797  /* get average grain potential in volt, always averaged over radius */
2798  PredQuan[i] = gv.bin[nd]->avdpot/radius.depth_x_fillfac;
2799  /* actually absolute error, potential can be zero! */
2800  RelError[i] = AssertQuantity[i] - PredQuan[i];
2801  }
2802 
2803  else
2804  {
2805  fprintf( ioMONITOR,
2806  " monitor error: lgCheckMonitors received an insane chAssertType=%s, impossible\n",
2807  chAssertType[i] );
2808  ShowMe();
2809  cdEXIT(EXIT_FAILURE);
2810  }
2811 
2812  if( chAssertLimit[i] == '=' )
2813  {
2814  /* predicted quantity should be within error of expected */
2815  if( fabs(RelError[i]) > fabs(AssertError[i]) )
2816  {
2817  lg1OK[i] = false;
2818  lgMonitorsOK = false;
2819  }
2820  }
2821  else if( chAssertLimit[i] == '<' )
2822  {
2823  /* expected is an upper limit, so PredQuan/AssertQuantity should
2824  * be less than one, and so RelError =1-PredQuan[i]/AssertQuantity[i]
2825  * should be >= 0
2826  * in case of iterations or zones, iter < iterations should
2827  * trigger botched monitor when iter == iterations */
2828  if( RelError[i] <= 0. )
2829  {
2830  lg1OK[i] = false;
2831  lgMonitorsOK = false;
2832  }
2833  }
2834  else if( chAssertLimit[i] == '>' )
2835  {
2836  /* expected is a lower limit, so PredQuan/AssertQuantity should
2837  * be greater than one, and so RelError should be negative */
2838  if( RelError[i] >= 0. )
2839  {
2840  lg1OK[i] = false;
2841  lgMonitorsOK = false;
2842  }
2843  }
2844  }
2845 
2846  /* only print summary if we are talking */
2847  if( called.lgTalk && nAsserts>0 )
2848  {
2849  time_t now;
2850 
2851  /* First disambiguate any line identifications */
2852  if( lgDisambiguate )
2853  {
2854  /* change significant figures of WL for this printout */
2855  long sigfigsav = LineSave.sig_figs;
2856  double relint1, relint2, absint1;
2857 
2858  LineSave.sig_figs = LineSave.sig_figs_max;
2859 
2860  fprintf( ioMONITOR, "=============Line Disambiguation=======================================================\n" );
2861  fprintf( ioMONITOR, " Wavelengths || Intensities \n" );
2862  fprintf( ioMONITOR, "Label line match1 match2 match3 || asserted match1 match2 match3\n" );
2863 
2864  for( i=0; i<nAsserts; ++i )
2865  {
2866  if( ipDisambiguate[i][1] > 0 )
2867  {
2868  fprintf( ioMONITOR , "%-*s ", NCHLAB-1, chAssertLineLabel[i].c_str() );
2869  prt_wl( ioMONITOR , wavelength[i] );
2870  fprintf( ioMONITOR , " " );
2871  prt_wl( ioMONITOR , LineSave.lines[ipDisambiguate[i][0]].wavelength() );
2872  fprintf( ioMONITOR , " " );
2873  prt_wl( ioMONITOR , LineSave.lines[ipDisambiguate[i][1]].wavelength() );
2874  fprintf( ioMONITOR , " " );
2875  if( ipDisambiguate[i][2] > 0 )
2876  {
2877  prt_wl( ioMONITOR , LineSave.lines[ipDisambiguate[i][2]].wavelength() );
2878  cdLine_ip( ipDisambiguate[i][2], &relint2, &absint1, iLineType[i] );
2879  }
2880  else
2881  {
2882  fprintf( ioMONITOR , "--------" );
2883  relint2 = 0.0;
2884  }
2885  fprintf( ioMONITOR , " ||" );
2886 
2887  cdLine_ip( ipDisambiguate[i][1], &relint1, &absint1, iLineType[i] );
2888 
2889  if( lgPrtSciNot )
2890  {
2891  fprintf( ioMONITOR , " %10.3e %10.3e %10.3e %10.3e\n",
2892  AssertQuantity[i],
2893  PredQuan[i] ,
2894  relint1,
2895  relint2 );
2896  }
2897  else
2898  {
2899  fprintf( ioMONITOR , " %10.4f %10.4f %10.4f %10.4f\n",
2900  AssertQuantity[i],
2901  PredQuan[i] ,
2902  relint1,
2903  relint2 );
2904  }
2905  }
2906  }
2907  fprintf( ioMONITOR, "\n" );
2908 
2909  /* revert to original significant figures */
2910  LineSave.sig_figs = sigfigsav;
2911  }
2912 
2913  /* write start of title and version number of code */
2914  fprintf( ioMONITOR, "=============Results of monitors: Cloudy %s ", t_version::Inst().chVersion );
2915 
2916  /* usually print date and time info - do not if "no times" command entered,
2917  * which set this flag false */
2918  if( prt.lgPrintTime )
2919  {
2920  /* now add date of this run */
2921  now = time(NULL);
2922  /* now print this time at the end of the string. the system put cr at the end of the string */
2923  fprintf(ioMONITOR,"%s", ctime(&now) );
2924  }
2925  else
2926  {
2927  fprintf(ioMONITOR,"\n" );
2928  }
2929 
2930  if( lgMonitorsOK )
2931  {
2932  fprintf( ioMONITOR, " No errors were found. Summary follows.\n");
2933  }
2934  else
2935  {
2936  fprintf( ioMONITOR, " Errors were found. Summary follows.\n");
2937  }
2938 
2939  fprintf( ioMONITOR,
2940  " %-*s%*s computed asserted Rel Err Set err type \n",
2941  NCHLAB, "Label", LineSave.wl_length, "line" );
2942  /* now print a summary */
2943  for( i=0; i<nAsserts; ++i )
2944  {
2945  PrtOneMonitor( ioMONITOR, chAssertType[i], chAssertLineLabel[i].c_str(),
2946  wavelength[i], iLineType[i], PredQuan[i], chAssertLimit[i], AssertQuantity[i],
2947  RelError[i], AssertError[i], lg1OK[i], lgQuantityLog[i], lgFound[i] );
2948  }
2949  fprintf( ioMONITOR , " \n");
2950 
2951  /* NB - in following, perl scripts detect these strings - be careful if they
2952  * are changed - the scripts may no longer detect major errors */
2953  if( !lgMonitorsOK && lgBigBotch )
2954  {
2955  /* there were big botches */
2956  fprintf( ioMONITOR, " BIG BOTCHED MONITORS!!! Big Botched Monitors!!! \n");
2957  }
2958  else if( !lgMonitorsOK )
2959  {
2960  fprintf( ioMONITOR, " BOTCHED MONITORS!!! Botched Monitors!!! \n");
2961  }
2962 
2963  if( MonitorResults.nSumErrorCaseMonitor>0 )
2964  {
2965  fprintf(ioMONITOR,"\n The mean of the %li monitor Case A, B relative "
2966  "residuals is %.2f\n\n" ,
2967  MonitorResults.nSumErrorCaseMonitor,
2968  MonitorResults.SumErrorCaseMonitor /MonitorResults.nSumErrorCaseMonitor );
2969  }
2970 
2971  /* explain how we were compiled, but only if printing time */
2972  if( prt.lgPrintTime )
2973  {
2974  fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo );
2975  }
2976  }
2977  return lgMonitorsOK;
2978 }
2979 
2980 STATIC void prtLineType( FILE *ioMONITOR, const int iLineType )
2981 {
2982  switch( iLineType )
2983  {
2984  case -1:
2985  /* The monitor is not an atomic transition. */
2986  fprintf( ioMONITOR, " " );
2987  break;
2988  case 0:
2989  fprintf( ioMONITOR, " intr " );
2990  break;
2991  case 1:
2992  fprintf( ioMONITOR, " emer " );
2993  break;
2994  case 2:
2995  fprintf( ioMONITOR, " intr cumu" );
2996  break;
2997  case 3:
2998  fprintf( ioMONITOR, " emer cumu" );
2999  break;
3000  default:
3001  fprintf( ioMONITOR, "ERROR: Unrecognized line type: %d\n",
3002  iLineType );
3003  cdEXIT( EXIT_FAILURE );
3004  break;
3005  }
3006 }
3007 
3008 void PrtOneMonitor( FILE *ioMONITOR, const char* chAssertType, const string chAssertLineLabel,
3009  const realnum wavelength, const int iLineType, const double PredQuan, const char chAssertLimit, const double AssertQuantity,
3010  const double RelError, const double AssertError, const bool lg1OK, const bool lgQuantityLog, const bool lgFound )
3011 {
3012  DEBUG_ENTRY( "PrtOneMonitor()" );
3013 
3014  double prtPredQuan, prtAssertQuantity;
3015  // this is option to print log of quantity rather than linear. default is
3016  // linear, and log only for a few such as ionization to see small numbers
3017  if( lgQuantityLog )
3018  {
3019  if( PredQuan > 0. )
3020  prtPredQuan = log10( MAX2( SMALLDOUBLE , PredQuan ) );
3021  else
3022  prtPredQuan = -37.;
3023  prtAssertQuantity = log10( MAX2( SMALLDOUBLE , AssertQuantity ) );
3024  }
3025  else
3026  {
3027  prtPredQuan = PredQuan;
3028  prtAssertQuantity = AssertQuantity;
3029  }
3030 
3031  /* start of line, flag problems */
3032  if( lg1OK )
3033  {
3034  // the ChkMonitor is a unique label so that we can grep on all output
3035  // and see what happened, without picking up input stream
3036  double relative = fabs( RelError / SDIV( fabs(AssertError)));
3037 
3038  if( relative < 0.25 || chAssertLimit != '=' )
3039  {
3040  fprintf( ioMONITOR, " ChkMonitor ");
3041  }
3042  else if( relative < 0.50 )
3043  {
3044  fprintf( ioMONITOR, " ChkMonitor - ");
3045  }
3046  else if( relative < 0.75 )
3047  {
3048  fprintf( ioMONITOR, " ChkMonitor -- ");
3049  }
3050  else if( relative < 0.90 )
3051  {
3052  fprintf( ioMONITOR, " ChkMonitor --- ");
3053  }
3054  else if( relative < 0.95 )
3055  {
3056  fprintf( ioMONITOR, " ChkMonitor ---- ");
3057  }
3058  else if( relative < 0.98 )
3059  {
3060  fprintf( ioMONITOR, " ChkMonitor ----- ");
3061  }
3062  else
3063  {
3064  fprintf( ioMONITOR, " ChkMonitor ------ ");
3065  }
3066  }
3067  else
3068  {
3069  fprintf( ioMONITOR, " ChkMonitor botch>>");
3070  }
3071 
3072  fprintf( ioMONITOR , "%-*s ", NCHLAB-1, chAssertLineLabel.c_str() );
3073 
3074  /* special formatting for the emission lines */
3075  if( strcmp( chAssertType, "Ll" )==0 || strcmp( chAssertType, "Lr" )==0 ||
3076  strcmp( chAssertType, "Lb" )==0 )
3077  {
3078  prt_wl( ioMONITOR , wavelength );
3079  }
3080  else
3081  {
3082  fprintf( ioMONITOR , "%*i", LineSave.wl_length, (int)wavelength );
3083  }
3084 
3085  const char* format = " %10.4f %c %10.4f %7.3f %7.3f ";
3086  if( lgPrtSciNot )
3087  format = " %10.3e %c %10.3e %7.3f %7.3f ";
3088 
3089  fprintf( ioMONITOR, format,
3090  prtPredQuan,
3091  chAssertLimit,
3092  prtAssertQuantity,
3093  RelError,
3094  AssertError);
3095 
3096  prtLineType( ioMONITOR, iLineType );
3097 
3098  // if botched and the botch is > 3 sigma, say BIG BOTCH,
3099  // the lg1OK is needed since some tests (number of zones, etc)
3100  // are limits, not the quantity, and if limit is large the
3101  // miss will be big too
3102 
3103  /* >>chng 02 nov 27, added lgFound so don't get big botch when line simply missing */
3104  if( !lg1OK && (fabs(RelError) > 3.*AssertError) && lgFound )
3105  {
3106  /* >>chng 20 aug 03, discern between scope of big botch */
3107  if( fabs(RelError) > 9.*AssertError )
3108  fprintf( ioMONITOR , " <<BIG BIG (>9sig) BOTCH!!\n");
3109  else
3110  fprintf( ioMONITOR , " <<BIG (3sig) BOTCH!!\n");
3111  lgBigBotch = true;
3112  }
3113  else
3114  {
3115  fprintf( ioMONITOR , "\n");
3116  }
3117 
3118  return;
3119 }
monitor_results.h
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:150
cdLine_ip
void cdLine_ip(long int ipLine, double *relint, double *absint)
Definition: cddrive.cpp:1110
t_radius::Radius
double Radius
Definition: radius.h:31
t_atmdat::WaveLengthCaseB
realnum WaveLengthCaseB[8][25][24]
Definition: atmdat.h:355
t_radius::depth
double depth
Definition: radius.h:31
t_thermal::htot
double htot
Definition: thermal.h:169
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
lgSpaceAllocated
static bool lgSpaceAllocated
Definition: monitor_results.cpp:44
InitMonitorResults
void InitMonitorResults(void)
Definition: monitor_results.cpp:99
thermal
t_thermal thermal
Definition: thermal.cpp:6
chIonLbl
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:233
colden
t_colden colden
Definition: colden.cpp:5
t_iso_sp::st
qList st
Definition: iso.h:482
Parser::FFmtRead
double FFmtRead(void)
Definition: parser.cpp:472
WavlenErrorGet
realnum WavlenErrorGet(realnum wavelength, long sig_figs)
Definition: lines_service.cpp:196
exp10
double exp10(double x)
Definition: cddefines.h:1368
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:65
ForcePass
double ForcePass(char chAssertLimit1)
Definition: monitor_results.cpp:84
lgQuantityLog
static vector< bool > lgQuantityLog
Definition: monitor_results.cpp:74
blend_iterator
map< std::string, std::vector< TransitionProxy > >::iterator blend_iterator
Definition: transition.h:668
get_ptr
T * get_ptr(T *v)
Definition: cddefines.h:1109
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
multi_arr::alloc
void alloc()
Definition: container_classes.h:1113
get_error_ratio
double get_error_ratio(double pred, double assert)
Definition: monitor_results.cpp:1696
struc
t_struc struc
Definition: struc.cpp:6
t_version::lgRelease
bool lgRelease
Definition: version.h:28
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:246
MonitorResults
t_monitorresults MonitorResults
Definition: monitor_results.cpp:39
t_cpu::i
t_cpu_i & i()
Definition: cpu.h:419
t_hcmap::temap
double * temap
Definition: hcmap.h:42
t_monitorresults::nSumErrorCaseMonitor
long int nSumErrorCaseMonitor
Definition: monitor_results.h:32
ParseMonitorResults
void ParseMonitorResults(Parser &p)
Definition: monitor_results.cpp:116
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:250
cdIonFrac
int cdIonFrac(const char *chLabel, long int IonStage, double *fracin, const char *chWeight, bool lgDensity)
Definition: cddrive.cpp:909
conv
t_conv conv
Definition: conv.cpp:5
t_optimize::nOptimiz
long int nOptimiz
Definition: optimize.h:250
t_conv::nTotalIoniz_start
long int nTotalIoniz_start
Definition: conv.h:164
Parser::GetQuote
int GetQuote(string &chLabel)
Definition: parser.cpp:213
chAssertLimit
static vector< char > chAssertLimit
Definition: monitor_results.cpp:64
LineSave
t_LineSave LineSave
Definition: lines.cpp:10
Parser::nMatchErase
bool nMatchErase(const char *chKey)
Definition: parser.h:173
cdTemp
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition: cddrive.cpp:1311
qStateConstProxy::DepartCoef
double DepartCoef() const
Definition: quantumstate.h:311
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
Parser::getWave
double getWave()
Definition: parser.cpp:379
nzone
long int nzone
Definition: cddefines.cpp:14
t_called::lgTalk
bool lgTalk
Definition: called.h:12
molezone
Definition: mole.h:378
ErrorDefaultPerformance
static realnum ErrorDefaultPerformance
Definition: monitor_results.cpp:53
Param
static multi_arr< double, 2 > Param
Definition: monitor_results.cpp:71
MIN2
#define MIN2(a, b)
Definition: cddefines.h:803
Parser
Definition: parser.h:43
cap4
void cap4(char *chCAP, const char *chLab)
Definition: service.cpp:261
strAssertSpecies
static vector< string > strAssertSpecies
Definition: monitor_results.cpp:78
t_LineSave::wl_length
int wl_length
Definition: lines.h:114
t_LineSave::lines
vector< LinSv > lines
Definition: lines.h:132
trimTrailingWhiteSpace
void trimTrailingWhiteSpace(string &str)
Definition: service.cpp:153
Parser::StringError
NORETURN void StringError() const
Definition: parser.cpp:203
dense
t_dense dense
Definition: global.cpp:15
Singleton< t_version >::Inst
static t_version & Inst()
Definition: cddefines.h:209
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
prt_line_err
void prt_line_err(FILE *ioOUT, const char *label, realnum wvlng)
Definition: prt.cpp:161
t_radius::depth_x_fillfac
double depth_x_fillfac
Definition: radius.h:80
t_rfield::EnergyIncidCont
realnum EnergyIncidCont
Definition: rfield.h:467
wind
Wind wind
Definition: wind.cpp:5
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
iteration
long int iteration
Definition: cddefines.cpp:16
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:554
diatomics::ortho_density
double ortho_density
Definition: h2_priv.h:326
lines_service.h
t_monitorresults::SumErrorCaseMonitor
double SumErrorCaseMonitor
Definition: monitor_results.h:31
cdH2_colden
double cdH2_colden(long iVib, long iRot)
Definition: mole_h2.cpp:2320
molezone::levels
qList * levels
Definition: mole.h:417
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:536
lgMonitorsOK
bool lgMonitorsOK
Definition: monitor_results.cpp:38
AssertError
static vector< double > AssertError
Definition: monitor_results.cpp:70
diatomics::para_density
double para_density
Definition: h2_priv.h:326
POW2
#define POW2
Definition: cddefines.h:979
multi_arr< double, 2 >
STATIC
#define STATIC
Definition: cddefines.h:118
nAsserts
static long nAsserts
Definition: monitor_results.cpp:75
t_monitorresults
Definition: monitor_results.h:30
AssertQuantity
static vector< double > AssertQuantity
Definition: monitor_results.cpp:70
pressure
t_pressure pressure
Definition: pressure.cpp:9
rfield
t_rfield rfield
Definition: rfield.cpp:9
input_readvector
void input_readvector(const char *chFile, vector< double > &vec, bool *lgError)
Definition: input.cpp:246
realnum
float realnum
Definition: cddefines.h:124
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
prtLineType
STATIC void prtLineType(FILE *ioMONITOR, const int iLineType)
Definition: monitor_results.cpp:2980
qList::size
size_t size() const
Definition: quantumstate.h:121
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
t_pressure::RadBetaMax
realnum RadBetaMax
Definition: pressure.h:96
getLevelsGeneric
const molezone * getLevelsGeneric(const char *chLabel, bool lgValidate, vector< long > &LevelList)
Definition: generic_state.cpp:103
Parser::NoNumb
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:345
ErrorDefault
static realnum ErrorDefault
Definition: monitor_results.cpp:50
Parser::GetElem
long int GetElem(void) const
Definition: parser.cpp:321
h2
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
assertBlends
static vector< std::vector< TransitionProxy > * > assertBlends
Definition: monitor_results.cpp:61
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1011
t_prt::lgPrintTime
bool lgPrintTime
Definition: prt.h:161
t_LineSave::sig_figs
long int sig_figs
Definition: lines.h:109
optimize
t_optimize optimize
Definition: optimize.cpp:6
t_elementnames::chElementNameShort
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
radius
t_radius radius
Definition: radius.cpp:5
prt
t_prt prt
Definition: prt.cpp:14
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:159
t_hcmap::nmap
long int nmap
Definition: hcmap.h:39
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:160
chAssertLineLabel
static vector< string > chAssertLineLabel
Definition: monitor_results.cpp:55
qStateConstProxy::status
enum level_status status() const
Definition: quantumstate.h:359
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
lgPrtSciNot
bool lgPrtSciNot
Definition: monitor_results.cpp:38
PrtOneMonitor
void PrtOneMonitor(FILE *ioMONITOR, const char *chAssertType, const string chAssertLineLabel, const realnum wavelength, const int iLineType, const double PredQuan, const char chAssertLimit, const double AssertQuantity, const double RelError, const double AssertError, const bool lg1OK, const bool lgQuantityLog, const bool lgFound)
Definition: monitor_results.cpp:3008
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1067
multi_arr::reserve
void reserve(size_type i1)
Definition: container_classes.h:1077
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:64
t_hcmap::hmap
double * hmap
Definition: hcmap.h:45
t_hmi::hmidep
double hmidep
Definition: hmi.h:42
cdColm
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition: cddrive.cpp:592
t_optimize::lgOptimize
bool lgOptimize
Definition: optimize.h:257
lgBigBotch
bool lgBigBotch
Definition: monitor_results.cpp:38
Parser::GetParam
bool GetParam(const char *chKey, double *val)
Definition: parser.h:154
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
powpq
double powpq(double x, int p, int q)
Definition: service.cpp:630
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:350
generic_state.h
t_cpu_i::lgMPI
bool lgMPI() const
Definition: cpu.h:391
wavelength
static vector< realnum > wavelength
Definition: monitor_results.cpp:76
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
AssertQuantity2
static vector< double > AssertQuantity2
Definition: monitor_results.cpp:70
t_rfield::EnergyDiffCont
realnum EnergyDiffCont
Definition: rfield.h:467
t_dense::eden
double eden
Definition: dense.h:201
lgInitDone
static bool lgInitDone
Definition: monitor_results.cpp:42
Parser::lgEOL
bool lgEOL(void) const
Definition: parser.h:113
NCHLAB
const int NCHLAB
Definition: cddefines.h:304
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
lgCheckMonitors
bool lgCheckMonitors(FILE *ioMONITOR)
Definition: monitor_results.cpp:1703
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
blends
map< std::string, std::vector< TransitionProxy > > blends
Definition: transition.cpp:36
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:33
t_hcmap::cmap
double * cmap
Definition: hcmap.h:48
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
t_atmdat::CaseBWlHeI
vector< realnum > CaseBWlHeI
Definition: atmdat.h:358
hcmap
t_hcmap hcmap
Definition: hcmap.cpp:23
t_version::lgReleaseBranch
bool lgReleaseBranch
Definition: version.h:25
Parser::PrintLine
int PrintLine(FILE *fp) const
Definition: parser.h:206
t_LineSave::nsum
long int nsum
Definition: lines.h:87
t_struc::testr
realnum * testr
Definition: struc.h:25
NASSERTS
static const int NASSERTS
Definition: monitor_results.cpp:47
t_LineSave::sig_figs_max
static const long sig_figs_max
Definition: lines.h:110
caps
void caps(char *chCard)
Definition: service.cpp:295
gv
GrainVar gv
Definition: grainvar.cpp:5
hmi
t_hmi hmi
Definition: hmi.cpp:5
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
cpu
static t_cpu cpu
Definition: cpu.h:427
ShowMe
void ShowMe(void)
Definition: service.cpp:205
iLineType
static vector< int > iLineType
Definition: monitor_results.cpp:59
Parser::GetRange
bool GetRange(const char *chKey, double *val1, double *val2)
Definition: parser.h:163
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:349
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:32
chAssertType
static char ** chAssertType
Definition: monitor_results.cpp:67
null_molezone
molezone * null_molezone
Definition: mole_species.cpp:51
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:529
Wind::windv
realnum windv
Definition: wind.h:18
lines_table
int lines_table()
Definition: parse_table.cpp:1727
ipLine
static long int * ipLine
Definition: prt_linesum.cpp:14
called
t_called called
Definition: called.cpp:4
t_hcmap::nMapAlloc
long int nMapAlloc
Definition: hcmap.h:53
GrainVar::rate_h2_form_grains_used_total
double rate_h2_form_grains_used_total
Definition: grainvar.h:574
t_struc::pressure
realnum * pressure
Definition: struc.h:25
t_thermal::ctot
double ctot
Definition: thermal.h:130
Generated by   doxygen 1.8.5