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

Generated for cloudy by doxygen 1.8.4