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