cloudy trunk
Loading...
Searching...
No Matches
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
37bool lgMonitorsOK , lgBigBotch, lgPrtSciNot;
38t_monitorresults MonitorResults;
39
40/* flag to remember that InitMonitorResults was called */
41static 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 */
46static const int NASSERTS = 200;
47
48/* default relative error for monitored physical quantities */
49static realnum ErrorDefault;
50
51/* default relative error for monitored performance metrics */
52static realnum ErrorDefaultPerformance;
53
54/* dim of 5 also appears in MALLOC below */
55static const int NCHAR = 5;
56static char **chAssertLineLabel;
57
58/* this will be = for equal, < or > for limit */
59static char *chAssertLimit;
60
61/* this will be a two character label identifying which type of monitor */
62static char **chAssertType;
63
64/* the values and error in the asserted quantity */
65static double *AssertQuantity ,*AssertQuantity2 ,*AssertError, **Param;
66
67/* this flag says where we print linear or log quantity */
68static int *lgQuantityLog;
69static long nAsserts=0;
70static realnum *wavelength;
71
72inline 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*/
87void 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 */
104void 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 */
1603bool 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}
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
AssertFeIIDep
void AssertFeIIDep(double *pred, double *BigError, double *StdDev)
Definition: atom_feii.cpp:2452
called
t_called called
Definition: called.cpp:5
nzone
long int nzone
Definition: cddefines.cpp:14
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
iteration
long int iteration
Definition: cddefines.cpp:16
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
MIN2
#define MIN2
Definition: cddefines.h:761
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
POW2
#define POW2
Definition: cddefines.h:929
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
cap4
void cap4(char *chCAP, const char *chLab)
Definition: service.cpp:240
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
realnum
float realnum
Definition: cddefines.h:103
MAX2
#define MAX2
Definition: cddefines.h:782
get_ptr
T * get_ptr(T *v)
Definition: cddefines.h:1079
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ShowMe
void ShowMe(void)
Definition: service.cpp:181
cdColm
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition: cddrive.cpp:636
cdLine_ip
void cdLine_ip(long int ipLine, double *relint, double *absint)
Definition: cddrive.cpp:1414
cdIonFrac
int cdIonFrac(const char *chLabel, long int IonStage, double *fracin, const char *chWeight, bool lgDensity)
Definition: cddrive.cpp:1072
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1228
cdTemp
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
Definition: cddrive.cpp:1602
cdH2_colden
double cdH2_colden(long iVib, long iRot)
Definition: mole_h2.cpp:2318
LineSv
LinSv * LineSv
Definition: cdinit.cpp:70
GrainVar::rate_h2_form_grains_used_total
double rate_h2_form_grains_used_total
Definition: grainvar.h:574
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
Parser
Definition: parser.h:32
Parser::GetElem
long int GetElem(void) const
Definition: parser.cpp:209
Parser::FFmtRead
double FFmtRead(void)
Definition: parser.cpp:353
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
Parser::nMatchErase
bool nMatchErase(const char *chKey)
Definition: parser.h:158
Parser::GetParam
bool GetParam(const char *chKey, double *val)
Definition: parser.h:139
Parser::GetRange
bool GetRange(const char *chKey, double *val1, double *val2)
Definition: parser.h:148
Parser::lgEOL
bool lgEOL(void) const
Definition: parser.h:98
Parser::NoNumb
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:233
Parser::getWave
double getWave()
Definition: parser.cpp:260
Parser::GetQuote
int GetQuote(char *chLabel, bool lgABORT)
Definition: parser.h:209
Parser::PrintLine
int PrintLine(FILE *fp) const
Definition: parser.h:204
Singleton< t_version >::Inst
static t_version & Inst()
Definition: cddefines.h:175
diatomics::para_density
double para_density
Definition: h2_priv.h:321
diatomics::ortho_density
double ortho_density
Definition: h2_priv.h:319
multi_arr
Definition: container_classes.h:942
size_t
t_cpu_i::lgMPI
bool lgMPI() const
Definition: cpu.h:322
t_cpu::i
t_cpu_i & i()
Definition: cpu.h:347
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
t_version::lgRelease
bool lgRelease
Definition: version.h:28
t_version::lgReleaseBranch
bool lgReleaseBranch
Definition: version.h:25
colden
t_colden colden
Definition: colden.cpp:5
ipCOL_HTOT
#define ipCOL_HTOT
Definition: colden.h:12
conv
t_conv conv
Definition: conv.cpp:5
cpu
static t_cpu cpu
Definition: cpu.h:355
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:195
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
dense
t_dense dense
Definition: dense.cpp:24
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
gv
GrainVar gv
Definition: grainvar.cpp:5
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
hcmap
t_hcmap hcmap
Definition: hcmap.cpp:21
hmi
t_hmi hmi
Definition: hmi.cpp:5
input_readvector
void input_readvector(const char *chFile, double vector[], long n, bool *lgEOF)
Definition: input.cpp:189
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
LineSave
t_LineSave LineSave
Definition: lines.cpp:5
lines_table
int lines_table()
Definition: parse_table.cpp:2042
WavlenErrorGet
realnum WavlenErrorGet(realnum wavelength)
Definition: lines_service.cpp:182
lines_service.h
lgSpaceAllocated
static bool lgSpaceAllocated
Definition: monitor_results.cpp:43
lgQuantityLog
static int * lgQuantityLog
Definition: monitor_results.cpp:68
NASSERTS
static const int NASSERTS
Definition: monitor_results.cpp:46
ForcePass
double ForcePass(char chAssertLimit1)
Definition: monitor_results.cpp:72
ErrorDefault
static realnum ErrorDefault
Definition: monitor_results.cpp:49
chAssertLineLabel
static char ** chAssertLineLabel
Definition: monitor_results.cpp:56
NCHAR
static const int NCHAR
Definition: monitor_results.cpp:55
lgPrtSciNot
bool lgPrtSciNot
Definition: monitor_results.cpp:37
chAssertType
static char ** chAssertType
Definition: monitor_results.cpp:62
lgInitDone
static bool lgInitDone
Definition: monitor_results.cpp:41
MonitorResults
t_monitorresults MonitorResults
Definition: monitor_results.cpp:38
wavelength
static realnum * wavelength
Definition: monitor_results.cpp:70
AssertQuantity2
static double * AssertQuantity2
Definition: monitor_results.cpp:65
AssertQuantity
static double * AssertQuantity
Definition: monitor_results.cpp:65
nAsserts
static long nAsserts
Definition: monitor_results.cpp:69
lgBigBotch
bool lgBigBotch
Definition: monitor_results.cpp:37
Param
static double ** Param
Definition: monitor_results.cpp:65
ErrorDefaultPerformance
static realnum ErrorDefaultPerformance
Definition: monitor_results.cpp:52
AssertError
static double * AssertError
Definition: monitor_results.cpp:65
lgMonitorsOK
bool lgMonitorsOK
Definition: monitor_results.cpp:37
InitMonitorResults
void InitMonitorResults(void)
Definition: monitor_results.cpp:87
chAssertLimit
static char * chAssertLimit
Definition: monitor_results.cpp:59
lgCheckMonitors
bool lgCheckMonitors(FILE *ioMONITOR)
Definition: monitor_results.cpp:1603
ParseMonitorResults
void ParseMonitorResults(Parser &p)
Definition: monitor_results.cpp:104
monitor_results.h
optimize
t_optimize optimize
Definition: optimize.cpp:5
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
SPEEDLIGHT
UNUSED const double SPEEDLIGHT
Definition: physconst.h:100
pressure
t_pressure pressure
Definition: pressure.cpp:5
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
prt
t_prt prt
Definition: prt.cpp:10
ipLine
static long int * ipLine
Definition: prt_linesum.cpp:15
radius
t_radius radius
Definition: radius.cpp:5
rfield
t_rfield rfield
Definition: rfield.cpp:8
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
struc
t_struc struc
Definition: struc.cpp:6
Wind::windv
realnum windv
Definition: wind.h:18
t_LineSave::sig_figs
long int sig_figs
Definition: lines.h:91
t_LineSave::nsum
long int nsum
Definition: lines.h:62
t_atmdat::CaseBWlHeI
vector< realnum > CaseBWlHeI
Definition: atmdat.h:209
t_atmdat::WaveLengthCaseB
realnum WaveLengthCaseB[8][25][24]
Definition: atmdat.h:206
t_called::lgTalk
bool lgTalk
Definition: called.h:12
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:38
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:166
t_conv::nTotalIoniz_start
long int nTotalIoniz_start
Definition: conv.h:171
t_dense::eden
double eden
Definition: dense.h:190
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_elementnames::chElementNameShort
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
t_hcmap::nMapAlloc
long int nMapAlloc
Definition: hcmap.h:53
t_hcmap::hmap
realnum * hmap
Definition: hcmap.h:45
t_hcmap::cmap
realnum * cmap
Definition: hcmap.h:48
t_hcmap::temap
realnum * temap
Definition: hcmap.h:42
t_hcmap::nmap
long int nmap
Definition: hcmap.h:39
t_hmi::hmidep
double hmidep
Definition: hmi.h:33
t_monitorresults
Definition: monitor_results.h:30
t_monitorresults::nSumErrorCaseMonitor
long int nSumErrorCaseMonitor
Definition: monitor_results.h:32
t_monitorresults::SumErrorCaseMonitor
double SumErrorCaseMonitor
Definition: monitor_results.h:31
t_optimize::nOptimiz
long int nOptimiz
Definition: optimize.h:246
t_optimize::lgOptimize
bool lgOptimize
Definition: optimize.h:253
t_pressure::RadBetaMax
realnum RadBetaMax
Definition: pressure.h:136
t_prt::lgPrintTime
bool lgPrintTime
Definition: prt.h:120
t_radius::depth_x_fillfac
double depth_x_fillfac
Definition: radius.h:74
t_radius::Radius
double Radius
Definition: radius.h:25
t_radius::depth
double depth
Definition: radius.h:38
t_rfield::EnergyDiffCont
realnum EnergyDiffCont
Definition: rfield.h:485
t_rfield::EnergyIncidCont
realnum EnergyIncidCont
Definition: rfield.h:484
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
t_struc::testr
realnum * testr
Definition: struc.h:25
t_struc::pressure
realnum * pressure
Definition: struc.h:33
t_thermal::ctot
double ctot
Definition: thermal.h:112
t_thermal::htot
double htot
Definition: thermal.h:149
thermal
t_thermal thermal
Definition: thermal.cpp:5
wind
Wind wind
Definition: wind.cpp:5
Generated by doxygen 1.9.5