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 
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 */
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 */
51 
52 /* default relative error for monitored performance metrics */
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;
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*/
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
112 }
113 
114 /*======================================================================*/
115 /*ParseMonitorResults parse the monitor command */
117 {
118  long nelem,
119  n2;
120 
121  DEBUG_ENTRY( "ParseMonitorResults()" );
122 
123  if( !lgInitDone )
124  {
125  fprintf( ioQQQ, " ParseMonitorResults called before InitAsserResults\n" );
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
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" );
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" );
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 */
300 
301  /* now get ionization stage, which will be saved into wavelength */
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" );
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.*/
325  if( p.lgEOL() )
326  p.NoNumb("ionization fraction");
327  }
328 
329  /* optional error, default available */
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 */
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" );
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" );
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" );
397  }
398 
399  /* not meaningful */
400  wavelength[nAsserts] = 0;
401 
402  /* now get log of molecular fraction */
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 */
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 */
418  if( p.lgEOL() )
419  {
420  /* default error was set above */
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" );
485  }
486 
487  /* copy string into array */
488  chAssertLineLabel[nAsserts] = chLabel;
489 
490  /* now get line wavelength */
492 
493  /* now get intensity or luminosity -
494  * rel intensity is linear and intensity or luminosity are log */
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" );
511  }
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" );
521  }
522 
523  /* optional error, default available */
525  if( p.lgEOL() )
526  {
527  /* default error was set above */
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 */
539  if( p.lgEOL() )
540  /* default error was set above */
541  AssertError[nAsserts] = ErrorDefault;
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");
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");
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);
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);
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" );
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 */
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 */
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" );
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" );
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" );
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 */
734  }
735  else
736  {
737  /* do the ground state */
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" );
762  }
763 
764  /* now get temperature for AssertQuantity2 array*/
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 */
775  exp10(AssertQuantity2[nAsserts] );
776  }
777 
778  /* print the temperature in the wavelength column */
780 
781  /* heating or cooling, both log, put into error */
783  if( p.lgEOL() )
784  {
785  p.NoNumb("heating/cooling");
786  }
787 
788  /* AssertQuantity array will have heating or cooling */
790 
791  /* optional error, default available (cannot do before loop since we
792  * do not know how many numbers are on line */
794  if( p.lgEOL() )
795  {
796  /* default error was set above */
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" );
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" );
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 */
929  if( p.lgEOL() )
930  {
931  p.NoNumb("column density");
932  }
933  /* entered as log, but we will compare with linear */
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 */
940  if( p.lgEOL() )
941  {
942  /* default error was set above */
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" );
974  }
975  /* now past the 2 in h2, get the real number */
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. )
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 */
990  if( p.lgEOL() )
991  {
992  /* default error was set above */
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 */
1009  /* the potential itself, in volt, always linear */
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 */
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 */
1085  }
1086  else
1087  {
1088  fprintf( ioQQQ,
1089  " I could not identify an element name on this line.\n");
1090  fprintf( ioQQQ, " Sorry.\n" );
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  {
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" );
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 */
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 */
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*/
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 */
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" );
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  }
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 */
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.;
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.;
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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  {
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 */
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 */
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() )
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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 */
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.;
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" );
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 );
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 */
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;
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 
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 */
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;
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");
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 */
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() );
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 */
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" );
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" );
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();
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 
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 }
bool nMatch(const char *chKey) const
Definition: parser.h:140
void cdLine_ip(long int ipLine, double *relint, double *absint)
Definition: cddrive.cpp:1112
double Radius
Definition: radius.h:31
realnum WaveLengthCaseB[8][25][24]
Definition: atmdat.h:355
double depth
Definition: radius.h:31
double htot
Definition: thermal.h:169
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
t_atmdat atmdat
Definition: atmdat.cpp:6
static bool lgSpaceAllocated
void InitMonitorResults(void)
t_thermal thermal
Definition: thermal.cpp:6
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:230
t_colden colden
Definition: colden.cpp:5
qList st
Definition: iso.h:482
double FFmtRead(void)
Definition: parser.cpp:438
realnum WavlenErrorGet(realnum wavelength, long sig_figs)
double exp10(double x)
Definition: cddefines.h:1383
const int ipHE_LIKE
Definition: iso.h:65
double ForcePass(char chAssertLimit1)
static vector< bool > lgQuantityLog
map< std::string, std::vector< TransitionProxy > >::iterator blend_iterator
Definition: transition.h:668
T * get_ptr(T *v)
Definition: cddefines.h:1113
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double get_error_ratio(double pred, double assert)
t_struc struc
Definition: struc.cpp:6
bool lgRelease
Definition: version.h:28
const realnum SMALLFLOAT
Definition: cpu.h:246
t_monitorresults MonitorResults
t_cpu_i & i()
Definition: cpu.h:415
double * temap
Definition: hcmap.h:42
long int nSumErrorCaseMonitor
void ParseMonitorResults(Parser &p)
const double SMALLDOUBLE
Definition: cpu.h:250
int cdIonFrac(const char *chLabel, long int IonStage, double *fracin, const char *chWeight, bool lgDensity)
Definition: cddrive.cpp:912
t_conv conv
Definition: conv.cpp:5
long int nOptimiz
Definition: optimize.h:250
long int nTotalIoniz_start
Definition: conv.h:164
int GetQuote(string &chLabel)
Definition: parser.cpp:184
static vector< char > chAssertLimit
t_LineSave LineSave
Definition: lines.cpp:9
bool nMatchErase(const char *chKey)
Definition: parser.h:163
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition: cddrive.cpp:1296
double DepartCoef() const
Definition: quantumstate.h:338
FILE * ioQQQ
Definition: cddefines.cpp:7
double getWave()
Definition: parser.cpp:345
long int nzone
Definition: cddefines.cpp:14
bool lgTalk
Definition: called.h:12
Definition: mole.h:378
static realnum ErrorDefaultPerformance
static multi_arr< double, 2 > Param
#define MIN2(a, b)
Definition: cddefines.h:807
Definition: parser.h:42
void cap4(char *chCAP, const char *chLab)
Definition: service.cpp:264
static vector< string > strAssertSpecies
int wl_length
Definition: lines.h:114
vector< LinSv > lines
Definition: lines.h:132
void trimTrailingWhiteSpace(string &str)
Definition: service.cpp:155
NORETURN void StringError() const
Definition: parser.cpp:174
t_dense dense
Definition: global.cpp:15
static t_version & Inst()
Definition: cddefines.h:209
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
void prt_line_err(FILE *ioOUT, const char *label, realnum wvlng)
Definition: prt.cpp:161
double depth_x_fillfac
Definition: radius.h:79
realnum EnergyIncidCont
Definition: rfield.h:465
Wind wind
Definition: wind.cpp:5
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
long int iteration
Definition: cddefines.cpp:16
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
#define MALLOC(exp)
Definition: cddefines.h:556
double ortho_density
Definition: h2_priv.h:326
double SumErrorCaseMonitor
double cdH2_colden(long iVib, long iRot)
Definition: mole_h2.cpp:2312
qList * levels
Definition: mole.h:417
long int n_HighestResolved_max
Definition: iso.h:536
bool lgMonitorsOK
static vector< double > AssertError
double para_density
Definition: h2_priv.h:326
#define POW2
Definition: cddefines.h:983
#define STATIC
Definition: cddefines.h:118
static long nAsserts
static vector< double > AssertQuantity
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
void input_readvector(const char *chFile, vector< double > &vec, bool *lgError)
Definition: input.cpp:202
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
STATIC void prtLineType(FILE *ioMONITOR, const int iLineType)
size_t size() const
Definition: quantumstate.h:131
#define cdEXIT(FAIL)
Definition: cddefines.h:484
realnum RadBetaMax
Definition: pressure.h:96
const molezone * getLevelsGeneric(const char *chLabel, bool lgValidate, vector< long > &LevelList)
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:318
static realnum ErrorDefault
long int GetElem(void) const
Definition: parser.cpp:294
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
static vector< std::vector< TransitionProxy > * > assertBlends
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1015
bool lgPrintTime
Definition: prt.h:161
long int sig_figs
Definition: lines.h:109
t_optimize optimize
Definition: optimize.cpp:6
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
t_radius radius
Definition: radius.cpp:5
t_prt prt
Definition: prt.cpp:14
long int nTotalIoniz
Definition: conv.h:159
long int nmap
Definition: hcmap.h:39
bool lgElmtOn[LIMELM]
Definition: dense.h:160
static vector< string > chAssertLineLabel
enum level_status status() const
Definition: quantumstate.h:400
#define ASSERT(exp)
Definition: cddefines.h:617
bool lgPrtSciNot
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)
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1069
void reserve(size_type i1)
const int ipH_LIKE
Definition: iso.h:64
double * hmap
Definition: hcmap.h:45
double hmidep
Definition: hmi.h:42
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition: cddrive.cpp:595
bool lgOptimize
Definition: optimize.h:257
bool lgBigBotch
bool GetParam(const char *chKey, double *val)
Definition: parser.h:144
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double powpq(double x, int p, int q)
Definition: service.cpp:726
const int ipHELIUM
Definition: cddefines.h:349
bool lgMPI() const
Definition: cpu.h:388
static vector< realnum > wavelength
vector< GrainBin * > bin
Definition: grainvar.h:585
static vector< double > AssertQuantity2
realnum EnergyDiffCont
Definition: rfield.h:465
double eden
Definition: dense.h:201
static bool lgInitDone
bool lgEOL(void) const
Definition: parser.h:103
const int NCHLAB
Definition: cddefines.h:303
#define MAX2(a, b)
Definition: cddefines.h:828
bool lgCheckMonitors(FILE *ioMONITOR)
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
map< std::string, std::vector< TransitionProxy > > blends
Definition: transition.cpp:36
realnum ** csupra
Definition: secondaries.h:33
double * cmap
Definition: hcmap.h:48
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
vector< realnum > CaseBWlHeI
Definition: atmdat.h:358
t_hcmap hcmap
Definition: hcmap.cpp:23
bool lgReleaseBranch
Definition: version.h:25
int PrintLine(FILE *fp) const
Definition: parser.h:200
long int nsum
Definition: lines.h:87
realnum * testr
Definition: struc.h:25
static const int NASSERTS
static const long sig_figs_max
Definition: lines.h:110
void caps(char *chCard)
Definition: service.cpp:304
GrainVar gv
Definition: grainvar.cpp:5
t_hmi hmi
Definition: hmi.cpp:5
t_secondaries secondaries
Definition: secondaries.cpp:5
static t_cpu cpu
Definition: cpu.h:423
void ShowMe(void)
Definition: service.cpp:205
static vector< int > iLineType
bool GetRange(const char *chKey, double *val1, double *val2)
Definition: parser.h:153
const int ipHYDROGEN
Definition: cddefines.h:348
realnum colden[NCOLD]
Definition: colden.h:32
static char ** chAssertType
molezone * null_molezone
long int numLevels_local
Definition: iso.h:529
realnum windv
Definition: wind.h:18
int lines_table()
static long int * ipLine
Definition: prt_linesum.cpp:14
t_called called
Definition: called.cpp:4
long int nMapAlloc
Definition: hcmap.h:53
double rate_h2_form_grains_used_total
Definition: grainvar.h:576
realnum * pressure
Definition: struc.h:25
double ctot
Definition: thermal.h:130