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-2017 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  /* look for element label within quotes,
1078  * will be null terminated */
1079  else if( (nelem = p.GetElem()) >= 0 )
1080  {
1081  /* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem])
1082  * into array that will be used to get ionization after calculation */
1083  chAssertLineLabel[nAsserts] = elementnames.chElementNameShort[nelem];
1084  trimTrailingWhiteSpace( chAssertLineLabel[nAsserts] );
1085  }
1086  else
1087  {
1088  fprintf( ioQQQ,
1089  " I could not identify an element name on this line.\n");
1090  fprintf( ioQQQ, " Sorry.\n" );
1091  cdEXIT(EXIT_FAILURE);
1092  }
1093 
1094  /* now get ionization stage (or grain type), which will be saved into wavelength */
1095  if( chAssertLineLabel[nAsserts] == "face" )
1096  {
1097  /* this option is temperature at illuminated face so no element */
1098  wavelength[nAsserts] = 0;
1099  }
1100  else if( chAssertLineLabel[nAsserts] == "H2" )
1101  {
1102  int j;
1103  /* this option is temperature of H2 so element is zero */
1104  wavelength[nAsserts] = 0;
1105  /* first find the 2 in H2 */
1106  j = (int)p.FFmtRead();
1107  if( j != 2 )
1108  TotalInsanity();
1109  }
1110  else
1111  {
1112  wavelength[nAsserts] = (realnum)p.FFmtRead();
1113  if( p.lgEOL() )
1114  {
1115  p.NoNumb("ionization stage");
1116  }
1117  /* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */
1118  if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 )
1119  {
1120  fprintf( ioQQQ,
1121  " The ionization stage is inappropriate for this element.\n");
1122  fprintf( ioQQQ, " Sorry.\n" );
1123  cdEXIT(EXIT_FAILURE);
1124  }
1125  }
1126 
1127  if( p.nMatch( "GRID" ) )
1128  {
1129  AssertQuantity[nAsserts] = AssertVector[optimize.nOptimiz];
1130  }
1131  else
1132  {
1133  /* now get temperature, log if number is <= 10, else linear */
1134  AssertQuantity[nAsserts] = p.FFmtRead();
1135  if( p.lgEOL() )
1136  p.NoNumb("temperature");
1137  }
1138 
1139  /* optional error, default available */
1140  AssertError[nAsserts] = p.FFmtRead();
1141  if( p.lgEOL() )
1142  AssertError[nAsserts] = ErrorDefault;
1143 
1144  /* now make sure we end up with positive linear temperature
1145  * number is log if <=10 unless linear keyword appears */
1146  if( AssertQuantity[nAsserts] <= 10. && !p.nMatch( "LINE" ) )
1147  {
1148  /* log since negative or 0 */
1149  AssertQuantity[nAsserts] =
1150  exp10(AssertQuantity[nAsserts] );
1151  /* entered as a log, so print as a log too */
1152  lgQuantityLog[nAsserts] = true;
1153  }
1154  }
1155 
1156  /* monitor log of helium hydrogen ionization correction factor */
1157  else if( p.nMatch("HHEI") )
1158  {
1159  /* this flag will mean H-He icf */
1160  strcpy( chAssertType[nAsserts] , "hh" );
1161  /* say that this is zone numbering */
1162  chAssertLineLabel[nAsserts] = "HHei" ;
1163 
1164  /* now get the ionization correction factor, it is always the linear
1165  * quantity itself, since can be positive or negative*/
1166  AssertQuantity[nAsserts] = p.FFmtRead();
1167  if( p.lgEOL() )
1168  {
1169  p.NoNumb("ionization correction factor");
1170  }
1171  /* will not use this */
1172  wavelength[nAsserts] = 0;
1173 
1174  /* optional error, default available (cannot do before loop since we
1175  * do not know how many numbers are on line */
1176  AssertError[nAsserts] = p.FFmtRead();
1177  if( p.lgEOL() )
1178  {
1179  /* default error was set above */
1180  AssertError[nAsserts] = ErrorDefault;
1181  }
1182  }
1183 
1184  /* this large H2 molecule */
1185  else if( p.nMatch(" H2 ") )
1186  {
1187  /* ortho to para ratio for last computed zone */
1188  if( p.nMatch("ORTH") )
1189  {
1190  /* this flag will mean ortho to para ratio */
1191  strcpy( chAssertType[nAsserts] , "or" );
1192  /* say that this is ortho to para density ratio */
1193  chAssertLineLabel[nAsserts] = "orth" ;
1194  }
1195  else
1196  {
1197  fprintf( ioQQQ,
1198  " I could not identify a second keyword on this line.\n");
1199  fprintf( ioQQQ, " Sorry.\n" );
1200  cdEXIT(EXIT_FAILURE);
1201  }
1202 
1203  /* now get the first number, which better be the 2 in H2 */
1204  n2 = (long)p.FFmtRead();
1205  if( p.lgEOL() || n2 != 2 )
1206  {
1207  p.NoNumb("the 2 in H2 ?!");
1208  }
1209  AssertQuantity[nAsserts] = p.FFmtRead();
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  /* monitor we are running in MPI mode */
1224  else if( p.nMatch(" MPI") )
1225  {
1226  /* this flag will mean number of zones */
1227  strcpy( chAssertType[nAsserts] , "mp" );
1228  /* say that this is zone numbering */
1229  chAssertLineLabel[nAsserts] = "mpi" ;
1230 
1231  wavelength[nAsserts] = 0.;
1232  AssertQuantity[nAsserts] = 1.;
1233  AssertError[nAsserts] = ErrorDefault;
1234  }
1235 
1236  /* monitor number of zones */
1237  else if( p.nMatch("NZON") )
1238  {
1239  /* this flag will mean number of zones */
1240  strcpy( chAssertType[nAsserts] , "z " );
1241  /* say that this is zone numbering */
1242  chAssertLineLabel[nAsserts] = "zone" ;
1243 
1244  /* now get number of zones */
1245  wavelength[nAsserts] = 0.;
1246  AssertQuantity[nAsserts] = p.FFmtRead();
1247  if( p.lgEOL() )
1248  p.NoNumb("zone number");
1249 
1250  /* optional error */
1251  AssertError[nAsserts] = p.FFmtRead();
1252  if( p.lgEOL() )
1253  AssertError[nAsserts] = ErrorDefault;
1254  }
1255 
1256  /* monitor (probably upper limit to) error in pressure across model */
1257  else if( p.nMatch("PRES") && p.nMatch("ERRO") )
1258  {
1259  /* this flag indicates ratio of standard deviation to the mean pressure */
1260  strcpy( chAssertType[nAsserts] , "pr" );
1261  /* say that this is error in pressure */
1262  chAssertLineLabel[nAsserts] = "pres" ;
1263 
1264  /* now get the pressure error, which will be saved into wavelength
1265  * in nearly all cases this is limit to error */
1266  wavelength[nAsserts] = 0;
1267  AssertQuantity[nAsserts] = (double)p.FFmtRead();
1268  if( p.lgEOL() )
1269  {
1270  p.NoNumb("pressure error");
1271  }
1272  else if( AssertQuantity[nAsserts] <= 0.)
1273  {
1274  /* number <= 0 is log of error */
1275  AssertQuantity[nAsserts] = exp10(AssertQuantity[nAsserts]);
1276  }
1277 
1278  /* optional error, default available (cannot do before loop since we
1279  * do not know how many numbers are on line */
1280  AssertError[nAsserts] = p.FFmtRead();
1281  if( p.lgEOL() )
1282  {
1283  /* default error was set above */
1284  AssertError[nAsserts] = ErrorDefault;
1285  }
1286  }
1287 
1288  else if( p.nMatch("PRADMAX") )
1289  {
1290  /* monitor pradmax - max ratio of rad to gas pressure */
1291  /* this flag indicates ratio of rad to gas pressure */
1292  strcpy( chAssertType[nAsserts] , "RM" );
1293  /* say that this is error in pressure */
1294  chAssertLineLabel[nAsserts] = "Prad" ;
1295 
1296  /* now get the pressure error, which will be saved into wavelength
1297  * in nearly all cases this is limit to error */
1298  wavelength[nAsserts] = 0;
1299  AssertQuantity[nAsserts] = (double)p.FFmtRead();
1300  if( p.lgEOL() )
1301  {
1302  p.NoNumb("PRADMAX");
1303  }
1304  else if( AssertQuantity[nAsserts] <= 0.)
1305  {
1306  /* number <= 0 is log of error */
1307  AssertQuantity[nAsserts] = exp10(AssertQuantity[nAsserts]);
1308  }
1309 
1310  /* optional error, default available (cannot do before loop since we
1311  * do not know how many numbers are on line */
1312  AssertError[nAsserts] = p.FFmtRead();
1313  if( p.lgEOL() )
1314  {
1315  /* default error was set above */
1316  AssertError[nAsserts] = ErrorDefault;
1317  }
1318  }
1319 
1320  /* monitor secondary ionization rate, csupra */
1321  else if( p.nMatch("CSUP") )
1322  {
1323  /* this flag will mean secondary ionization, entered as log */
1324  strcpy( chAssertType[nAsserts] , "sc" );
1325  /* say that this is sec ioniz */
1326  chAssertLineLabel[nAsserts] = "sion" ;
1327 
1328  /* now get rate, saved into monitor quantity */
1329  AssertQuantity[nAsserts] = p.FFmtRead();
1330  if( p.lgEOL() )
1331  {
1332  p.NoNumb("secondary ionization rate");
1333  }
1334  /* entered as log, make linear */
1335  AssertQuantity[nAsserts] = exp10( AssertQuantity[nAsserts] );
1336 
1337  /* no wavelength */
1338  wavelength[nAsserts] = 0;
1339 
1340  /* optional error, default available (cannot do before loop since we
1341  * do not know how many numbers are on line */
1342  AssertError[nAsserts] = p.FFmtRead();
1343  if( p.lgEOL() )
1344  {
1345  /* default error was set above */
1346  AssertError[nAsserts] = ErrorDefault;
1347  }
1348 
1349  /* we want to print the log of eden, not linear value */
1350  lgQuantityLog[nAsserts] = true;
1351  }
1352 
1353  /* monitor heating rate, erg/cm3/s, htot */
1354  else if( p.nMatch("HTOT") )
1355  {
1356  /* this flag will mean heating, entered as log */
1357  strcpy( chAssertType[nAsserts] , "ht" );
1358 
1359  /* say that this is heating rate */
1360  chAssertLineLabel[nAsserts] = "htot" ;
1361 
1362  /* now get rate, saved into monitor quantity */
1363  AssertQuantity[nAsserts] = p.FFmtRead();
1364  if( p.lgEOL() )
1365  {
1366  p.NoNumb("heating rate");
1367  }
1368  /* entered as log, make linear */
1369  AssertQuantity[nAsserts] = exp10( AssertQuantity[nAsserts] );
1370 
1371  /* no wavelength */
1372  wavelength[nAsserts] = 0;
1373 
1374  /* optional error, default available (cannot do before loop since we
1375  * do not know how many numbers are on line */
1376  AssertError[nAsserts] = p.FFmtRead();
1377  if( p.lgEOL() )
1378  {
1379  /* default error was set above */
1380  AssertError[nAsserts] = ErrorDefault;
1381  }
1382 
1383  /* we want to print the log of the heating, not linear value */
1384  lgQuantityLog[nAsserts] = true;
1385  }
1386 
1387  /* monitor cooling rate, erg/cm3/s, ctot */
1388  else if( p.nMatch("CTOT") )
1389  {
1390  /* this flag will mean cooling */
1391  strcpy( chAssertType[nAsserts] , "ct" );
1392 
1393  /* say that this is cooling rate */
1394  chAssertLineLabel[nAsserts] = "ctot";
1395 
1396  /* Look for GRID command */
1397  if( p.nMatch( "GRID" ) )
1398  {
1399  AssertQuantity[nAsserts] = AssertVector[optimize.nOptimiz];
1400  }
1401  else
1402  {
1403  AssertQuantity[nAsserts] = p.FFmtRead();
1404  /* now get rate, saved into monitor quantity */
1405  if( p.lgEOL() )
1406  {
1407  p.NoNumb("cooling rate");
1408  }
1409  }
1410  /* entered as log, make linear */
1411  AssertQuantity[nAsserts] = exp10( AssertQuantity[nAsserts] );
1412 
1413  /* no wavelength */
1414  wavelength[nAsserts] = 0;
1415 
1416  /* optional error, default available (cannot do before loop since we
1417  * do not know how many numbers are on line */
1418  AssertError[nAsserts] = p.FFmtRead();
1419  if( p.lgEOL() )
1420  {
1421  /* default error was set above */
1422  AssertError[nAsserts] = ErrorDefault;
1423  }
1424 
1425  /* we want to print the log of the heating, not linear value */
1426  lgQuantityLog[nAsserts] = true;
1427  }
1428 
1429  /* monitor number of iterations per zone, a test of convergence */
1430  else if( p.nMatch("ITRZ") )
1431  {
1432  /* this flag will mean number of iterations per zone */
1433  strcpy( chAssertType[nAsserts] , "iz" );
1434  /* say that this is iterations per zone */
1435  chAssertLineLabel[nAsserts] = "itrz" ;
1436 
1437  /* now get quantity */
1438  AssertQuantity[nAsserts] = p.FFmtRead();
1439  if( p.lgEOL() )
1440  {
1441  p.NoNumb("iterations per zone");
1442  }
1443  /* wavelength is meaningless */
1444  wavelength[nAsserts] = 0;
1445 
1446  /* optional error, default available */
1447  AssertError[nAsserts] = p.FFmtRead();
1448  if( p.lgEOL() )
1449  AssertError[nAsserts] = ErrorDefaultPerformance;
1450  }
1451 
1452  /* monitor electron density of the last zone */
1453  else if( p.nMatch("EDEN") )
1454  {
1455  /* this flag will mean electron density of the last zone */
1456  strcpy( chAssertType[nAsserts] , "e " );
1457  /* say that this is electron density */
1458  chAssertLineLabel[nAsserts] = "eden" ;
1459 
1460  /* now get electron density, which is a log */
1461  AssertQuantity[nAsserts] =
1462  exp10( p.FFmtRead() );
1463  if( p.lgEOL() )
1464  {
1465  p.NoNumb(" electron density of the last zone");
1466  }
1467 
1468  /* optional error, default available (cannot do before loop since we
1469  * do not know how many numbers are on line */
1470  AssertError[nAsserts] = p.FFmtRead();
1471  if( p.lgEOL() )
1472  {
1473  /* default error was set above */
1474  AssertError[nAsserts] = ErrorDefault;
1475  }
1476  wavelength[nAsserts] = 0;
1477 
1478  /* we want to print the log of eden, not linear value */
1479  lgQuantityLog[nAsserts] = true;
1480  }
1481 
1482  /* monitor energy density - Tu - of last zone */
1483  else if( p.nMatch(" TU ") )
1484  {
1485  /* this flag will mean energy density of the last zone */
1486  strcpy( chAssertType[nAsserts] , "Tu" );
1487  chAssertLineLabel[nAsserts] = "Tu" ;
1488 
1489  /* get temperature */
1490  AssertQuantity[nAsserts] = p.FFmtRead();
1491  if( p.lgEOL() )
1492  p.NoNumb("energy density of last zone");
1493 
1494  /* now make sure we end up with positive linear temperature
1495  * number is log if <=10 unless linear keyword appears */
1496  lgQuantityLog[nAsserts] = false;
1497  if( AssertQuantity[nAsserts] <= 10. && !p.nMatch( "LINE" ) )
1498  {
1499  /* log since negative or 0 */
1500  AssertQuantity[nAsserts] =
1501  exp10(AssertQuantity[nAsserts] );
1502  /* entered as a log, so print as a log too */
1503  }
1504 
1505  /* optional error, default available (cannot do before loop since we
1506  * do not know how many numbers are on line */
1507  AssertError[nAsserts] = p.FFmtRead();
1508  if( p.lgEOL() )
1509  /* default error was set above */
1510  AssertError[nAsserts] = ErrorDefault;
1511  wavelength[nAsserts] = 0;
1512  }
1513 
1514  /* monitor thickness or depth of model */
1515  else if( p.nMatch("THIC") || p.nMatch("DEPT") )
1516  {
1517  /* this flag will mean thickness or depth */
1518  strcpy( chAssertType[nAsserts] , "th" );
1519  /* say that this is thickness */
1520  chAssertLineLabel[nAsserts] = "thic" ;
1521 
1522  /* now get thickness, which is a log */
1523  AssertQuantity[nAsserts] = exp10( p.FFmtRead() );
1524  if( p.lgEOL() )
1525  {
1526  p.NoNumb("thickness or depth of model");
1527  }
1528 
1529  /* optional error, default available (cannot do before loop since we
1530  * do not know how many numbers are on line */
1531  AssertError[nAsserts] = p.FFmtRead();
1532  if( p.lgEOL() )
1533  {
1534  /* default error was set above */
1535  AssertError[nAsserts] = ErrorDefault;
1536  }
1537  wavelength[nAsserts] = 0;
1538 
1539  /* we want to print the log of eden, not linear value */
1540  lgQuantityLog[nAsserts] = true;
1541  }
1542 
1543  /* monitor outer radius of model */
1544  else if( p.nMatch("RADI") )
1545  {
1546  /* this flag will mean radius */
1547  strcpy( chAssertType[nAsserts] , "ra" );
1548  /* say that this is radius */
1549  chAssertLineLabel[nAsserts] = "radi" ;
1550 
1551  /* now get radius, which is a log */
1552  AssertQuantity[nAsserts] = exp10( p.FFmtRead() );
1553  if( p.lgEOL() )
1554  {
1555  p.NoNumb("outer radius");
1556  }
1557 
1558  /* optional error, default available (cannot do before loop since we
1559  * do not know how many numbers are on line */
1560  AssertError[nAsserts] = p.FFmtRead();
1561  if( p.lgEOL() )
1562  {
1563  /* default error was set above */
1564  AssertError[nAsserts] = ErrorDefault;
1565  }
1566  wavelength[nAsserts] = 0;
1567 
1568  /* we want to print the log of radius, not linear value */
1569  lgQuantityLog[nAsserts] = true;
1570  }
1571 
1572  /* monitor number of iterations */
1573  else if( p.nMatch("NITE") )
1574  {
1575  /* this flag will mean number of iterations */
1576  strcpy( chAssertType[nAsserts] , "Z " );
1577  /* say that this is iteration numbering */
1578  chAssertLineLabel[nAsserts] = "iter" ;
1579 
1580  wavelength[nAsserts] = 0.f;
1581 
1582  /* now get number of iterations, which will be saved into wavelength */
1583  AssertQuantity[nAsserts] = p.FFmtRead();
1584  if( p.lgEOL() )
1585  {
1586  p.NoNumb("number of iterations");
1587  }
1588 
1589  /* optional error, default available (cannot do before loop since we
1590  * do not know how many numbers are on line */
1591  AssertError[nAsserts] = p.FFmtRead();
1592  if( p.lgEOL() )
1593  {
1594  /* default error was set above */
1595  AssertError[nAsserts] = ErrorDefault;
1596  }
1597  }
1598 
1599  /* monitor terminal velocity, at end of calculation
1600  * input in km/s and saved that way, even though code uses cm/s */
1601  else if( p.nMatch("VELO") )
1602  {
1603  /* this flag will mean velocity */
1604  strcpy( chAssertType[nAsserts] , "Ve" );
1605  /* say that this is velocity */
1606  chAssertLineLabel[nAsserts] = "vel" ;
1607  wavelength[nAsserts] = 0;
1608 
1609  /* now get velocity */
1610  AssertQuantity[nAsserts] = p.FFmtRead();
1611  if( p.lgEOL() )
1612  p.NoNumb("terminal velocity");
1613 
1614  /* optional error, default available (cannot do before loop since we
1615  * do not know how many numbers are on line */
1616  AssertError[nAsserts] = p.FFmtRead();
1617  if( p.lgEOL() )
1618  {
1619  /* default error was set above */
1620  AssertError[nAsserts] = ErrorDefault;
1621  }
1622  }
1623  /* monitor nothing - a pacifier */
1624  else if( p.nMatch("NOTH") )
1625  {
1626  strcpy( chAssertType[nAsserts] , "NO" );
1627  chAssertLineLabel[nAsserts] = "noth" ;
1628  wavelength[nAsserts] = 0;
1629  AssertQuantity[nAsserts] = 0.;
1630  AssertError[nAsserts] = ErrorDefault;
1631  }
1632  else
1633  {
1634  /* did not recognize a command */
1635  fprintf( ioQQQ,
1636  " Unrecognized command. The line image was\n");
1637  p.PrintLine(ioQQQ);
1638  fprintf( ioQQQ,
1639  " The options I know about are: ionization, line, departure coefficient, map, column, "
1640  "temperature, nzone, csupre, htot, itrz, eden, thickness, niter, \n");
1641  fprintf( ioQQQ, " Sorry.\n" );
1642  cdEXIT(EXIT_FAILURE);
1643  }
1644 
1645  /* increment number of asserts and confirm that limit not exceeded */
1646  ++nAsserts;
1647  if( nAsserts >= NASSERTS )
1648  {
1649  fprintf(ioQQQ,
1650  " ParseMonitorResults: too many asserts, limit is NASSERT=%d\n",
1651  NASSERTS );
1652  cdEXIT(EXIT_FAILURE);
1653  }
1654  return;
1655 }
1656 
1657 inline double get_error_ratio ( double pred, double assert )
1658 {
1659  return 1. - safe_div(pred,assert,1.);
1660 }
1661 
1662 /*============================================================================*/
1663 /*lgCheckMonitors, checks asserts, last thing cloudy calls, returns true if all are ok, false if problems */
1664 bool lgCheckMonitors(
1665  /* this is the file we will write this to, usually standard output,
1666  * but can be save */
1667  FILE * ioMONITOR )
1668 {
1669  double PredQuan[NASSERTS] , RelError[NASSERTS];
1670  /* this will be true if the quantity was found, and false if not. Used to prevent
1671  * big botch flag when quantity not found (as in removed old he atom) */
1672  bool lgFound[NASSERTS];
1673  double relint , absint;
1674  bool lg1OK[NASSERTS];
1675  long i,j;
1676  /* This structure is for reporting another close match for asserts of line
1677  * intensities only. The zeroth, first, and second elements for each monitor are,
1678  * respectively, the first, second, and third matches the code finds, if any.
1679  * A negative number means no match. A positive number indicates the pointer
1680  * in the line stack of that match. */
1681  /* use multi_arr here to prevent bogus array bounds violations being reported by pgCC */
1682  multi_arr<long,2> ipDisambiguate(NASSERTS,3);
1683  long lgDisambiguate = false;
1684  char chCaps[NCHLAB], chFind[NCHLAB];
1685  realnum errorwave;
1686 
1687  DEBUG_ENTRY( "lgCheckMonitors()" );
1688 
1689  /* this is a global variable in monitor_results.h, and can be checked by
1690  * other routines to see if asserts are ok - (most runs will not use asserts) */
1691  lgMonitorsOK = true;
1692 
1693  /* will be used if there were big botched monitors */
1694  lgBigBotch = false;
1695 
1696  /* the optimize*.in and stars_oppim*.in tests in the test suite include
1697  * asserts while optimizing. We do not want to test the asserts during
1698  * the optimization process, since we will find blown asserts and report
1699  * major problems. We do want to test asserts on the final model however.
1700  * Printout will usually be turned off in all except the final model,
1701  * so do not to the monitor tests if we are optimizing but not printing */
1702  if( !called.lgTalk && optimize.lgOptimize )
1703  {
1704  /* just return */
1705  return true;
1706  }
1707 
1708  /*fprintf(ioQQQ , "DEBUG grid %li\n", optimize.nOptimiz );*/
1709 
1710  /* this will usually just return, but with table lines will check
1711  * existence of some lines */
1712  if( lines_table() )
1713  {
1714  lgBigBotch = true;
1715  lgMonitorsOK = false;
1716  }
1717 
1718  /* first check all asserts, there probably are none */
1719  for(i=0; i<nAsserts; ++i )
1720  {
1721  lg1OK[i] = true;
1722  PredQuan[i] = 0.;
1723  RelError[i] = 0.;
1724  for(j=0; j<3; ++j )
1725  ipDisambiguate[i][j] = -1;
1726 
1727  /* this flag is set false if we don't find the requested quantity */
1728  lgFound[i] = true;
1729 
1730  /* which type of monitor? */
1731  /* is it intensity? */
1732  if( strcmp( chAssertType[i] , "Lr")==0 )
1733  {
1734  /* this is an intensity, get the line, returns false if could not find it */
1735  ipDisambiguate[i][0] = cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint, iLineType[i] );
1736  if( ipDisambiguate[i][0] <= 0 )
1737  {
1738  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1739  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1740 
1741  fprintf( ioMONITOR,
1742  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1743  /* go to next line */
1744  lg1OK[i] = false;
1745  RelError[i] = 100000.;
1746  PredQuan[i] = 0;
1747  lg1OK[i] = false;
1748  lgMonitorsOK = false;
1749  lgFound[i] = false;
1750  continue;
1751  }
1752  else
1753  {
1754  /********* LINE DISAMBIGUATION *************/
1755  /* Here we look for lines with same label and small wavelength
1756  * differences so that we can disambiguate below */
1757 
1758  /* change chLabel to all caps */
1759  strcpy(chFind, chAssertLineLabel[i].c_str());
1760  caps(chFind);
1761 
1762  /* get the error associated with specified significant figures */
1763  errorwave = WavlenErrorGet( wavelength[i], LineSave.sig_figs );
1764 
1765  /* go through rest of line stack to look for close matches */
1766  for( j=1; j < LineSave.nsum; j++ )
1767  {
1768  /* don't bother with this one, we've already identified it. */
1769  if( j==ipDisambiguate[i][0] )
1770  continue;
1771 
1772  /* change chLabel to all caps to be like input chALab */
1773  cap4(chCaps, LineSave.lines[j].chALab());
1774 
1775  /* look for wavelengths within 3 error bars.
1776  * For example, for a line entered in Angstroms with
1777  * four significant figures, the error bar is 0.5 Ang.
1778  * So here we will find any lines within 1.5 Angstroms
1779  * of the
1780  * asserted wavelength. check wavelength and chLabel for a match */
1781  if( fabs(LineSave.lines[j].wavelength()-wavelength[i]) < 3.f*errorwave )
1782  {
1783  /* now see if labels agree */
1784  if( strcmp(chCaps,chFind) == 0 )
1785  {
1786  double relint1, absint1, current_error;
1787 
1788  cdLine_ip( j, &relint1, &absint1, iLineType[i] );
1789 
1790  if (AssertError[i] > 0.)
1791  current_error = fabs(1. - relint1/AssertQuantity[i]);
1792  else
1793  current_error = relint1 - AssertQuantity[i];
1794 
1795  if( current_error < 2.*fabs(AssertError[i]) ||
1796  current_error < 2.*fabs(RelError[i]) )
1797  {
1798  lgDisambiguate = true;
1799  /* if second match (element 1) is already set,
1800  * this is third match (element 2). Set and break out. */
1801  if( ipDisambiguate[i][1] > 0 )
1802  {
1803  ipDisambiguate[i][2] = j;
1804  break;
1805  }
1806  else
1807  {
1808  ipDisambiguate[i][1] = j;
1809  }
1810  }
1811  }
1812  }
1813  }
1814  }
1815 
1816  PredQuan[i] = relint;
1817  if (AssertError[i] > 0.0)
1818  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1819  else
1820  RelError[i] = PredQuan[i]-AssertQuantity[i] ;
1821  }
1822 
1823  else if( strcmp( chAssertType[i], "Lb" ) == 0 )
1824  {
1825  if( cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint, iLineType[i] ) <= 0 )
1826  {
1827  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1828  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1829 
1830  fprintf( ioMONITOR,
1831  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1832  /* go to next line */
1833  RelError[i] = 10000000.;
1834  PredQuan[i] = 0;
1835  lg1OK[i] = false;
1836  lgFound[i] = false;
1837  lgMonitorsOK = false;
1838  continue;
1839  }
1840 
1841  double relint_cb = 0.,
1842  absint_cb = 0.;
1843  if( cdLine( "Ca B", wavelength[i], &relint_cb, &absint_cb, iLineType[i] ) <= 0 )
1844  {
1845  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1846  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1847 
1848  fprintf( ioMONITOR,
1849  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1850  /* go to next line */
1851  RelError[i] = 10000000.;
1852  PredQuan[i] = 0;
1853  lg1OK[i] = false;
1854  lgFound[i] = false;
1855  lgMonitorsOK = false;
1856  continue;
1857  }
1858 
1859  PredQuan[i] = absint / absint_cb;
1860  if (AssertError[i] > 0.0)
1861  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1862  else
1863  RelError[i] = PredQuan[i] - AssertQuantity[i];
1864  // printf("i=%ld\t %s\t %g\t pred = %g\t asse = %g\t err = %g\n",
1865  }
1866 
1867  /*this is line luminosity */
1868  else if( strcmp( chAssertType[i] , "Ll")==0 )
1869  {
1870 // long indice=cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint );
1871 // fprintf(ioQQQ,"index line %ld %s %g \n", indice, LineSave.lines[indice].chALab(), LineSave.lines[indice].wavelength());
1872  /* this is a luminosity, get the line, returns false if could not find it */
1873  if( cdLine( chAssertLineLabel[i].c_str(), wavelength[i], &relint, &absint, iLineType[i] ) <= 0 )
1874  //if (indice <= 0 )
1875  {
1876  fprintf( ioMONITOR, " monitor error: lgCheckMonitors could not find line ");
1877  prt_line_err( ioMONITOR, chAssertLineLabel[i].c_str(), wavelength[i] );
1878 
1879  fprintf( ioMONITOR,
1880  " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1881  /* go to next line */
1882  RelError[i] = 10000000.;
1883  PredQuan[i] = 0;
1884  lg1OK[i] = false;
1885  lgFound[i] = false;
1886  lgMonitorsOK = false;
1887  continue;
1888  }
1889  PredQuan[i] = absint;
1890  if (AssertError[i] > 0.0)
1891  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1892  else
1893  RelError[i] = PredQuan[i] - AssertQuantity[i];
1894  // printf("i=%ld\t %s\t %g\t pred = %g\t asse = %g\t err = %g\n",
1895  // i, chAssertLineLabel[i].c_str(), wavelength[i], PredQuan[i], AssertQuantity[i], RelError[i] );
1896  }
1897  else if( strcmp( chAssertType[i] , "hh" ) == 0)
1898  {
1899  double hfrac , hefrac;
1900  /* get H ionization fraction, returns false if could not find it */
1901  if( cdIonFrac(
1902  /* four char string, null terminated, giving the element name */
1903  "hydr",
1904  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
1905  1,
1906  /* will be fractional ionization */
1907  &hfrac,
1908  /* how to weight the average, must be "VOLUME" or "RADIUS" */
1909  "VOLUME" ,
1910  /* do not want extra factor of density */
1911  false) )
1912  {
1913  fprintf( ioMONITOR,
1914  " monitor error: lgCheckMonitors could not find h ionization fraction \n");
1915  lg1OK[i] = false;
1916  RelError[i] = 0;
1917  PredQuan[i] = 0;
1918  lgFound[i] = false;
1919  continue;
1920  }
1921  if( cdIonFrac(
1922  /* four char string, null terminated, giving the element name */
1923  "heli",
1924  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
1925  1,
1926  /* will be fractional ionization */
1927  &hefrac,
1928  /* how to weight the average, must be "VOLUME" or "RADIUS" */
1929  "VOLUME" ,
1930  /* do not want extra factor of density */
1931  false) )
1932  {
1933  fprintf( ioMONITOR,
1934  " monitor error: lgCheckMonitors could not find h ionization fraction \n");
1935  lg1OK[i] = false;
1936  RelError[i] = 0;
1937  PredQuan[i] = 0;
1938  lgFound[i] = false;
1939  continue;
1940  }
1941  /* the helium hydrogen ionization correction factor */
1942  PredQuan[i] = hefrac-hfrac;
1943  /* two icf's in difference, no need to div by mean since already on scale with unity */
1944  RelError[i] = fabs(AssertQuantity[i] - (hefrac-hfrac) );
1945  }
1946 
1947  else if( strcmp( chAssertType[i] , "mp" ) == 0)
1948  {
1949  PredQuan[i] = cpu.i().lgMPI() ? 1. : 0.;
1950  /* use absolute error */
1951  RelError[i] = AssertQuantity[i] - PredQuan[i];
1952  }
1953 
1954  else if( strcmp( chAssertType[i] , "z " ) == 0)
1955  {
1956  /* this is the number of zones */
1957  PredQuan[i] = (double)nzone;
1958  if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
1959  RelError[i] = ForcePass(chAssertLimit[i]);
1960  else
1961  {
1962  /* two integers in difference */
1963  if (AssertError[i] > 0.)
1964  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1965  else
1966  RelError[i] = PredQuan[i] - AssertQuantity[i];
1967  }
1968  }
1969 
1970  else if( strcmp( chAssertType[i] , "or" ) == 0)
1971  {
1972  /* ortho to para ratio for large H2 molecule in last zone */
1973  PredQuan[i] = h2.ortho_density / SDIV( h2.para_density );
1974 
1975  /* this is relative error */
1976  if (AssertError[i] > 0.)
1977  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1978  else
1979  RelError[i] = PredQuan[i] - AssertQuantity[i];
1980  }
1981 
1982  else if( strcmp( chAssertType[i] , "g2" ) == 0)
1983  {
1984  /* check Jura rate, rate per vol that H2 forms on grain surfaces */
1985  PredQuan[i] = gv.rate_h2_form_grains_used_total;
1986  /* this is relative error */
1987  if (AssertError[i] > 0.)
1988  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
1989  else
1990  RelError[i] = PredQuan[i] - AssertQuantity[i];
1991  }
1992 
1993  else if( strcmp( chAssertType[i] , "RM" ) == 0)
1994  {
1995  /* check Jura rate, rate per vol that H2 forms on grain surfaces */
1996  PredQuan[i] = pressure.RadBetaMax;
1997  /* this is relative error */
1998  if (AssertError[i] > 0.)
1999  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2000  else
2001  RelError[i] = PredQuan[i] - AssertQuantity[i];
2002  }
2003 
2004  else if( strcmp( chAssertType[i] , "pr" ) == 0)
2005  {
2006  /* standard deviation of the pressure */
2007  double sumx=0., sumx2=0., average;
2008  long int n;
2009  /* do sums to form standard deviation */
2010  for( n=0; n<nzone; n++ )
2011  {
2012  sumx += struc.pressure[n];
2013  sumx2 += POW2(struc.pressure[n]);
2014  }
2015  if( nzone>1 )
2016  {
2017  /* this is average */
2018  average = sumx/nzone;
2019  /* this is abs std */
2020  sumx = sqrt( (sumx2-POW2(sumx)/nzone)/(nzone-1) );
2021  /* save the relative std */
2022  PredQuan[i] = sumx / average;
2023  }
2024  else
2025  {
2026  PredQuan[i] = 0.;
2027  }
2028 
2029  // this is already relative error, do not need 1-ratio
2030  RelError[i] = PredQuan[i];
2031  }
2032 
2033  else if( strcmp( chAssertType[i] , "iz") == 0 )
2034  {
2035  /* this is number of iterations per zone, a test of convergence properties */
2036  if( nzone > 0 )
2037  PredQuan[i] = (double)(conv.nTotalIoniz-conv.nTotalIoniz_start)/(double)(nzone);
2038  else
2039  /* something big so monitor will botch. */
2040  PredQuan[i] = 1e10;
2041 
2042  if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
2043  RelError[i] = ForcePass(chAssertLimit[i]);
2044  else
2045  {
2046  /* this is relative error */
2047  if (AssertError[i] > 0.)
2048  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2049  else
2050  RelError[i] = PredQuan[i]- AssertQuantity[i];
2051  }
2052  }
2053 
2054  else if( strcmp( chAssertType[i] , "e ") == 0 )
2055  {
2056  /* this is electron density of the last zone */
2057  PredQuan[i] = dense.eden;
2058  /* this is relative error */
2059  if (AssertError[i] > 0.)
2060  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2061  else
2062  RelError[i] = PredQuan[i]- AssertQuantity[i];
2063  }
2064 
2065  else if( strcmp( chAssertType[i] , "Tu") == 0 )
2066  {
2067  /* this is radiation energy density of the last zone */
2068  PredQuan[i] = powpq((rfield.EnergyIncidCont+rfield.EnergyDiffCont)/
2069  (4.*STEFAN_BOLTZ),1,4);
2070  /* this is relative error */
2071  if (AssertError[i] > 0.)
2072  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2073  else
2074  RelError[i] = PredQuan[i]- AssertQuantity[i];
2075  }
2076 
2077  else if( strcmp( chAssertType[i] , "th") == 0 )
2078  {
2079  /* this is thickness */
2080  PredQuan[i] = radius.depth;
2081  /* this is relative error */
2082  if (AssertError[i] > 0.)
2083  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2084  else
2085  RelError[i] = PredQuan[i]- AssertQuantity[i];
2086  }
2087 
2088  else if( strcmp( chAssertType[i] , "ra") == 0 )
2089  {
2090  /* this is thickness */
2091  PredQuan[i] = radius.Radius;
2092  /* this is relative error */
2093  if (AssertError[i] > 0.)
2094  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2095  else
2096  RelError[i] = PredQuan[i]- AssertQuantity[i];
2097  }
2098 
2099  else if( strcmp( chAssertType[i] , "Ve") == 0 )
2100  {
2101  /* this is final velocity of wind in km/s (code uses cm/s) */
2102  PredQuan[i] = wind.windv/1e5;
2103  /* this is relative error */
2104  if (AssertError[i] > 0.)
2105  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2106  else
2107  RelError[i] = PredQuan[i]- AssertQuantity[i];
2108  }
2109 
2110  else if( strcmp( chAssertType[i] , "NO") == 0 )
2111  {
2112  /* this is nothing */
2113  PredQuan[i] = 0;
2114  /* this is relative error */
2115  RelError[i] = 0.;
2116  }
2117 
2118  else if( strcmp( chAssertType[i] , "sc") == 0 )
2119  {
2120  /* this is secondary ionization rate */
2121  PredQuan[i] = secondaries.csupra[ipHYDROGEN][0];
2122  /* this is relative error */
2123  if (AssertError[i] > 0.)
2124  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2125  else
2126  RelError[i] = PredQuan[i]- AssertQuantity[i];
2127  }
2128 
2129  else if( strcmp( chAssertType[i] , "ht") == 0 )
2130  {
2131  /* this is heating rate */
2132  PredQuan[i] = thermal.htot;
2133  /* this is relative error */
2134  if (AssertError[i] > 0.)
2135  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2136  else
2137  RelError[i] = PredQuan[i]- AssertQuantity[i];
2138  }
2139 
2140  else if( strcmp( chAssertType[i] , "ct") == 0 )
2141  {
2142  /* this is cooling rate */
2143  PredQuan[i] = thermal.ctot;
2144  /* this is relative error */
2145  if (AssertError[i] > 0.)
2146  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2147  else
2148  RelError[i] = PredQuan[i]- AssertQuantity[i];
2149  }
2150 
2151  else if( strcmp( chAssertType[i] , "Z ") == 0 )
2152  {
2153  /* this is the number of iterations */
2154  PredQuan[i] = (double)iteration;
2155  if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
2156  RelError[i] = ForcePass(chAssertLimit[i]);
2157  else
2158  {
2159  /* two integers in difference */
2160  if (AssertError[i] > 0.)
2161  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2162  else
2163  RelError[i] = PredQuan[i]- AssertQuantity[i];
2164  }
2165  }
2166 
2167  else if( strcmp( chAssertType[i] , "CB") == 0 )
2168  {
2169  long int nISOCaseB = (long)Param[i][0];
2170  long int nelemCaseB = (long)Param[i][1];
2171  string chElemLabelCaseB = chIonLbl( nelemCaseB+1, nelemCaseB+1-nISOCaseB );
2172 
2173  /* sets lowest quantum number index */
2174  int iCase;
2175  if( "Ca A" == chAssertLineLabel[i])
2176  iCase = 0;
2177  else if( "Ca B" == chAssertLineLabel[i] )
2178  iCase = 1;
2179  else
2180  TotalInsanity();
2181 
2182  iLineType[i] = 0;
2183 
2184  RelError[i] = 0.;
2185  long nHighestPrinted = iso_sp[nISOCaseB][nelemCaseB].n_HighestResolved_max;
2186  if( nISOCaseB == ipH_LIKE )
2187  {
2188  fprintf(ioMONITOR," Species nHi nLo Wl Computed Asserted error\n");
2189  /* limit of 10 is because that is all we printed and saved in prt_lines_hydro
2190  * wavelengths will come out of atmdat.WaveLengthCaseB - first index is
2191  * nelem on C scale, H is 0, second two are configurations of line on
2192  * physics scale, so Ha is 3-2 */
2193  for( long int ipLo=1+iCase; ipLo< MIN2(10,nHighestPrinted-1); ++ipLo )
2194  {
2195  for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi )
2196  {
2197  /* monitor the line */
2198  realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo];
2199  /* range option to restrict wavelength coverage */
2200  if( wl < Param[i][2] || wl > Param[i][3] )
2201  continue;
2202 
2203  double relint, absint, CBrelint, CBabsint;
2204  /* find the predicted line intensity */
2205  cdLine( chAssertLineLabel[i].c_str() , wl , &CBrelint , &CBabsint, iLineType[i] );
2206  if( CBrelint < Param[i][4] )
2207  continue;
2208  double error;
2209  /* now find the Case B intensity - may not all be present */
2210  if( cdLine( chElemLabelCaseB.c_str(), wl, &relint, &absint, iLineType[i] ) > 0 )
2211  {
2212  if (AssertError[i] > 0.)
2213  error = (CBabsint - absint)/MAX2(CBabsint , absint);
2214  else
2215  error = (CBabsint - absint);
2216  double RelativeError = fabs(error / AssertError[i]);
2217  /* start of line, flag problems */
2218  if( RelativeError < 1. )
2219  {
2220  if( RelativeError < 0.25 )
2221  {
2222  fprintf( ioMONITOR, " ChkMonitor ");
2223  }
2224  else if( RelativeError < 0.50 )
2225  {
2226  fprintf( ioMONITOR, " ChkMonitor - ");
2227  }
2228  else if( RelativeError < 0.75 )
2229  {
2230  fprintf( ioMONITOR, " ChkMonitor -- ");
2231  }
2232  else if( RelativeError < 0.90 )
2233  {
2234  fprintf( ioMONITOR, " ChkMonitor --- ");
2235  }
2236  else if( RelativeError < 0.95 )
2237  {
2238  fprintf( ioMONITOR, " ChkMonitor ---- ");
2239  }
2240  else if( RelativeError < 0.98 )
2241  {
2242  fprintf( ioMONITOR, " ChkMonitor ----- ");
2243  }
2244  else
2245  {
2246  fprintf( ioMONITOR, " ChkMonitor ------ ");
2247  }
2248 
2249  }
2250  else
2251  {
2252  fprintf( ioMONITOR, " ChkMonitor botch>>");
2253  }
2254  fprintf(ioMONITOR," %s %3li %3li ",
2255  chElemLabelCaseB.c_str() , ipHi , ipLo );
2256  prt_wl(ioMONITOR, wl );
2257  fprintf(ioMONITOR," %.2e %.2e %10.3f",
2258  log10(absint) , log10(CBabsint) , error );
2259  }
2260  else
2261  TotalInsanity();
2262  if( fabs(error) > fabs(AssertError[i]) )
2263  fprintf(ioMONITOR , " botch \n");
2264  else
2265  fprintf(ioMONITOR , "\n");
2266 
2267  PredQuan[i] = 0;
2268  AssertQuantity[i] = 0.;
2269  RelError[i] = MAX2( RelError[i] , fabs(error) );
2270 
2271  /* save sum which we will report later */
2272  MonitorResults.SumErrorCaseMonitor += RelError[i];
2273  ++MonitorResults.nSumErrorCaseMonitor;
2274 
2275  }
2276  }
2277  fprintf(ioMONITOR,"\n");
2278  }
2279  else if( nISOCaseB == ipHE_LIKE )
2280  {
2281  if( !dense.lgElmtOn[ipHELIUM] )
2282  {
2283  fprintf(ioQQQ,"PROBLEM monitor case B for a He is requested but He is not "
2284  "included.\n");
2285  fprintf(ioQQQ,"Do not turn off He if you want to monitor its spectrum.\n");
2286  cdEXIT(EXIT_FAILURE);
2287  }
2288 
2289  /* do He I as special case */
2290  fprintf(ioMONITOR," Wl Computed Asserted error\n");
2291  for( unsigned int ipLine=0; ipLine< atmdat.CaseBWlHeI.size(); ++ipLine )
2292  {
2293  /* monitor the line */
2294  realnum wl = atmdat.CaseBWlHeI[ipLine];
2295  /* range option to restrict wavelength coverage */
2296  if( wl < Param[i][2] || wl > Param[i][3] )
2297  continue;
2298  double relint , absint,CBrelint , CBabsint;
2299  cdLine( chAssertLineLabel[i].c_str(), wl, &CBrelint, &CBabsint, iLineType[i] );
2300  if( CBrelint < Param[i][4] )
2301  continue;
2302  double error;
2303  if( cdLine( chElemLabelCaseB.c_str(), wl, &relint, &absint, iLineType[i] ) > 0)
2304  {
2305  if (AssertError[i] > 0.0)
2306  error = (CBabsint - absint)/MAX2(CBabsint , absint);
2307  else
2308  error = (CBabsint - absint);
2309  double RelativeError = fabs(error / AssertError[i]);
2310  /* start of line, flag problems */
2311  if( RelativeError < 1. )
2312  {
2313  if( RelativeError < 0.25 )
2314  {
2315  fprintf( ioMONITOR, " ChkMonitor ");
2316  }
2317  else if( RelativeError < 0.50 )
2318  {
2319  fprintf( ioMONITOR, " ChkMonitor - ");
2320  }
2321  else if( RelativeError < 0.75 )
2322  {
2323  fprintf( ioMONITOR, " ChkMonitor -- ");
2324  }
2325  else if( RelativeError < 0.90 )
2326  {
2327  fprintf( ioMONITOR, " ChkMonitor --- ");
2328  }
2329  else if( RelativeError < 0.95 )
2330  {
2331  fprintf( ioMONITOR, " ChkMonitor ---- ");
2332  }
2333  else if( RelativeError < 0.98 )
2334  {
2335  fprintf( ioMONITOR, " ChkMonitor ----- ");
2336  }
2337  else
2338  {
2339  fprintf( ioMONITOR, " ChkMonitor ------ ");
2340  }
2341 
2342  }
2343  else
2344  {
2345  fprintf( ioMONITOR, " ChkMonitor botch>>");
2346  }
2347  prt_wl(ioMONITOR, wl );
2348  fprintf(ioMONITOR," %.2e %.2e %10.3f",
2349  absint , CBabsint , error );
2350  }
2351  else
2352  TotalInsanity();
2353  if( fabs(error) > fabs(AssertError[i]) )
2354  fprintf(ioMONITOR , " botch \n");
2355  else
2356  fprintf(ioMONITOR , "\n");
2357 
2358  PredQuan[i] = 0;
2359  AssertQuantity[i] = 0.;
2360  RelError[i] = MAX2( RelError[i] , fabs(error) );
2361 
2362  /* save sum which we will report later */
2363  MonitorResults.SumErrorCaseMonitor += RelError[i];
2364  ++MonitorResults.nSumErrorCaseMonitor;
2365  }
2366  fprintf(ioMONITOR,"\n");
2367  }
2368  else
2369  TotalInsanity();
2370  }
2371 
2372  // departure coefficients for a species given in quotes
2373  else if( strcmp( chAssertType[i] , "DC") == 0 )
2374  {
2375  string this_species = strAssertSpecies[i];
2376  if( this_species.find( '[' ) == string::npos )
2377  {
2378  // keyword EXCITED appeared
2379  if( AssertQuantity2[i] == 1 )
2380  this_species += "[2:]";
2381  else
2382  this_species += "[:]";
2383  }
2384 
2385  vector<long> speciesLevels;
2386  const molezone *sp = getLevelsGeneric( this_species.c_str(), false, speciesLevels );
2387  if( sp == null_molezone )
2388  {
2389  fprintf( ioQQQ, "PROBLEM Could not find species between quotes: \"%s\".\n",
2390  strAssertSpecies[i].c_str() );
2391  cdEXIT(EXIT_FAILURE);
2392  }
2393 
2394  if( sp->levels == NULL )
2395  {
2396  fprintf( ioQQQ, "WARNING Species '%s' has no internal structure."
2397  " Cannot compute departure coefficient\n",
2398  strAssertSpecies[i].c_str() );
2399  // Just report unity?
2400  PredQuan[i] = 1.;
2401  RelError[i] = 0.;
2402  }
2403  else if( speciesLevels.size() > 0 )
2404  {
2405  ASSERT( sp->levels->size() > 0 );
2406  RelError[i] = 0.;
2407  PredQuan[i] = 0.;
2408 
2409  long numPrintLevels = 0;
2410  for( vector<long>::const_iterator ilvl = speciesLevels.begin(); ilvl != speciesLevels.end(); ilvl++ )
2411  {
2412  const qStateConstProxy& st = (*sp->levels)[ *ilvl ];
2413  if( st.status() == LEVEL_INACTIVE )
2414  continue;
2415  ++numPrintLevels;
2416  PredQuan[i] += st.DepartCoef();
2417  }
2418  ASSERT( numPrintLevels > 0 );
2419  PredQuan[i] /= (double)(numPrintLevels);
2420  RelError[i] = AssertQuantity[i] - PredQuan[i];
2421 
2422 #if 0
2423  if( 0 && fp_equal( Param[i][1], 1. ) && PredQuan[i]==0. )
2424  {
2425  // this should only happen if either the present stage or the parent has zero density.
2426  ASSERT( dense.xIonDense[nelem][nelem-ipISO]==0. || dense.xIonDense[nelem][nelem+1-ipISO]==0. );
2427  PredQuan[i] = AssertQuantity[i];
2428  RelError[i] = 0.;
2429  }
2430 #endif
2431  }
2432  else
2433  {
2434  fprintf( ioQQQ, "WARNING Requested level(s) in '%s' do not exist\n",
2435  strAssertSpecies[i].c_str() );
2436  PredQuan[i] = 0.;
2437  RelError[i] = AssertQuantity[i];
2438  }
2439  }
2440 
2441  /* departure coefficients for something in isoelectronic sequences */
2442  else if( strcmp( chAssertType[i] , "DI") == 0 )
2443  {
2444  /* this is departure coefficient for XX-like ion given by wavelength */
2445  /* stored number was element number on C scale */
2446  long ipISO = (long)Param[i][0];
2447  long nelem = (long)wavelength[i];
2448  if( !dense.lgElmtOn[nelem] )
2449  {
2450  fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",nelem);
2451  PredQuan[i] = 0.;
2452  RelError[i] = 0.;
2453  }
2454  else
2455  {
2456  RelError[i] = 0.;
2457  PredQuan[i] = 0.;
2458  long numPrintLevels = iso_sp[ipISO][nelem].numLevels_local - (long)AssertQuantity2[i];
2459  for( long n=(long)AssertQuantity2[i]; n<numPrintLevels+(long)AssertQuantity2[i]; ++n )
2460  {
2461  PredQuan[i] += iso_sp[ipISO][nelem].st[n].DepartCoef();
2462  }
2463  ASSERT( numPrintLevels > 0 );
2464  PredQuan[i] /= (double)(numPrintLevels);
2465  RelError[i] = AssertQuantity[i] - PredQuan[i];
2466 
2467  if( fp_equal( Param[i][1], 1. ) && PredQuan[i]==0. )
2468  {
2469  // this should only happen if either the present stage or the parent has zero density.
2470  ASSERT( dense.xIonDense[nelem][nelem-ipISO]==0. || dense.xIonDense[nelem][nelem+1-ipISO]==0. );
2471  PredQuan[i] = AssertQuantity[i];
2472  RelError[i] = 0.;
2473  }
2474  }
2475  }
2476 
2477  /* this is H- departure coefficient */
2478  else if( strcmp( chAssertType[i] , "d-") == 0 )
2479  {
2480  PredQuan[i] = hmi.hmidep;
2481  RelError[i] = AssertQuantity[i] - hmi.hmidep;
2482  }
2483 
2484  /* this would be ionization fraction */
2485  else if( (strcmp( chAssertType[i] , "f ") == 0) ||
2486  (strcmp( chAssertType[i] , "F ") == 0) )
2487  {
2488  char chWeight[7];
2489  if( strcmp( chAssertType[i] , "F ") == 0 )
2490  {
2491  strcpy( chWeight , "VOLUME" );
2492  }
2493  else
2494  {
2495  /* this is default case, Fraction over radius */
2496  strcpy( chWeight , "RADIUS" );
2497  }
2498  /* get ionization fraction, returns false if could not find it */
2499  if( cdIonFrac(
2500  /* four char string, null terminated, giving the element name */
2501  chAssertLineLabel[i].c_str(),
2502  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
2503  (long)wavelength[i],
2504  /* will be fractional ionization */
2505  &relint,
2506  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2507  chWeight ,
2508  /* do not want extra factor of density */
2509  false) )
2510  {
2511  fprintf( ioMONITOR,
2512  " monitor error: lgCheckMonitors could not find a line with label %s %f \n",
2513  chAssertLineLabel[i].c_str() , wavelength[i] );
2514  /* go to next line */
2515  lg1OK[i] = false;
2516  RelError[i] = 0;
2517  PredQuan[i] = 0;
2518  lgFound[i] = false;
2519  continue;
2520  }
2521  /* this is ionization fraction */
2522  PredQuan[i] = relint;
2523  if (AssertError[i] > 0.)
2524  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2525  else
2526  RelError[i] = PredQuan[i]- AssertQuantity[i];
2527  }
2528 
2529  /* this would be column density of several molecules */
2530  else if( strcmp( chAssertType[i] , "cd") == 0)
2531  {
2532  /* for H2 column density - total or for a state? */
2533  if( ( chAssertLineLabel[i] == "H2" ) && (Param[i][0] >= 0.) )
2534  {
2535  /* this branch get state specific column density */
2536  /* get total H2 column density */
2537  if( (relint = cdH2_colden( (long)Param[i][0] , (long)Param[i][1] ) ) < 0. )
2538  {
2539  fprintf(ioQQQ," PROBLEM lgCheckMonitors did not find v=%li, J=%li for H2 column density.\n",
2540  (long)Param[i][0] , (long)Param[i][1] );
2541  lg1OK[i] = false;
2542  RelError[i] = 0;
2543  PredQuan[i] = 0;
2544  lgFound[i] = false;
2545  continue;
2546  }
2547  }
2548  else
2549  {
2550  /* get ionization fraction, returns 0 if all ok */
2551  if( cdColm(
2552  /* four char string, null terminated, giving the element name */
2553  chAssertLineLabel[i].c_str(),
2554  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number,
2555  * zero for molecule*/
2556  (long)wavelength[i],
2557  /* will be fractional ionization */
2558  &relint) )
2559  {
2560  fprintf( ioMONITOR,
2561  " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
2562  chAssertLineLabel[i].c_str() , wavelength[i] );
2563  /* go to next line */
2564  lg1OK[i] = false;
2565  RelError[i] = 0;
2566  PredQuan[i] = 0;
2567  lgFound[i] = false;
2568  continue;
2569  }
2570  }
2571  /* this is ionization fraction */
2572  PredQuan[i] = relint;
2573  if (AssertError[i] > 0.)
2574  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2575  else
2576  RelError[i] = PredQuan[i]- AssertQuantity[i];
2577  }
2578 
2579  /* this would be molecular fraction of CO or H2 */
2580  else if( (strcmp( chAssertType[i] , "MF") == 0) || (strcmp( chAssertType[i] , "mf") == 0) )
2581  {
2582  /* get molecular fraction, returns 0 if all ok */
2583  relint = 0.;
2584  if( chAssertLineLabel[i] == "H2" )
2585  {
2586  /* get total H2 column density */
2587  if( cdColm("H2 " , 0,
2588  /* will be fractional ionization */
2589  &relint) )
2590  TotalInsanity();
2591 
2592  relint = relint / colden.colden[ipCOL_HTOT];
2593  }
2594  else
2595  {
2596  fprintf( ioMONITOR,
2597  " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
2598  chAssertLineLabel[i].c_str() , wavelength[i] );
2599  /* go to next line */
2600  lg1OK[i] = false;
2601  RelError[i] = 0;
2602  PredQuan[i] = 0;
2603  continue;
2604  }
2605  /* this is ionization fraction */
2606  PredQuan[i] = relint;
2607  if (AssertError[i] > 0.)
2608  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2609  else
2610  RelError[i] = PredQuan[i]- AssertQuantity[i];
2611  }
2612 
2613  /* check heating/cooling at some temperature in a thermal map */
2614  else if( strcmp( chAssertType[i] , "mh") == 0 || strcmp( chAssertType[i] , "mc") == 0)
2615  {
2616  /* check heating or cooling (stored in error array) at temperature in monitor results */
2617  /* check that map was done, and arrays have nmap elements */
2618  if( hcmap.nMapAlloc == 0 )
2619  {
2620  /* this happens if map not done and space for h/c not allocated */
2621  fprintf( ioMONITOR,
2622  " monitor error: lgCheckMonitors cannot check map since map not done.\n");
2623  /* go to next line */
2624  lg1OK[i] = false;
2625  RelError[i] = 0;
2626  PredQuan[i] = 0;
2627  continue;
2628  }
2629  /* now check that requested temperature is within the range of the map we computed */
2630  if( AssertQuantity2[i]< hcmap.temap[0] || AssertQuantity2[i]> hcmap.temap[hcmap.nmap-1] )
2631  {
2632  fprintf( ioMONITOR,
2633  " monitor error: lgCheckMonitors cannot check map since temperature not within range.\n");
2634  /* go to next line */
2635  lg1OK[i] = false;
2636  RelError[i] = 0;
2637  PredQuan[i] = 0;
2638  continue;
2639  }
2640 
2641  /* we should have valid data - find closest temperature >- requested temperature */
2642  j = 0;
2643  while( AssertQuantity2[i]>hcmap.temap[j]*1.001 && j < hcmap.nmap )
2644  {
2645  ++j;
2646  }
2647 
2648  /* j points to correct cell in heating cooling array */
2649  /* we will not interpolate, just use this value, and clobber te to prove it*/
2650  if( strcmp( chAssertType[i] , "mh") == 0 )
2651  {
2652  /* heating */
2653  PredQuan[i] = hcmap.hmap[j];
2654  chAssertLineLabel[i] = "MapH" ;
2655  }
2656  else if( strcmp( chAssertType[i] , "mc") == 0)
2657  {
2658  /* cooling */
2659  PredQuan[i] = hcmap.cmap[j];
2660  chAssertLineLabel[i] = "MapC" ;
2661  }
2662  if (AssertError[i] > 0.)
2663  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2664  else
2665  RelError[i] = PredQuan[i]- AssertQuantity[i];
2666  }
2667 
2668  /* this will be an average temperature */
2669  else if( (strcmp( chAssertType[i] , "t ") == 0) ||
2670  (strcmp( chAssertType[i] , "T ") == 0) )
2671  {
2672  char chWeight[7];
2673  if( strcmp( chAssertType[i] , "T ") == 0 )
2674  {
2675  strcpy( chWeight , "VOLUME" );
2676  }
2677  else
2678  {
2679  /* this is default case, Fraction over radius */
2680  strcpy( chWeight , "RADIUS" );
2681  }
2682 
2683  /* options are average Te for ion, temp at ill face, or temp for grain */
2684  if( chAssertLineLabel[i] == "GTem" )
2685  {
2686  size_t nd;
2687  /* the minus one is because the grain types are counted from one,
2688  * but stuffed into the c array, that counts from zero */
2689  nd = (size_t)wavelength[i]-1;
2690  if( nd >= gv.bin.size() )
2691  {
2692  fprintf( ioQQQ, "Illegal grain number found: %f\n" , wavelength[i] );
2693  fprintf( ioQQQ, "Use 1 for first grain that is turned on, " );
2694  fprintf( ioQQQ, "2 for second, etc....\n" );
2695  fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" );
2696  cdEXIT(EXIT_FAILURE);
2697  }
2698  relint = gv.bin[nd]->avdust/radius.depth_x_fillfac;
2699  }
2700  else if( chAssertLineLabel[i] == "face" )
2701  {
2702  /* this is the temperature at the illuminated face */
2703  relint = struc.testr[0];
2704  }
2705  else
2706  {
2707  /* get temperature, returns false if could not find it */
2708  if( cdTemp(
2709  /* four char string, null terminated, giving the element name */
2710  chAssertLineLabel[i].c_str(),
2711  /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */
2712  (long)wavelength[i],
2713  /* will be mean temperatue */
2714  &relint,
2715  /* how to weight the average, must be "VOLUME" or "RADIUS" */
2716  chWeight ) )
2717  {
2718  fprintf( ioMONITOR,
2719  " monitor error: lgCheckMonitors could not find an ion with label %s ion %li \n",
2720  chAssertLineLabel[i].c_str() , (long)wavelength[i] );
2721  /* go to next line */
2722  lg1OK[i] = false;
2723  RelError[i] = 0;
2724  PredQuan[i] = 0;
2725  lgFound[i] = false;
2726  continue;
2727  }
2728  }
2729  /* this is the temperature */
2730  PredQuan[i] = relint;
2731  if (AssertError[i] > 0.)
2732  RelError[i] = get_error_ratio( PredQuan[i], AssertQuantity[i] );
2733  else
2734  RelError[i] = PredQuan[i]- AssertQuantity[i];
2735  }
2736 
2737  /* this would be grain potential in volt */
2738  else if( strcmp( chAssertType[i], "gp" ) == 0 )
2739  {
2740  /* the minus one is because the grain types are counted from one,
2741  * but stuffed into the c array, that counts from zero */
2742  size_t nd = (size_t)wavelength[i]-1;
2743  if( nd >= gv.bin.size() )
2744  {
2745  fprintf( ioQQQ, "Illegal grain number found: %g\n" , wavelength[i] );
2746  fprintf( ioQQQ, "Use 1 for first grain that is turned on, " );
2747  fprintf( ioQQQ, "2 for second, etc....\n" );
2748  fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" );
2749  cdEXIT(EXIT_FAILURE);
2750  }
2751 
2752  /* get average grain potential in volt, always averaged over radius */
2753  PredQuan[i] = gv.bin[nd]->avdpot/radius.depth_x_fillfac;
2754  /* actually absolute error, potential can be zero! */
2755  RelError[i] = AssertQuantity[i] - PredQuan[i];
2756  }
2757 
2758  else
2759  {
2760  fprintf( ioMONITOR,
2761  " monitor error: lgCheckMonitors received an insane chAssertType=%s, impossible\n",
2762  chAssertType[i] );
2763  ShowMe();
2764  cdEXIT(EXIT_FAILURE);
2765  }
2766 
2767  if( chAssertLimit[i] == '=' )
2768  {
2769  /* predicted quantity should be within error of expected */
2770  if( fabs(RelError[i]) > fabs(AssertError[i]) )
2771  {
2772  lg1OK[i] = false;
2773  lgMonitorsOK = false;
2774  }
2775  }
2776  else if( chAssertLimit[i] == '<' )
2777  {
2778  /* expected is an upper limit, so PredQuan/AssertQuantity should
2779  * be less than one, and so RelError =1-PredQuan[i]/AssertQuantity[i]
2780  * should be >= 0
2781  * in case of iterations or zones, iter < iterations should
2782  * trigger botched monitor when iter == iterations */
2783  if( RelError[i] <= 0. )
2784  {
2785  lg1OK[i] = false;
2786  lgMonitorsOK = false;
2787  }
2788  }
2789  else if( chAssertLimit[i] == '>' )
2790  {
2791  /* expected is a lower limit, so PredQuan/AssertQuantity should
2792  * be greater than one, and so RelError should be negative */
2793  if( RelError[i] >= 0. )
2794  {
2795  lg1OK[i] = false;
2796  lgMonitorsOK = false;
2797  }
2798  }
2799  }
2800 
2801  /* only print summary if we are talking */
2802  if( called.lgTalk && nAsserts>0 )
2803  {
2804  time_t now;
2805 
2806  /* First disambiguate any line identifications */
2807  if( lgDisambiguate )
2808  {
2809  /* change significant figures of WL for this printout */
2810  long sigfigsav = LineSave.sig_figs;
2811  double relint1, relint2, absint1;
2812 
2813  LineSave.sig_figs = LineSave.sig_figs_max;
2814 
2815  fprintf( ioMONITOR, "=============Line Disambiguation=======================================================\n" );
2816  fprintf( ioMONITOR, " Wavelengths || Intensities \n" );
2817  fprintf( ioMONITOR, "Label line match1 match2 match3 || asserted match1 match2 match3\n" );
2818 
2819  for( i=0; i<nAsserts; ++i )
2820  {
2821  if( ipDisambiguate[i][1] > 0 )
2822  {
2823  fprintf( ioMONITOR , "%-*s ", NCHLAB-1, chAssertLineLabel[i].c_str() );
2824  prt_wl( ioMONITOR , wavelength[i] );
2825  fprintf( ioMONITOR , " " );
2826  prt_wl( ioMONITOR , LineSave.lines[ipDisambiguate[i][0]].wavelength() );
2827  fprintf( ioMONITOR , " " );
2828  prt_wl( ioMONITOR , LineSave.lines[ipDisambiguate[i][1]].wavelength() );
2829  fprintf( ioMONITOR , " " );
2830  if( ipDisambiguate[i][2] > 0 )
2831  {
2832  prt_wl( ioMONITOR , LineSave.lines[ipDisambiguate[i][2]].wavelength() );
2833  cdLine_ip( ipDisambiguate[i][2], &relint2, &absint1, iLineType[i] );
2834  }
2835  else
2836  {
2837  fprintf( ioMONITOR , "--------" );
2838  relint2 = 0.0;
2839  }
2840  fprintf( ioMONITOR , " ||" );
2841 
2842  cdLine_ip( ipDisambiguate[i][1], &relint1, &absint1, iLineType[i] );
2843 
2844  if( lgPrtSciNot )
2845  {
2846  fprintf( ioMONITOR , " %10.3e %10.3e %10.3e %10.3e\n",
2847  AssertQuantity[i],
2848  PredQuan[i] ,
2849  relint1,
2850  relint2 );
2851  }
2852  else
2853  {
2854  fprintf( ioMONITOR , " %10.4f %10.4f %10.4f %10.4f\n",
2855  AssertQuantity[i],
2856  PredQuan[i] ,
2857  relint1,
2858  relint2 );
2859  }
2860  }
2861  }
2862  fprintf( ioMONITOR, "\n" );
2863 
2864  /* revert to original significant figures */
2865  LineSave.sig_figs = sigfigsav;
2866  }
2867 
2868  /* write start of title and version number of code */
2869  fprintf( ioMONITOR, "=============Results of monitors: Cloudy %s ", t_version::Inst().chVersion );
2870 
2871  /* usually print date and time info - do not if "no times" command entered,
2872  * which set this flag false */
2873  if( prt.lgPrintTime )
2874  {
2875  /* now add date of this run */
2876  now = time(NULL);
2877  /* now print this time at the end of the string. the system put cr at the end of the string */
2878  fprintf(ioMONITOR,"%s", ctime(&now) );
2879  }
2880  else
2881  {
2882  fprintf(ioMONITOR,"\n" );
2883  }
2884 
2885  if( lgMonitorsOK )
2886  {
2887  fprintf( ioMONITOR, " No errors were found. Summary follows.\n");
2888  }
2889  else
2890  {
2891  fprintf( ioMONITOR, " Errors were found. Summary follows.\n");
2892  }
2893 
2894  fprintf( ioMONITOR,
2895  " %-*s%*s computed asserted Rel Err Set err type \n",
2896  NCHLAB, "Label", LineSave.wl_length, "line" );
2897  /* now print a summary */
2898  for( i=0; i<nAsserts; ++i )
2899  {
2900  PrtOneMonitor( ioMONITOR, chAssertType[i], chAssertLineLabel[i].c_str(),
2901  wavelength[i], iLineType[i], PredQuan[i], chAssertLimit[i], AssertQuantity[i],
2902  RelError[i], AssertError[i], lg1OK[i], lgQuantityLog[i], lgFound[i] );
2903  }
2904  fprintf( ioMONITOR , " \n");
2905 
2906  /* NB - in following, perl scripts detect these strings - be careful if they
2907  * are changed - the scripts may no longer detect major errors */
2908  if( !lgMonitorsOK && lgBigBotch )
2909  {
2910  /* there were big botches */
2911  fprintf( ioMONITOR, " BIG BOTCHED MONITORS!!! Big Botched Monitors!!! \n");
2912  }
2913  else if( !lgMonitorsOK )
2914  {
2915  fprintf( ioMONITOR, " BOTCHED MONITORS!!! Botched Monitors!!! \n");
2916  }
2917 
2918  if( MonitorResults.nSumErrorCaseMonitor>0 )
2919  {
2920  fprintf(ioMONITOR,"\n The mean of the %li monitor Case A, B relative "
2921  "residuals is %.2f\n\n" ,
2922  MonitorResults.nSumErrorCaseMonitor,
2923  MonitorResults.SumErrorCaseMonitor /MonitorResults.nSumErrorCaseMonitor );
2924  }
2925 
2926  /* explain how we were compiled, but only if printing time */
2927  if( prt.lgPrintTime )
2928  {
2929  fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo );
2930  }
2931  }
2932  return lgMonitorsOK;
2933 }
2934 
2935 STATIC void prtLineType( FILE *ioMONITOR, const int iLineType )
2936 {
2937  switch( iLineType )
2938  {
2939  case -1:
2940  /* The monitor is not an atomic transition. */
2941  fprintf( ioMONITOR, " " );
2942  break;
2943  case 0:
2944  fprintf( ioMONITOR, " intr " );
2945  break;
2946  case 1:
2947  fprintf( ioMONITOR, " emer " );
2948  break;
2949  case 2:
2950  fprintf( ioMONITOR, " intr cumu" );
2951  break;
2952  case 3:
2953  fprintf( ioMONITOR, " emer cumu" );
2954  break;
2955  default:
2956  fprintf( ioMONITOR, "ERROR: Unrecognized line type: %d\n",
2957  iLineType );
2958  cdEXIT( EXIT_FAILURE );
2959  break;
2960  }
2961 }
2962 
2963 void PrtOneMonitor( FILE *ioMONITOR, const char* chAssertType, const string chAssertLineLabel,
2964  const realnum wavelength, const int iLineType, const double PredQuan, const char chAssertLimit, const double AssertQuantity,
2965  const double RelError, const double AssertError, const bool lg1OK, const bool lgQuantityLog, const bool lgFound )
2966 {
2967  DEBUG_ENTRY( "PrtOneMonitor()" );
2968 
2969  double prtPredQuan, prtAssertQuantity;
2970  // this is option to print log of quantity rather than linear. default is
2971  // linear, and log only for a few such as ionization to see small numbers
2972  if( lgQuantityLog )
2973  {
2974  if( PredQuan > 0. )
2975  prtPredQuan = log10( MAX2( SMALLDOUBLE , PredQuan ) );
2976  else
2977  prtPredQuan = -37.;
2978  prtAssertQuantity = log10( MAX2( SMALLDOUBLE , AssertQuantity ) );
2979  }
2980  else
2981  {
2982  prtPredQuan = PredQuan;
2983  prtAssertQuantity = AssertQuantity;
2984  }
2985 
2986  /* start of line, flag problems */
2987  if( lg1OK )
2988  {
2989  // the ChkMonitor is a unique label so that we can grep on all output
2990  // and see what happened, without picking up input stream
2991  double relative = fabs( RelError / SDIV( fabs(AssertError)));
2992 
2993  if( relative < 0.25 || chAssertLimit != '=' )
2994  {
2995  fprintf( ioMONITOR, " ChkMonitor ");
2996  }
2997  else if( relative < 0.50 )
2998  {
2999  fprintf( ioMONITOR, " ChkMonitor - ");
3000  }
3001  else if( relative < 0.75 )
3002  {
3003  fprintf( ioMONITOR, " ChkMonitor -- ");
3004  }
3005  else if( relative < 0.90 )
3006  {
3007  fprintf( ioMONITOR, " ChkMonitor --- ");
3008  }
3009  else if( relative < 0.95 )
3010  {
3011  fprintf( ioMONITOR, " ChkMonitor ---- ");
3012  }
3013  else if( relative < 0.98 )
3014  {
3015  fprintf( ioMONITOR, " ChkMonitor ----- ");
3016  }
3017  else
3018  {
3019  fprintf( ioMONITOR, " ChkMonitor ------ ");
3020  }
3021  }
3022  else
3023  {
3024  fprintf( ioMONITOR, " ChkMonitor botch>>");
3025  }
3026 
3027  fprintf( ioMONITOR , "%-*s ", NCHLAB-1, chAssertLineLabel.c_str() );
3028 
3029  /* special formatting for the emission lines */
3030  if( strcmp( chAssertType, "Ll" )==0 || strcmp( chAssertType, "Lr" )==0 ||
3031  strcmp( chAssertType, "Lb" )==0 )
3032  {
3033  prt_wl( ioMONITOR , wavelength );
3034  }
3035  else
3036  {
3037  fprintf( ioMONITOR , "%*i", LineSave.wl_length, (int)wavelength );
3038  }
3039 
3040  const char* format = " %10.4f %c %10.4f %7.3f %7.3f ";
3041  if( lgPrtSciNot )
3042  format = " %10.3e %c %10.3e %7.3f %7.3f ";
3043 
3044  fprintf( ioMONITOR, format,
3045  prtPredQuan,
3046  chAssertLimit,
3047  prtAssertQuantity,
3048  RelError,
3049  AssertError);
3050 
3051  prtLineType( ioMONITOR, iLineType );
3052 
3053  // if botched and the botch is > 3 sigma, say BIG BOTCH,
3054  // the lg1OK is needed since some tests (number of zones, etc)
3055  // are limits, not the quantity, and if limit is large the
3056  // miss will be big too
3057 
3058  /* >>chng 02 nov 27, added lgFound so don't get big botch when line simply missing */
3059  if( !lg1OK && (fabs(RelError) > 3.*AssertError) && lgFound )
3060  {
3061  fprintf( ioMONITOR , " <<BIG BOTCH!!\n");
3062  lgBigBotch = true;
3063  }
3064  else
3065  {
3066  fprintf( ioMONITOR , "\n");
3067  }
3068 
3069  return;
3070 }
monitor_results.h
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:140
cdLine_ip
void cdLine_ip(long int ipLine, double *relint, double *absint)
Definition: cddrive.cpp:1112
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:230
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:438
WavlenErrorGet
realnum WavlenErrorGet(realnum wavelength, long sig_figs)
Definition: lines_service.cpp:196
exp10
double exp10(double x)
Definition: cddefines.h:1383
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:1113
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
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:1657
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:415
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:912
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:184
chAssertLimit
static vector< char > chAssertLimit
Definition: monitor_results.cpp:64
LineSave
t_LineSave LineSave
Definition: lines.cpp:9
Parser::nMatchErase
bool nMatchErase(const char *chKey)
Definition: parser.h:163
cdTemp
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition: cddrive.cpp:1296
qStateConstProxy::DepartCoef
double DepartCoef() const
Definition: quantumstate.h:338
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
Parser::getWave
double getWave()
Definition: parser.cpp:345
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:807
Parser
Definition: parser.h:42
cap4
void cap4(char *chCAP, const char *chLab)
Definition: service.cpp:264
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:155
Parser::StringError
NORETURN void StringError() const
Definition: parser.cpp:174
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:79
t_rfield::EnergyIncidCont
realnum EnergyIncidCont
Definition: rfield.h:465
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:858
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:556
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:2312
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:983
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:202
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:2935
qList::size
size_t size() const
Definition: quantumstate.h:131
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:484
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:318
ErrorDefault
static realnum ErrorDefault
Definition: monitor_results.cpp:50
Parser::GetElem
long int GetElem(void) const
Definition: parser.cpp:294
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:1015
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:400
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:617
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:2963
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1069
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:595
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:144
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
powpq
double powpq(double x, int p, int q)
Definition: service.cpp:726
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:349
generic_state.h
t_cpu_i::lgMPI
bool lgMPI() const
Definition: cpu.h:388
wavelength
static vector< realnum > wavelength
Definition: monitor_results.cpp:76
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:585
AssertQuantity2
static vector< double > AssertQuantity2
Definition: monitor_results.cpp:70
t_rfield::EnergyDiffCont
realnum EnergyDiffCont
Definition: rfield.h:465
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:103
NCHLAB
const int NCHLAB
Definition: cddefines.h:303
MAX2
#define MAX2(a, b)
Definition: cddefines.h:828
lgCheckMonitors
bool lgCheckMonitors(FILE *ioMONITOR)
Definition: monitor_results.cpp:1664
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
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:1006
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:200
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:304
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:423
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:153
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:348
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:1681
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:576
t_struc::pressure
realnum * pressure
Definition: struc.h:25
t_thermal::ctot
double ctot
Definition: thermal.h:130
Generated by   doxygen 1.8.5