00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "input.h"
00008 #include "conv.h"
00009 #include "optimize.h"
00010 #include "iso.h"
00011 #include "called.h"
00012 #include "atmdat.h"
00013 #include "hcmap.h"
00014 #include "thermal.h"
00015 #include "pressure.h"
00016 #include "struc.h"
00017 #include "wind.h"
00018 #include "h2.h"
00019 #include "colden.h"
00020 #include "dense.h"
00021 #include "lines.h"
00022 #include "secondaries.h"
00023 #include "radius.h"
00024 #include "version.h"
00025 #include "hmi.h"
00026 #include "prt.h"
00027 #include "grainvar.h"
00028 #include "atomfeii.h"
00029 #include "cddrive.h"
00030 #include "elementnames.h"
00031 #include "monitor_results.h"
00032 #include "taulines.h"
00033 #include "grid.h"
00034 #include "lines_service.h"
00035 #include "parser.h"
00036
00037
00038 static bool lgInitDone=false ,
00039
00040 lgSpaceAllocated=false;
00041
00042
00043 static const int NASSERTS = 200;
00044
00045
00046 static realnum ErrorDefault;
00047
00048
00049 static realnum ErrorDefaultPerformance;
00050
00051
00052 static const int NCHAR = 5;
00053 static char **chAssertLineLabel;
00054
00055
00056 static char *chAssertLimit;
00057
00058
00059 static char **chAssertType;
00060
00061
00062 static double *AssertQuantity ,*AssertQuantity2 ,*AssertError, **Param;
00063
00064
00065 static int *lgQuantityLog;
00066 static long nAsserts=0;
00067 static realnum *wavelength;
00068
00069 inline double ForcePass( char chAssertLimit1 )
00070 {
00071
00072 if( chAssertLimit1 == '=' )
00073 return 0.;
00074 else if( chAssertLimit1 == '<' )
00075 return 1.;
00076 else if( chAssertLimit1 == '>' )
00077 return -1.;
00078 else
00079 TotalInsanity();
00080 }
00081
00082
00083
00084 void InitMonitorResults(void)
00085 {
00086
00087
00088
00089 lgInitDone = true;
00090 nAsserts = 0;
00091
00092
00093
00094 ErrorDefault = 0.05f;
00095
00096 ErrorDefaultPerformance = 0.2f;
00097 }
00098
00099
00100
00101 void ParseMonitorResults(Parser &p)
00102 {
00103 long i ,
00104 nelem,
00105 n2;
00106 char chLabel[INPUT_LINE_LENGTH];
00107
00108 DEBUG_ENTRY( "ParseMonitorResults()" );
00109
00110 if( !lgInitDone )
00111 {
00112 fprintf( ioQQQ, " ParseMonitorResults called before InitAsserResults\n" );
00113 cdEXIT(EXIT_FAILURE);
00114 }
00115
00116
00117 if( !lgSpaceAllocated )
00118 {
00119
00120
00121 lgSpaceAllocated = true;
00122
00123
00124 chAssertLineLabel = ((char **)MALLOC(NASSERTS*sizeof(char *)));
00125
00126
00127 chAssertType = ((char **)MALLOC(NASSERTS*sizeof(char *)));
00128
00129
00130 Param = ((double **)MALLOC(NASSERTS*sizeof(double * )));
00131
00132
00133 for( i=0; i<NASSERTS; ++i )
00134 {
00135 chAssertLineLabel[i] = ((char *)MALLOC(NCHAR*sizeof(char )));
00136 strcpy(chAssertLineLabel[i],"unkn" );
00137
00138
00139 chAssertType[i] = ((char *)MALLOC(3*sizeof(char )));
00140 strcpy(chAssertType[i],"un" );
00141
00142
00143 Param[i] = ((double *)MALLOC(5*sizeof(double )));
00144 }
00145
00146
00147 AssertQuantity = ((double *)MALLOC(NASSERTS*sizeof(double )));
00148
00149 AssertQuantity2 = ((double *)MALLOC(NASSERTS*sizeof(double )));
00150
00151
00152 AssertError = ((double *)MALLOC(NASSERTS*sizeof(double )));
00153
00154
00155 wavelength = ((realnum *)MALLOC(NASSERTS*sizeof(realnum )));
00156
00157
00158 lgQuantityLog = ((int *)MALLOC(NASSERTS*sizeof(int )));
00159
00160
00161 chAssertLimit = ((char *)MALLOC(NASSERTS*sizeof(char )));
00162 }
00163
00164
00165
00166 vector<double> AssertVector;
00167 if( p.nMatch( "GRID" ) )
00168 {
00169 ASSERT( optimize.nOptimiz >= 0 );
00170 AssertVector.resize( optimize.nOptimiz+1 );
00171
00172
00173 p.GetQuote( chLabel, true );
00174
00175 bool lgEOF;
00176 input_readvector( chLabel, get_ptr(AssertVector), optimize.nOptimiz+1, &lgEOF );
00177 if( lgEOF )
00178 fprintf(ioQQQ,"PROBLEM the file %s does not have enough values. sorry\n", chLabel );
00179 }
00180
00181
00182 strcpy( chAssertType[nAsserts] , "un" );
00183
00184
00185
00186 lgQuantityLog[nAsserts] = false;
00187
00188
00189 if( p.nMatch("<" ) )
00190 {
00191 chAssertLimit[nAsserts] = '<';
00192 }
00193 else if( p.nMatch(">" ) )
00194 {
00195 chAssertLimit[nAsserts] = '>';
00196 }
00197 else
00198 {
00199 chAssertLimit[nAsserts] = '=';
00200 }
00201
00202
00203
00204 if( p.nMatch(" SET" ) )
00205 {
00206
00207
00208
00209
00210
00211 if( nAsserts >0 )
00212 fprintf(ioQQQ," The default monitor error is being changed after"
00213 " some asserts were entered. \n This will only affect asserts"
00214 " that come after this command.\n");
00215 --nAsserts;
00216
00217 if( p.nMatch("ERRO" ) )
00218 {
00219 if( p.nMatch("PERF" ) )
00220 {
00221
00222 ErrorDefaultPerformance =
00223 (realnum)p.FFmtRead();
00224 if( p.lgEOL() )
00225 p.NoNumb("error");
00226 }
00227
00228 ErrorDefault =
00229 (realnum)p.FFmtRead();
00230 if( p.lgEOL() )
00231 p.NoNumb("error");
00232 }
00233 else
00234 {
00235
00236 fprintf( ioQQQ,
00237 " I could not identify an option on this ASSERT SET XXX command.\n");
00238 fprintf( ioQQQ, " Sorry.\n" );
00239 cdEXIT(EXIT_FAILURE);
00240 }
00241 }
00242
00243
00244 else if( p.nMatch("IONI" ) )
00245 {
00246
00247
00248
00249
00250
00251
00252 if( p.nMatch("VOLU" ) )
00253 {
00254 strcpy( chAssertType[nAsserts] , "F " );
00255 }
00256 else
00257 {
00258
00259 strcpy( chAssertType[nAsserts] , "f " );
00260 }
00261
00262
00263 if( (nelem = p.GetElem()) < 0 )
00264 {
00265 fprintf( ioQQQ,
00266 " I could not identify an element name on this line.\n");
00267 fprintf( ioQQQ, " Sorry.\n" );
00268 cdEXIT(EXIT_FAILURE);
00269 }
00270 ASSERT( nelem>= 0);
00271 ASSERT( nAsserts>= 0);
00272
00273
00274 strcpy( chAssertLineLabel[nAsserts], elementnames.chElementNameShort[nelem] );
00275
00276
00277 wavelength[nAsserts] = (realnum)p.FFmtRead();
00278 if( p.lgEOL() )
00279 {
00280 p.NoNumb("ionization stage");
00281 }
00282
00283 if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 )
00284 {
00285 fprintf( ioQQQ,
00286 " The ionization stage is inappropriate for this element.\n");
00287 fprintf( ioQQQ, " Sorry.\n" );
00288 cdEXIT(EXIT_FAILURE);
00289 }
00290
00291 if( p.nMatch( "GRID" ) )
00292 {
00293 AssertQuantity[nAsserts] = AssertVector[optimize.nOptimiz];
00294 }
00295 else
00296 {
00297
00298
00299 AssertQuantity[nAsserts] = p.FFmtRead();
00300 if( p.lgEOL() )
00301 p.NoNumb("ionization fraction");
00302 }
00303
00304
00305 AssertError[nAsserts] = p.FFmtRead();
00306 if( p.lgEOL() )
00307 AssertError[nAsserts] = ErrorDefault;
00308
00309
00310
00311 if( AssertQuantity[nAsserts] <= 0. )
00312 {
00313
00314 AssertQuantity[nAsserts] =
00315 pow(10.,AssertQuantity[nAsserts] );
00316
00317 lgQuantityLog[nAsserts] = true;
00318 }
00319 else if( AssertQuantity[nAsserts] > 1. )
00320 {
00321 fprintf( ioQQQ,
00322 " The ionization fraction must be less than one.\n");
00323 fprintf( ioQQQ, " Sorry.\n" );
00324 cdEXIT(EXIT_FAILURE);
00325 }
00326
00327
00328 if( fabs(AssertQuantity[nAsserts]) <= SMALLDOUBLE )
00329 {
00330 fprintf( ioQQQ,
00331 " The ionization ionization fraction is too small, or zero. Check input\n");
00332 fprintf( ioQQQ, " Sorry.\n" );
00333 cdEXIT(EXIT_FAILURE);
00334 }
00335 }
00336
00337
00338 else if( p.nMatch("MOLE" )&&p.nMatch("FRAC" ) )
00339 {
00340
00341
00342
00343
00344
00345
00347 if( p.nMatch("VOLU" ) )
00348 {
00349 strcpy( chAssertType[nAsserts] , "MF" );
00350 }
00351 else
00352 {
00353
00354 strcpy( chAssertType[nAsserts] , "mf" );
00355 }
00356
00357 if( p.nMatchErase(" H2 " ) )
00358 {
00359 strcpy( chAssertLineLabel[nAsserts], "H2 " );
00360
00361 }
00362 else if( p.nMatch(" CO " ) )
00363 {
00364 strcpy( chAssertLineLabel[nAsserts], "CO " );
00365 }
00366 else
00367 {
00368 fprintf( ioQQQ,
00369 " I could not identify CO or H2 on this line.\n");
00370 fprintf( ioQQQ, " Sorry.\n" );
00371 cdEXIT(EXIT_FAILURE);
00372 }
00373
00374
00375 wavelength[nAsserts] = 0;
00376
00377
00378 AssertQuantity[nAsserts] = p.FFmtRead();
00379 if( p.lgEOL() )
00380 {
00381 p.NoNumb("molecular fraction");
00382 }
00383 if( AssertQuantity[nAsserts] <= 0. )
00384 {
00385
00386 AssertQuantity[nAsserts] =
00387 pow(10.,AssertQuantity[nAsserts] );
00388 }
00389
00390
00391
00392 AssertError[nAsserts] = p.FFmtRead();
00393 if( p.lgEOL() )
00394 {
00395
00396 AssertError[nAsserts] = ErrorDefault;
00397 }
00398
00399 lgQuantityLog[nAsserts] = true;
00400 }
00401
00402
00403 else if( p.nMatch(" LINE " ) )
00404 {
00405 if( p.nMatch("LUMI") || p.nMatch("INTE"))
00406 {
00407
00408 strcpy( chAssertType[nAsserts] , "Ll" );
00409
00410
00411 lgQuantityLog[nAsserts] = true;
00412 }
00413 else
00414 {
00415
00416 strcpy( chAssertType[nAsserts] , "Lr" );
00417
00418
00419 lgQuantityLog[nAsserts] = false;
00420 }
00421
00422
00423
00424 p.GetQuote( chLabel , true );
00425
00426
00427 if( chLabel[4] != '\0' )
00428 {
00429 fprintf( ioQQQ,
00430 " The label must be exactly 4 char long, between double quotes.\n");
00431 fprintf( ioQQQ, " Sorry.\n" );
00432 cdEXIT(EXIT_FAILURE);
00433 }
00434
00435
00436 strcpy( chAssertLineLabel[nAsserts], chLabel );
00437
00438
00439 wavelength[nAsserts] = (realnum)p.getWave();
00440
00441
00442
00443 AssertQuantity[nAsserts] = p.FFmtRead();
00444 if( p.lgEOL() )
00445 {
00446 p.NoNumb("intensity/luminosity");
00447 }
00448
00449 if( lgQuantityLog[nAsserts] )
00450 {
00451 if( AssertQuantity[nAsserts] > DBL_MAX_10_EXP ||
00452 AssertQuantity[nAsserts] < -DBL_MAX_10_EXP )
00453 {
00454 fprintf( ioQQQ,
00455 " The asserted quantity is a log, but is too large or small, value is %e.\n",
00456 AssertQuantity[nAsserts] );
00457 fprintf( ioQQQ, " I would crash if I tried to use it.\n Sorry.\n" );
00458 cdEXIT(EXIT_FAILURE);
00459 }
00460 AssertQuantity[nAsserts] =
00461 pow(10.,AssertQuantity[nAsserts] );
00462 }
00463 if( AssertQuantity[nAsserts]<= 0. )
00464 {
00465 fprintf( ioQQQ,
00466 " The relative intensity must be positive, and was %e.\n",AssertQuantity[nAsserts] );
00467 fprintf( ioQQQ, " Sorry.\n" );
00468 cdEXIT(EXIT_FAILURE);
00469 }
00470
00471
00472 AssertError[nAsserts] = p.FFmtRead();
00473 if( p.lgEOL() )
00474 {
00475
00476 AssertError[nAsserts] = ErrorDefault;
00477 }
00478 }
00479
00480
00481 else if( p.nMatch("CASE" ) )
00482 {
00483
00484 strcpy( chAssertType[nAsserts] , "CB" );
00485
00486 AssertError[nAsserts] = p.FFmtRead();
00487 if( p.lgEOL() )
00488
00489 AssertError[nAsserts] = ErrorDefault;
00490 AssertQuantity[nAsserts] = 0;
00491 wavelength[nAsserts] = 0.;
00492
00493
00494
00495 if( p.GetParam("FAINT",&Param[nAsserts][4]) )
00496 {
00497 if( p.lgEOL() )
00498 {
00499
00500 fprintf(ioQQQ," The monitor Case B faint option must have a number,"
00501 " the relative intensity of the fainest line to monitor.\n");
00502 cdEXIT(EXIT_FAILURE);
00503 }
00504
00505 if( Param[nAsserts][4]<=0. )
00506 Param[nAsserts][4] = pow(10., Param[nAsserts][4] );
00507 }
00508 else
00509 {
00510
00511 Param[nAsserts][4] = SMALLFLOAT;
00512 }
00513
00514
00515 if( p.GetRange("RANG",&Param[nAsserts][2],&Param[nAsserts][3]) )
00516 {
00517 if( p.lgEOL() )
00518 {
00519
00520 fprintf(ioQQQ," The monitor Case B range option must have two numbers,"
00521 " the lower and upper limit to the wavelengths in Angstroms.\n");
00522 fprintf(ioQQQ," There must be a total of three numbers on the line,"
00523 " the relative error followed by the lower and upper limits to the "
00524 "wavelength in Angstroms.\n");
00525 cdEXIT(EXIT_FAILURE);
00526 }
00527 if( Param[nAsserts][2]>Param[nAsserts][3])
00528 {
00529
00530 double sav = Param[nAsserts][3];
00531 Param[nAsserts][3] = Param[nAsserts][2];
00532 Param[nAsserts][2] = sav;
00533 }
00534 }
00535 else
00536 {
00537
00538 Param[nAsserts][2] = 0.;
00539 Param[nAsserts][3] = 1e30;
00540 }
00541
00542 if( p.nMatch("H-LI" ) )
00543 {
00544
00545 if( (nelem = p.GetElem()) < 0 )
00546 {
00547
00548 fprintf(ioQQQ, "monitor H-like case B did not find an element on this line, sorry\n");
00549 p.PrintLine(ioQQQ);
00550 cdEXIT(EXIT_FAILURE);
00551 }
00552 if( nelem>7 )
00553 {
00554
00555 fprintf(ioQQQ, "monitor H-like cannot do elements heavier than O, sorry\n");
00556 p.PrintLine(ioQQQ);
00557 cdEXIT(EXIT_FAILURE);
00558 }
00559 Param[nAsserts][0] = ipH_LIKE;
00560 Param[nAsserts][1] = nelem;
00561
00562 strcpy( chAssertLineLabel[nAsserts], "Ca ");
00563 if( p.nMatch("CASE A " ) )
00564 strcat( chAssertLineLabel[nAsserts] , "A");
00565 else
00566 strcat( chAssertLineLabel[nAsserts] , "B");
00567 }
00568 else if( p.nMatch("HE-L") )
00569 {
00570
00571 Param[nAsserts][0] = ipHE_LIKE;
00572 Param[nAsserts][1] = ipHELIUM;
00573
00574 strcpy( chAssertLineLabel[nAsserts] , "+Col");
00575 }
00576 else
00577 {
00578
00579 fprintf( ioQQQ,
00580 " I could not identify an iso-sequence on this Case A/B command.\n");
00581 fprintf( ioQQQ, " Sorry.\n" );
00582 cdEXIT(EXIT_FAILURE);
00583 }
00584 }
00585 else if( p.nMatch("DEPA") )
00586 {
00587
00588 AssertQuantity[nAsserts] = p.FFmtRead();
00589 if( p.lgEOL() )
00590 p.NoNumb("average departure coefficient");
00591
00592
00593 AssertError[nAsserts] = p.FFmtRead();
00594 if( p.lgEOL() )
00595
00596 AssertError[nAsserts] = ErrorDefault;
00597
00598
00599 if( p.nMatch("H-LI" ) )
00600 {
00601
00602 strcpy( chAssertLineLabel[nAsserts] , "dHyd" );
00603
00604 strcpy( chAssertType[nAsserts] , "D " );
00605
00606 if( (wavelength[nAsserts] = (realnum)p.GetElem()) < 0 )
00607 {
00608 fprintf( ioQQQ,
00609 " I could not identify an element name on this line.\n");
00610 fprintf( ioQQQ, " Sorry.\n" );
00611 cdEXIT(EXIT_FAILURE);
00612 }
00613 }
00614
00615
00616 else if( p.nMatch("HE-L" ) )
00617 {
00618
00619 strcpy( chAssertLineLabel[nAsserts] , "d He" );
00620
00621 strcpy( chAssertType[nAsserts] , "De" );
00622
00623 if( (wavelength[nAsserts] = (realnum)p.GetElem()) < 0 )
00624 {
00625 fprintf( ioQQQ,
00626 " I could not identify an element name on this line.\n");
00627 fprintf( ioQQQ, " Sorry.\n" );
00628 cdEXIT(EXIT_FAILURE);
00629 }
00630 }
00631
00632
00633 else if( p.nMatch("FEII") || p.nMatch("FE II" ) )
00634 {
00635
00636 strcpy( chAssertLineLabel[nAsserts] , "d Fe" );
00637
00638 strcpy( chAssertType[nAsserts] , "d " );
00639
00640 wavelength[nAsserts] = 2;
00641 }
00642
00643
00644 else if( p.nMatch("HMIN" ) )
00645 {
00646
00647 strcpy( chAssertLineLabel[nAsserts] , "d H-" );
00648
00649 strcpy( chAssertType[nAsserts] , "d-" );
00650
00651 wavelength[nAsserts] = -1;
00652 }
00653 else
00654 {
00655 fprintf( ioQQQ,
00656 " There must be a second key: FEII, H-LIke, HMINus, or HE-Like followed by element.\n");
00657 fprintf( ioQQQ, " Sorry.\n" );
00658 cdEXIT(EXIT_FAILURE);
00659 }
00660
00661
00662 if( p.nMatch("EXCI" ) )
00663 {
00664
00665 AssertQuantity2[nAsserts] = 1.;
00666 }
00667 else
00668 {
00669
00670 AssertQuantity2[nAsserts] = 0.;
00671 }
00672 }
00673
00674
00675 else if( p.nMatch(" MAP" ) )
00676 {
00677
00678
00679
00680 if( p.nMatch("HEAT" ) )
00681 {
00682 strcpy( chAssertType[nAsserts] , "mh" );
00683 }
00684 else if( p.nMatch("COOL" ) )
00685 {
00686 strcpy( chAssertType[nAsserts] , "mc" );
00687 }
00688 else
00689 {
00690 fprintf( ioQQQ,
00691 " There must be a second key, HEATing or COOLing.\n");
00692 fprintf( ioQQQ, " Sorry.\n" );
00693 cdEXIT(EXIT_FAILURE);
00694 }
00695
00696
00697 AssertQuantity2[nAsserts] = p.FFmtRead();
00698 if( p.lgEOL() )
00699 {
00700 p.NoNumb("temperature");
00701 }
00702
00703 if( AssertQuantity2[nAsserts] <= 10. )
00704 {
00705
00706 AssertQuantity2[nAsserts] =
00707 pow(10.,AssertQuantity2[nAsserts] );
00708 }
00709
00710
00711 wavelength[nAsserts] = (realnum)AssertQuantity2[nAsserts];
00712
00713
00714 AssertQuantity[nAsserts] = p.FFmtRead();
00715 if( p.lgEOL() )
00716 {
00717 p.NoNumb("heating/cooling");
00718 }
00719
00720
00721 AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts]);
00722
00723
00724
00725 AssertError[nAsserts] = p.FFmtRead();
00726 if( p.lgEOL() )
00727 {
00728
00729 AssertError[nAsserts] = ErrorDefault;
00730 }
00731
00732
00733 lgQuantityLog[nAsserts] = true;
00734 }
00735
00736
00737 else if( p.nMatch("COLU" ) )
00738 {
00739
00740 strcpy( chAssertType[nAsserts] , "cd" );
00741
00742
00743 wavelength[nAsserts] = 0;
00744
00745 char chLabel[INPUT_LINE_LENGTH];
00746
00747 if ( p.GetQuote(chLabel,false) == 0 )
00748 {
00749
00750 for (long i=strlen(chLabel);i<NCHAR-1;++i)
00751 chLabel[i] = ' ';
00752 chLabel[NCHAR-1] = '\0';
00753 strcpy( chAssertLineLabel[nAsserts], chLabel );
00754 }
00755 else if( p.nMatchErase(" H2 ") )
00756 {
00757 strcpy( chAssertLineLabel[nAsserts], "H2 " );
00758 }
00759 else if( p.nMatchErase("H3+ "))
00760 {
00761 strcpy( chAssertLineLabel[nAsserts], "H3+ " );
00762 }
00763 else if( p.nMatchErase("H2+ "))
00764 {
00765 strcpy( chAssertLineLabel[nAsserts], "H2+ " );
00766 }
00767 else if( p.nMatchErase(" H- "))
00768 {
00769 strcpy( chAssertLineLabel[nAsserts], "H- " );
00770 }
00771 else if( p.nMatchErase("H2G "))
00772 {
00773 strcpy( chAssertLineLabel[nAsserts], "H2g " );
00774 }
00775 else if( p.nMatchErase("H2* "))
00776 {
00777 strcpy( chAssertLineLabel[nAsserts], "H2* " );
00778 }
00779 else if( p.nMatchErase("HEH+"))
00780 {
00781 strcpy( chAssertLineLabel[nAsserts], "HeH+" );
00782 }
00783 else if( p.nMatchErase(" O2 "))
00784 {
00785 strcpy( chAssertLineLabel[nAsserts], "O2 " );
00786 }
00787 else if( p.nMatchErase("H2O "))
00788 {
00789 strcpy( chAssertLineLabel[nAsserts], "H2O " );
00790 }
00791 else if( p.nMatchErase(" C2 "))
00792 {
00793 strcpy( chAssertLineLabel[nAsserts], "C2 " );
00794 }
00795 else if( p.nMatchErase(" C3 "))
00796 {
00797 strcpy( chAssertLineLabel[nAsserts], "C3 " );
00798 }
00799 else if( p.nMatch(" CO "))
00800 {
00801 strcpy( chAssertLineLabel[nAsserts], "CO " );
00802 }
00803 else if( p.nMatch("SIO "))
00804 {
00805 strcpy( chAssertLineLabel[nAsserts], "SiO " );
00806 }
00807 else if( p.nMatch(" OH "))
00808 {
00809 strcpy( chAssertLineLabel[nAsserts], "OH " );
00810 }
00811 else if( p.nMatch(" CN ") )
00812 {
00813 strcpy( chAssertLineLabel[nAsserts], "CN " );
00814 }
00815 else if( p.nMatch(" CH ") )
00816 {
00817 strcpy( chAssertLineLabel[nAsserts], "CH " );
00818 }
00819 else if( p.nMatch(" CH+") )
00820 {
00821 strcpy( chAssertLineLabel[nAsserts], "CH+ " );
00822 }
00823 else
00824 {
00825 fprintf( ioQQQ,
00826 " I could not identify H2, H3+, H2+, H2g, H2*, H2H+, CO, O2, SiO, OH, C2 or C3 or on this line.\n");
00827 fprintf( ioQQQ, " Sorry.\n" );
00828 cdEXIT(EXIT_FAILURE);
00829 }
00830
00831 if( p.nMatch( "LEVE" ) )
00832 {
00833 if (strncmp( chAssertLineLabel[nAsserts], "H2 ", NCHAR) != 0)
00834 {
00835 fprintf( ioQQQ, " LEVEL option only implemented for H2.\n");
00836 fprintf( ioQQQ, " Sorry.\n" );
00837 cdEXIT(EXIT_FAILURE);
00838 }
00839
00840
00841 Param[nAsserts][0] = p.FFmtRead();
00842 if( p.lgEOL() )
00843 p.NoNumb("level v" );
00844 Param[nAsserts][1] = p.FFmtRead();
00845 if( p.lgEOL() )
00846 p.NoNumb("level J" );
00847
00848 wavelength[nAsserts] = (realnum)(100.*Param[nAsserts][0] + Param[nAsserts][1]);
00849 }
00850 else
00851 {
00852
00853 Param[nAsserts][0] = -1.;
00854 Param[nAsserts][1] = -1.;
00855 }
00856
00857
00858
00859 AssertQuantity[nAsserts] = p.FFmtRead();
00860 if( p.lgEOL() )
00861 {
00862 p.NoNumb("column density");
00863 }
00864
00865 AssertQuantity[nAsserts] =
00866 pow(10.,AssertQuantity[nAsserts] );
00867
00868
00869
00870 AssertError[nAsserts] = p.FFmtRead();
00871 if( p.lgEOL() )
00872 {
00873
00874 AssertError[nAsserts] = ErrorDefault;
00875 }
00876
00877
00878
00879 if( p.nMatch( " LOG" ) )
00880 {
00881 AssertError[nAsserts] = pow(10. , AssertError[nAsserts] );
00882 }
00883
00884
00885 lgQuantityLog[nAsserts] = true;
00886 }
00887
00888
00889 else if( p.nMatch("GRAI") && p.nMatch(" H2 ") )
00890 {
00891
00892 strcpy( chAssertType[nAsserts] , "g2" );
00893
00894 strcpy( chAssertLineLabel[nAsserts], "R H2" );
00895
00896 nelem = (long int)p.FFmtRead();
00897 if( nelem!=2 )
00898 {
00899 fprintf( ioQQQ,
00900 " I did not find a 2, as in H2, as the first number on this line.\n");
00901 fprintf( ioQQQ,
00902 " The rate should be the second number.\n");
00903 fprintf( ioQQQ, " Sorry.\n" );
00904 cdEXIT(EXIT_FAILURE);
00905 }
00906
00907 AssertQuantity[nAsserts] = p.FFmtRead();
00908 if( p.lgEOL() )
00909 {
00910 p.NoNumb("grain H2 formation");
00911 }
00912
00913 if( AssertQuantity[nAsserts] <=0. )
00914 AssertQuantity[nAsserts] = pow(10.,AssertQuantity[nAsserts] );
00915
00916 wavelength[nAsserts] = 0;
00917
00918
00919
00920 AssertError[nAsserts] = p.FFmtRead();
00921 if( p.lgEOL() )
00922 {
00923
00924 AssertError[nAsserts] = ErrorDefault;
00925
00926 lgQuantityLog[nAsserts] = true;
00927 }
00928 }
00929
00930
00931 else if( p.nMatch( "GRAI" ) && p.nMatch( "POTE") )
00932 {
00933
00934 strcpy( chAssertType[nAsserts] , "gp" );
00935
00936 strcpy( chAssertLineLabel[nAsserts], "GPot" );
00937
00938
00939 wavelength[nAsserts] = (realnum)p.FFmtRead();
00940
00941 AssertQuantity[nAsserts] = p.FFmtRead();
00942
00943 if( p.lgEOL() )
00944 {
00945 p.NoNumb("grain potential");
00946 }
00947
00948
00949
00950 AssertError[nAsserts] = p.FFmtRead();
00951
00952 if( p.lgEOL() )
00953 {
00954
00955 AssertError[nAsserts] = ErrorDefault;
00956 }
00957 }
00958
00959
00960 else if( p.nMatch("TEMP") )
00961 {
00962
00963
00964
00965
00966
00967 if( p.nMatch("VOLU") )
00968 {
00969 strcpy( chAssertType[nAsserts] , "T ");
00970 }
00971 else
00972 {
00973
00974 strcpy( chAssertType[nAsserts] , "t ");
00975 }
00976
00977
00978
00979 if( p.nMatch("GRAI") )
00980 {
00981 strcpy( chAssertLineLabel[nAsserts], "GTem" );
00982
00983
00984 nelem = LONG_MAX-2;
00985
00986
00987
00988
00989
00990 }
00991
00992
00993 else if( p.nMatch("FACE") )
00994 {
00995 strcpy( chAssertLineLabel[nAsserts], "face" );
00996
00997 nelem = 0;
00998 }
00999
01000
01001 else if( p.nMatch(" H2 ") )
01002 {
01003 strcpy( chAssertLineLabel[nAsserts], "H2 " );
01004
01005 nelem = 0;
01006 }
01007
01008
01009
01010 else if( (nelem = p.GetElem()) >= 0 )
01011 {
01012
01013
01014 strcpy( chAssertLineLabel[nAsserts], elementnames.chElementNameShort[nelem] );
01015 }
01016 else
01017 {
01018 fprintf( ioQQQ,
01019 " I could not identify an element name on this line.\n");
01020 fprintf( ioQQQ, " Sorry.\n" );
01021 cdEXIT(EXIT_FAILURE);
01022 }
01023
01024
01025 if( strcmp( chAssertLineLabel[nAsserts], "face" )== 0)
01026 {
01027
01028 wavelength[nAsserts] = 0;
01029 }
01030 else if( strcmp( chAssertLineLabel[nAsserts], "H2 " )== 0)
01031 {
01032 int j;
01033
01034 wavelength[nAsserts] = 0;
01035
01036 j = (int)p.FFmtRead();
01037 if( j != 2 )
01038 TotalInsanity();
01039 ++i;
01040 }
01041 else
01042 {
01043 wavelength[nAsserts] = (realnum)p.FFmtRead();
01044 if( p.lgEOL() )
01045 {
01046 p.NoNumb("ionization stage");
01047 }
01048
01049 if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 )
01050 {
01051 fprintf( ioQQQ,
01052 " The ionization stage is inappropriate for this element.\n");
01053 fprintf( ioQQQ, " Sorry.\n" );
01054 cdEXIT(EXIT_FAILURE);
01055 }
01056 }
01057
01058 if( p.nMatch( "GRID" ) )
01059 {
01060 AssertQuantity[nAsserts] = AssertVector[optimize.nOptimiz];
01061 }
01062 else
01063 {
01064
01065 AssertQuantity[nAsserts] = p.FFmtRead();
01066 if( p.lgEOL() )
01067 p.NoNumb("temperature");
01068 }
01069
01070
01071 AssertError[nAsserts] = p.FFmtRead();
01072 if( p.lgEOL() )
01073 AssertError[nAsserts] = ErrorDefault;
01074
01075
01076
01077 if( AssertQuantity[nAsserts] <= 10. && !p.nMatch( "LINE" ) )
01078 {
01079
01080 AssertQuantity[nAsserts] =
01081 pow(10.,AssertQuantity[nAsserts] );
01082
01083 lgQuantityLog[nAsserts] = true;
01084 }
01085 }
01086
01087
01088 else if( p.nMatch("HHEI") )
01089 {
01090
01091 strcpy( chAssertType[nAsserts] , "hh" );
01092
01093 strcpy( chAssertLineLabel[nAsserts], "HHei" );
01094
01095
01096
01097 AssertQuantity[nAsserts] = p.FFmtRead();
01098 if( p.lgEOL() )
01099 {
01100 p.NoNumb("ionization correction factor");
01101 }
01102
01103 wavelength[nAsserts] = 0;
01104
01105
01106
01107 AssertError[nAsserts] = p.FFmtRead();
01108 if( p.lgEOL() )
01109 {
01110
01111 AssertError[nAsserts] = ErrorDefault;
01112 }
01113 }
01114
01115
01116 else if( p.nMatch(" H2 ") )
01117 {
01118
01119 if( p.nMatch("ORTH") )
01120 {
01121
01122 strcpy( chAssertType[nAsserts] , "or" );
01123
01124 strcpy( chAssertLineLabel[nAsserts], "orth" );
01125 }
01126 else
01127 {
01128 fprintf( ioQQQ,
01129 " I could not identify a second keyword on this line.\n");
01130 fprintf( ioQQQ, " Sorry.\n" );
01131 cdEXIT(EXIT_FAILURE);
01132 }
01133
01134
01135 n2 = (long)p.FFmtRead();
01136 if( p.lgEOL() || n2 != 2 )
01137 {
01138 p.NoNumb("the 2 in H2 ?!");
01139 }
01140 AssertQuantity[nAsserts] = p.FFmtRead();
01141
01142 wavelength[nAsserts] = 0;
01143
01144
01145
01146 AssertError[nAsserts] = p.FFmtRead();
01147 if( p.lgEOL() )
01148 {
01149
01150 AssertError[nAsserts] = ErrorDefault;
01151 }
01152 }
01153
01154
01155 else if( p.nMatch(" MPI") )
01156 {
01157
01158 strcpy( chAssertType[nAsserts] , "mp" );
01159
01160 strcpy( chAssertLineLabel[nAsserts], "mpi " );
01161
01162 wavelength[nAsserts] = 0.;
01163 AssertQuantity[nAsserts] = 1.;
01164 AssertError[nAsserts] = ErrorDefault;
01165 }
01166
01167
01168 else if( p.nMatch("NZON") )
01169 {
01170
01171 strcpy( chAssertType[nAsserts] , "z " );
01172
01173 strcpy( chAssertLineLabel[nAsserts], "zone" );
01174
01175
01176 wavelength[nAsserts] = 0.;
01177 AssertQuantity[nAsserts] = p.FFmtRead();
01178 if( p.lgEOL() )
01179 p.NoNumb("zone number");
01180
01181
01182 AssertError[nAsserts] = p.FFmtRead();
01183 if( p.lgEOL() )
01184 AssertError[nAsserts] = ErrorDefault;
01185 }
01186
01187
01188 else if( p.nMatch("PRES") && p.nMatch("ERRO") )
01189 {
01190
01191 strcpy( chAssertType[nAsserts] , "pr" );
01192
01193 strcpy( chAssertLineLabel[nAsserts], "pres" );
01194
01195
01196
01197 wavelength[nAsserts] = 0;
01198 AssertQuantity[nAsserts] = (double)p.FFmtRead();
01199 if( p.lgEOL() )
01200 {
01201 p.NoNumb("pressure error");
01202 }
01203 else if( AssertQuantity[nAsserts] <= 0.)
01204 {
01205
01206 AssertQuantity[nAsserts] = pow(10. , AssertQuantity[nAsserts]);
01207 }
01208
01209
01210
01211 AssertError[nAsserts] = p.FFmtRead();
01212 if( p.lgEOL() )
01213 {
01214
01215 AssertError[nAsserts] = ErrorDefault;
01216 }
01217 }
01218
01219 else if( p.nMatch("PRADMAX") )
01220 {
01221
01222
01223 strcpy( chAssertType[nAsserts] , "RM" );
01224
01225 strcpy( chAssertLineLabel[nAsserts], "Prad" );
01226
01227
01228
01229 wavelength[nAsserts] = 0;
01230 AssertQuantity[nAsserts] = (double)p.FFmtRead();
01231 if( p.lgEOL() )
01232 {
01233 p.NoNumb("PRADMAX");
01234 }
01235 else if( AssertQuantity[nAsserts] <= 0.)
01236 {
01237
01238 AssertQuantity[nAsserts] = pow(10. , AssertQuantity[nAsserts]);
01239 }
01240
01241
01242
01243 AssertError[nAsserts] = p.FFmtRead();
01244 if( p.lgEOL() )
01245 {
01246
01247 AssertError[nAsserts] = ErrorDefault;
01248 }
01249 }
01250
01251
01252 else if( p.nMatch("CSUP") )
01253 {
01254
01255 strcpy( chAssertType[nAsserts] , "sc" );
01256
01257 strcpy( chAssertLineLabel[nAsserts], "sion" );
01258
01259
01260 AssertQuantity[nAsserts] = p.FFmtRead();
01261 if( p.lgEOL() )
01262 {
01263 p.NoNumb("secondary ionization rate");
01264 }
01265
01266 AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts] );
01267
01268
01269 wavelength[nAsserts] = 0;
01270
01271
01272
01273 AssertError[nAsserts] = p.FFmtRead();
01274 if( p.lgEOL() )
01275 {
01276
01277 AssertError[nAsserts] = ErrorDefault;
01278 }
01279
01280
01281 lgQuantityLog[nAsserts] = true;
01282 }
01283
01284
01285 else if( p.nMatch("HTOT") )
01286 {
01287
01288 strcpy( chAssertType[nAsserts] , "ht" );
01289
01290
01291 strcpy( chAssertLineLabel[nAsserts], "heat" );
01292
01293
01294 AssertQuantity[nAsserts] = p.FFmtRead();
01295 if( p.lgEOL() )
01296 {
01297 p.NoNumb("heating rate");
01298 }
01299
01300 AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts] );
01301
01302
01303 wavelength[nAsserts] = 0;
01304
01305
01306
01307 AssertError[nAsserts] = p.FFmtRead();
01308 if( p.lgEOL() )
01309 {
01310
01311 AssertError[nAsserts] = ErrorDefault;
01312 }
01313
01314
01315 lgQuantityLog[nAsserts] = true;
01316 }
01317
01318
01319 else if( p.nMatch("ITRZ") )
01320 {
01321
01322 strcpy( chAssertType[nAsserts] , "iz" );
01323
01324 strcpy( chAssertLineLabel[nAsserts], "itrz" );
01325
01326
01327 AssertQuantity[nAsserts] = p.FFmtRead();
01328 if( p.lgEOL() )
01329 {
01330 p.NoNumb("iterations per zone");
01331 }
01332
01333 wavelength[nAsserts] = 0;
01334
01335
01336 AssertError[nAsserts] = p.FFmtRead();
01337 if( p.lgEOL() )
01338 AssertError[nAsserts] = ErrorDefaultPerformance;
01339 }
01340
01341
01342 else if( p.nMatch("EDEN") )
01343 {
01344
01345 strcpy( chAssertType[nAsserts] , "e " );
01346
01347 strcpy( chAssertLineLabel[nAsserts], "eden" );
01348
01349
01350 AssertQuantity[nAsserts] =
01351 pow(10., p.FFmtRead() );
01352 if( p.lgEOL() )
01353 {
01354 p.NoNumb(" electron density of the last zone");
01355 }
01356
01357
01358
01359 AssertError[nAsserts] = p.FFmtRead();
01360 if( p.lgEOL() )
01361 {
01362
01363 AssertError[nAsserts] = ErrorDefault;
01364 }
01365 wavelength[nAsserts] = 0;
01366
01367
01368 lgQuantityLog[nAsserts] = true;
01369 }
01370
01371
01372 else if( p.nMatch(" TU ") )
01373 {
01374
01375 strcpy( chAssertType[nAsserts] , "Tu" );
01376 strcpy( chAssertLineLabel[nAsserts], "Tu " );
01377
01378
01379 AssertQuantity[nAsserts] = p.FFmtRead();
01380 if( p.lgEOL() )
01381 p.NoNumb("energy density of last zone");
01382
01383
01384
01385 lgQuantityLog[nAsserts] = false;
01386 if( AssertQuantity[nAsserts] <= 10. && !p.nMatch( "LINE" ) )
01387 {
01388
01389 AssertQuantity[nAsserts] =
01390 pow(10.,AssertQuantity[nAsserts] );
01391
01392 }
01393
01394
01395
01396 AssertError[nAsserts] = p.FFmtRead();
01397 if( p.lgEOL() )
01398
01399 AssertError[nAsserts] = ErrorDefault;
01400 wavelength[nAsserts] = 0;
01401 }
01402
01403
01404 else if( p.nMatch("THIC") || p.nMatch("DEPT") )
01405 {
01406
01407 strcpy( chAssertType[nAsserts] , "th" );
01408
01409 strcpy( chAssertLineLabel[nAsserts], "thic" );
01410
01411
01412 AssertQuantity[nAsserts] = pow(10., p.FFmtRead() );
01413 if( p.lgEOL() )
01414 {
01415 p.NoNumb("thickness or depth of model");
01416 }
01417
01418
01419
01420 AssertError[nAsserts] = p.FFmtRead();
01421 if( p.lgEOL() )
01422 {
01423
01424 AssertError[nAsserts] = ErrorDefault;
01425 }
01426 wavelength[nAsserts] = 0;
01427
01428
01429 lgQuantityLog[nAsserts] = true;
01430 }
01431
01432
01433 else if( p.nMatch("RADI") )
01434 {
01435
01436 strcpy( chAssertType[nAsserts] , "ra" );
01437
01438 strcpy( chAssertLineLabel[nAsserts], "radi" );
01439
01440
01441 AssertQuantity[nAsserts] = pow(10., p.FFmtRead() );
01442 if( p.lgEOL() )
01443 {
01444 p.NoNumb("outer radius");
01445 }
01446
01447
01448
01449 AssertError[nAsserts] = p.FFmtRead();
01450 if( p.lgEOL() )
01451 {
01452
01453 AssertError[nAsserts] = ErrorDefault;
01454 }
01455 wavelength[nAsserts] = 0;
01456
01457
01458 lgQuantityLog[nAsserts] = true;
01459 }
01460
01461
01462 else if( p.nMatch("NITE") )
01463 {
01464
01465 strcpy( chAssertType[nAsserts] , "Z " );
01466
01467 strcpy( chAssertLineLabel[nAsserts], "iter" );
01468
01469 wavelength[nAsserts] = 0.f;
01470
01471
01472 AssertQuantity[nAsserts] = p.FFmtRead();
01473 if( p.lgEOL() )
01474 {
01475 p.NoNumb("number of iterations");
01476 }
01477
01478
01479
01480 AssertError[nAsserts] = p.FFmtRead();
01481 if( p.lgEOL() )
01482 {
01483
01484 AssertError[nAsserts] = ErrorDefault;
01485 }
01486 }
01487
01488
01489
01490 else if( p.nMatch("VELO") )
01491 {
01492
01493 strcpy( chAssertType[nAsserts] , "Ve" );
01494
01495 strcpy( chAssertLineLabel[nAsserts], "vel " );
01496 wavelength[nAsserts] = 0;
01497
01498
01499 AssertQuantity[nAsserts] = p.FFmtRead();
01500 if( p.lgEOL() )
01501 p.NoNumb("terminal velocity");
01502
01503
01504
01505 AssertError[nAsserts] = p.FFmtRead();
01506 if( p.lgEOL() )
01507 {
01508
01509 AssertError[nAsserts] = ErrorDefault;
01510 }
01511 }
01512
01513 else if( p.nMatch("NOTH") )
01514 {
01515 strcpy( chAssertType[nAsserts] , "NO" );
01516 strcpy( chAssertLineLabel[nAsserts], "noth" );
01517 wavelength[nAsserts] = 0;
01518 AssertQuantity[nAsserts] = 0.;
01519 AssertError[nAsserts] = ErrorDefault;
01520 }
01521 else
01522 {
01523
01524 fprintf( ioQQQ,
01525 " Unrecognized command. The line image was\n");
01526 p.PrintLine(ioQQQ);
01527 fprintf( ioQQQ,
01528 " The options I know about are: ionization, line, departure coefficient, map, column, "
01529 "temperature, nzone, csupre, htot, itrz, eden, thickness, niter, \n");
01530 fprintf( ioQQQ, " Sorry.\n" );
01531 cdEXIT(EXIT_FAILURE);
01532 }
01533
01534
01535 ++nAsserts;
01536 if( nAsserts >= NASSERTS )
01537 {
01538 fprintf(ioQQQ,
01539 " ParseMonitorResults: too many asserts, limit is NASSERT=%d\n",
01540 NASSERTS );
01541 cdEXIT(EXIT_FAILURE);
01542 }
01543 return;
01544 }
01545
01546
01547
01548 bool lgCheckMonitors(
01549
01550
01551 FILE * ioMONITOR )
01552 {
01553 double PredQuan[NASSERTS] , RelError[NASSERTS];
01554
01555
01556 bool lgFound[NASSERTS];
01557 double relint , absint;
01558 bool lg1OK[NASSERTS];
01559 long i,j;
01560
01561
01562
01563
01564
01565 long ipDisambiguate[NASSERTS][3];
01566 long lgDisambiguate = false;
01567 char chCaps[5], chFind[5];
01568 realnum errorwave;
01569
01570 DEBUG_ENTRY( "lgCheckMonitors()" );
01571
01572
01573
01574 lgMonitorsOK = true;
01575
01576
01577 lgBigBotch = false;
01578
01579
01580
01581
01582
01583
01584
01585 if( !called.lgTalk && optimize.lgOptimize )
01586 {
01587
01588 return true;
01589 }
01590
01591
01592
01593
01594
01595 if( lines_table() )
01596 {
01597 lgBigBotch = true;
01598 lgMonitorsOK = false;
01599 }
01600
01601
01602 for(i=0; i<nAsserts; ++i )
01603 {
01604 lg1OK[i] = true;
01605 PredQuan[i] = 0.;
01606 RelError[i] = 0.;
01607 ipDisambiguate[i][0] = -1;
01608 ipDisambiguate[i][1] = -1;
01609 ipDisambiguate[i][2] = -1;
01610
01611
01612 lgFound[i] = true;
01613
01614
01615
01616 if( strcmp( chAssertType[i] , "Lr")==0 )
01617 {
01618
01619 ipDisambiguate[i][0] = cdLine( chAssertLineLabel[i] , wavelength[i] , &relint,&absint );
01620 if( ipDisambiguate[i][0] <= 0 )
01621 {
01622 fprintf( ioMONITOR,
01623 " monitor error: lgCheckMonitors could not find a line with label %s ",
01624 chAssertLineLabel[i] );
01625 prt_wl( ioMONITOR, wavelength[i] );
01626 fprintf( ioMONITOR, "\n" );
01627
01628 fprintf( ioMONITOR,
01629 " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
01630
01631 lg1OK[i] = false;
01632 RelError[i] = 100000.;
01633 PredQuan[i] = 0;
01634 lg1OK[i] = false;
01635 lgMonitorsOK = false;
01636 lgFound[i] = false;
01637 continue;
01638 }
01639 else
01640 {
01641
01642
01643
01644
01645
01646 cap4(chFind, chAssertLineLabel[i]);
01647
01648
01649 errorwave = WavlenErrorGet( wavelength[i] );
01650
01651
01652 for( j=1; j < LineSave.nsum; j++ )
01653 {
01654
01655 if( j==ipDisambiguate[i][0] )
01656 continue;
01657
01658
01659 cap4(chCaps, LineSv[j].chALab);
01660
01661
01662
01663
01664
01665
01666
01667 if( fabs(LineSv[j].wavelength-wavelength[i]) < 3.f*errorwave )
01668 {
01669
01670 if( strcmp(chCaps,chFind) == 0 )
01671 {
01672 double relint1, absint1, current_error;
01673
01674 cdLine_ip( j, &relint1, &absint1 );
01675 current_error = fabs(1. - relint1/AssertQuantity[i]);
01676
01677 if( current_error < 2.*AssertError[i] ||
01678 current_error < 2.*fabs(RelError[i]) )
01679 {
01680 lgDisambiguate = true;
01681
01682
01683 if( ipDisambiguate[i][1] > 0 )
01684 {
01685 ipDisambiguate[i][2] = j;
01686 break;
01687 }
01688 else
01689 {
01690 ipDisambiguate[i][1] = j;
01691 }
01692 }
01693 }
01694 }
01695 }
01696 }
01697
01698 PredQuan[i] = relint;
01699 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
01700 }
01701
01702
01703 else if( strcmp( chAssertType[i] , "Ll")==0 )
01704 {
01705
01706 if( cdLine( chAssertLineLabel[i] , wavelength[i] , &relint,&absint )<=0 )
01707 {
01708 fprintf( ioMONITOR,
01709 " monitor error: lgCheckMonitors could not find a line with label %s ",
01710 chAssertLineLabel[i] );
01711 prt_wl( ioMONITOR, wavelength[i] );
01712 fprintf( ioMONITOR, "\n" );
01713
01714 fprintf( ioMONITOR,
01715 " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
01716
01717 lg1OK[i] = false;
01718 RelError[i] = 10000000.;
01719 PredQuan[i] = 0;
01720 lg1OK[i] = false;
01721 lgFound[i] = false;
01722 lgMonitorsOK = false;
01723 continue;
01724 }
01725
01726 PredQuan[i] = pow(10.,absint);
01727 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
01728
01729 }
01730 else if( strcmp( chAssertType[i] , "hh" ) == 0)
01731 {
01732 double hfrac , hefrac;
01733
01734 if( cdIonFrac(
01735
01736 "hydr",
01737
01738 1,
01739
01740 &hfrac,
01741
01742 "VOLUME" ,
01743
01744 false) )
01745 {
01746 fprintf( ioMONITOR,
01747 " monitor error: lgCheckMonitors could not find h ionization fraction \n");
01748 lg1OK[i] = false;
01749 RelError[i] = 0;
01750 PredQuan[i] = 0;
01751 lgFound[i] = false;
01752 continue;
01753 }
01754 if( cdIonFrac(
01755
01756 "heli",
01757
01758 1,
01759
01760 &hefrac,
01761
01762 "VOLUME" ,
01763
01764 false) )
01765 {
01766 fprintf( ioMONITOR,
01767 " monitor error: lgCheckMonitors could not find h ionization fraction \n");
01768 lg1OK[i] = false;
01769 RelError[i] = 0;
01770 PredQuan[i] = 0;
01771 lgFound[i] = false;
01772 continue;
01773 }
01774
01775 PredQuan[i] = hefrac-hfrac;
01776
01777 RelError[i] = fabs(AssertQuantity[i] - (hefrac-hfrac) );
01778 }
01779
01780 else if( strcmp( chAssertType[i] , "mp" ) == 0)
01781 {
01782 PredQuan[i] = cpu.lgMPI() ? 1. : 0.;
01783
01784 RelError[i] = AssertQuantity[i] - PredQuan[i];
01785 }
01786
01787 else if( strcmp( chAssertType[i] , "z " ) == 0)
01788 {
01789
01790 PredQuan[i] = (double)nzone;
01791 if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
01792 RelError[i] = ForcePass(chAssertLimit[i]);
01793 else
01794
01795 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
01796 }
01797
01798 else if( strcmp( chAssertType[i] , "or" ) == 0)
01799 {
01800
01801 PredQuan[i] = h2.ortho_density / SDIV( h2.para_density );
01802
01803
01804 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01805 }
01806
01807 else if( strcmp( chAssertType[i] , "g2" ) == 0)
01808 {
01809
01810 PredQuan[i] = gv.rate_h2_form_grains_used_total;
01811
01812 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01813 }
01814
01815 else if( strcmp( chAssertType[i] , "RM" ) == 0)
01816 {
01817
01818 PredQuan[i] = pressure.RadBetaMax;
01819
01820 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01821 }
01822
01823 else if( strcmp( chAssertType[i] , "pr" ) == 0)
01824 {
01825
01826 double sumx=0., sumx2=0., average;
01827 long int n;
01828
01829 for( n=0; n<nzone; n++ )
01830 {
01831 sumx += struc.pressure[n];
01832 sumx2 += POW2(struc.pressure[n]);
01833 }
01834 if( nzone>1 )
01835 {
01836
01837 average = sumx/nzone;
01838
01839 sumx = sqrt( (sumx2-POW2(sumx)/nzone)/(nzone-1) );
01840
01841 PredQuan[i] = sumx / average;
01842 }
01843 else
01844 {
01845 PredQuan[i] = 0.;
01846 }
01847
01848
01849 RelError[i] = PredQuan[i];
01850 }
01851
01852 else if( strcmp( chAssertType[i] , "iz") == 0 )
01853 {
01854
01855 if( nzone > 0 )
01856 PredQuan[i] = (double)(conv.nTotalIoniz-conv.nTotalIoniz_start)/(double)(nzone);
01857 else
01858
01859 PredQuan[i] = 1e10;
01860
01861 if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
01862 RelError[i] = ForcePass(chAssertLimit[i]);
01863 else
01864
01865 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01866 }
01867
01868 else if( strcmp( chAssertType[i] , "e ") == 0 )
01869 {
01870
01871 PredQuan[i] = dense.eden;
01872
01873 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01874 }
01875
01876 else if( strcmp( chAssertType[i] , "Tu") == 0 )
01877 {
01878
01879 PredQuan[i] = pow((rfield.EnergyIncidCont+rfield.EnergyDiffCont)/
01880 SPEEDLIGHT/7.56464e-15,0.25);
01881
01882 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01883 }
01884
01885 else if( strcmp( chAssertType[i] , "th") == 0 )
01886 {
01887
01888 PredQuan[i] = radius.depth;
01889
01890 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01891 }
01892
01893 else if( strcmp( chAssertType[i] , "ra") == 0 )
01894 {
01895
01896 PredQuan[i] = radius.Radius;
01897
01898 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01899 }
01900
01901 else if( strcmp( chAssertType[i] , "Ve") == 0 )
01902 {
01903
01904 PredQuan[i] = wind.windv/1e5;
01905
01906 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01907 }
01908
01909 else if( strcmp( chAssertType[i] , "NO") == 0 )
01910 {
01911
01912 PredQuan[i] = 0;
01913
01914 RelError[i] = 0.;
01915 }
01916
01917 else if( strcmp( chAssertType[i] , "sc") == 0 )
01918 {
01919
01920 PredQuan[i] = secondaries.csupra[ipHYDROGEN][0];
01921
01922 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01923 }
01924
01925 else if( strcmp( chAssertType[i] , "ht") == 0 )
01926 {
01927
01928 PredQuan[i] = thermal.htot;
01929
01930 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i];
01931 }
01932
01933 else if( strcmp( chAssertType[i] , "Z ") == 0 )
01934 {
01935
01936 PredQuan[i] = (double)iteration;
01937 if( t_version::Inst().lgRelease || t_version::Inst().lgReleaseBranch )
01938 RelError[i] = ForcePass(chAssertLimit[i]);
01939 else
01940
01941 RelError[i] = (double)(AssertQuantity[i] - PredQuan[i]);
01942 }
01943
01944 else if( strcmp( chAssertType[i] , "CB") == 0 )
01945 {
01946 long int nISOCaseB = (long)Param[i][0];
01947 long int nelemCaseB = (long)Param[i][1];
01948 char chElemLabelCaseB[5];
01949
01950 strcpy( chElemLabelCaseB , elementnames.chElementSym[nelemCaseB] );
01951 char chNumb[4];
01952 sprintf( chNumb , "%2li" , nelemCaseB+1-nISOCaseB );
01953 strcat( chElemLabelCaseB , chNumb );
01954
01955
01956 int iCase;
01957 if( strcmp( "Ca A", chAssertLineLabel[i] ) == 0 )
01958 iCase = 0;
01959 else if( strcmp( "Ca B", chAssertLineLabel[i] ) == 0 )
01960 iCase = 1;
01961 else if( strcmp( "+Col", chAssertLineLabel[i] ) == 0 )
01962 iCase = 1;
01963 else
01964 TotalInsanity();
01965
01966 RelError[i] = 0.;
01967 long nHighestPrinted = iso.n_HighestResolved_max[nISOCaseB][nelemCaseB];
01968 if( nISOCaseB == ipH_LIKE )
01969 {
01970 fprintf(ioMONITOR," Species nHi nLo Wl Computed Asserted error\n");
01971
01972
01973
01974
01975 for( long int ipLo=1+iCase; ipLo< MIN2(10,nHighestPrinted-1); ++ipLo )
01976 {
01977 for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi )
01978 {
01979
01980 realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo];
01981
01982 if( wl < Param[i][2] || wl > Param[i][3] )
01983 continue;
01984
01985 double relint , absint,CBrelint , CBabsint;
01986
01987 cdLine( chAssertLineLabel[i] , wl , &CBrelint , &CBabsint );
01988 if( CBrelint < Param[i][4] )
01989 continue;
01990 CBabsint = pow( 10., CBabsint );
01991 double error;
01992
01993 if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0)
01994 {
01995 absint = pow( 10., absint );
01996 error = (CBabsint - absint)/MAX2(CBabsint , absint);
01997 double RelativeError = fabs(error) / AssertError[i];
01998
01999 if( RelativeError < 1. )
02000 {
02001 if( RelativeError < 0.25 )
02002 {
02003 fprintf( ioMONITOR, " ChkMonitor ");
02004 }
02005 else if( RelativeError < 0.50 )
02006 {
02007 fprintf( ioMONITOR, " ChkMonitor - ");
02008 }
02009 else if( RelativeError < 0.75 )
02010 {
02011 fprintf( ioMONITOR, " ChkMonitor -- ");
02012 }
02013 else if( RelativeError < 0.90 )
02014 {
02015 fprintf( ioMONITOR, " ChkMonitor --- ");
02016 }
02017 else if( RelativeError < 0.95 )
02018 {
02019 fprintf( ioMONITOR, " ChkMonitor ---- ");
02020 }
02021 else if( RelativeError < 0.98 )
02022 {
02023 fprintf( ioMONITOR, " ChkMonitor ----- ");
02024 }
02025 else
02026 {
02027 fprintf( ioMONITOR, " ChkMonitor ------ ");
02028 }
02029
02030 }
02031 else
02032 {
02033 fprintf( ioMONITOR, " ChkMonitor botch>>");
02034 }
02035 fprintf(ioMONITOR," %s %3li %3li ",
02036 chElemLabelCaseB , ipHi , ipLo );
02037 prt_wl(ioMONITOR, wl );
02038 fprintf(ioMONITOR," %.2e %.2e %10.3f",
02039 absint , CBabsint , error );
02040 }
02041 else
02042 TotalInsanity();
02043 if( fabs(error) > AssertError[i] )
02044 fprintf(ioMONITOR , " botch \n");
02045 else
02046 fprintf(ioMONITOR , "\n");
02047
02048 PredQuan[i] = 0;
02049 AssertQuantity[i] = 0.;
02050 RelError[i] = MAX2( RelError[i] , fabs(error) );
02051
02052
02053 MonitorResults.SumErrorCaseMonitor += RelError[i];
02054 ++MonitorResults.nSumErrorCaseMonitor;
02055
02056 }
02057 }
02058 fprintf(ioMONITOR,"\n");
02059 }
02060 else if( nISOCaseB == ipHE_LIKE )
02061 {
02062 if( !dense.lgElmtOn[ipHELIUM] )
02063 {
02064 fprintf(ioQQQ,"PROBLEM monitor case B for a He is requested but He is not "
02065 "included.\n");
02066 fprintf(ioQQQ,"Do not turn off He if you want to monitor its spectrum.\n");
02067 cdEXIT(EXIT_FAILURE);
02068 }
02069 # if 0
02070 # define N_CASEB_HEI 11
02071 realnum CaseBWlHeI[N_CASEB_HEI]=
02072 { 10830.f, 3889.f, 3188.f, 5016.f, 3965.f, 7065.f, 5876.f, 4471.f,
02073 4026.f, 6678.f, 4922.f };
02074 # endif
02075
02076 fprintf(ioMONITOR," Wl Computed Asserted error\n");
02077 for( long int ipLine=0; ipLine< atmdat.nCaseBHeI ; ++ipLine )
02078 {
02079
02080 realnum wl = atmdat.CaseBWlHeI[ipLine];
02081
02082 if( wl < Param[i][2] || wl > Param[i][3] )
02083 continue;
02084 double relint , absint,CBrelint , CBabsint;
02085 cdLine( chAssertLineLabel[i] , wl , &CBrelint , &CBabsint );
02086 if( CBrelint < Param[i][4] )
02087 continue;
02088 CBabsint = pow( 10., CBabsint );
02089 double error;
02090 if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0)
02091 {
02092 absint = pow( 10., absint );
02093 error = (CBabsint - absint)/MAX2(CBabsint , absint);
02094 double RelativeError = fabs(error) / AssertError[i];
02095
02096 if( RelativeError < 1. )
02097 {
02098 if( RelativeError < 0.25 )
02099 {
02100 fprintf( ioMONITOR, " ChkMonitor ");
02101 }
02102 else if( RelativeError < 0.50 )
02103 {
02104 fprintf( ioMONITOR, " ChkMonitor - ");
02105 }
02106 else if( RelativeError < 0.75 )
02107 {
02108 fprintf( ioMONITOR, " ChkMonitor -- ");
02109 }
02110 else if( RelativeError < 0.90 )
02111 {
02112 fprintf( ioMONITOR, " ChkMonitor --- ");
02113 }
02114 else if( RelativeError < 0.95 )
02115 {
02116 fprintf( ioMONITOR, " ChkMonitor ---- ");
02117 }
02118 else if( RelativeError < 0.98 )
02119 {
02120 fprintf( ioMONITOR, " ChkMonitor ----- ");
02121 }
02122 else
02123 {
02124 fprintf( ioMONITOR, " ChkMonitor ------ ");
02125 }
02126
02127 }
02128 else
02129 {
02130 fprintf( ioMONITOR, " ChkMonitor botch>>");
02131 }
02132 prt_wl(ioMONITOR, wl );
02133 fprintf(ioMONITOR," %.2e %.2e %10.3f",
02134 absint , CBabsint , error );
02135 }
02136 else
02137 TotalInsanity();
02138 if( fabs(error) > AssertError[i] )
02139 fprintf(ioMONITOR , " botch \n");
02140 else
02141 fprintf(ioMONITOR , "\n");
02142
02143 PredQuan[i] = 0;
02144 AssertQuantity[i] = 0.;
02145 RelError[i] = MAX2( RelError[i] , fabs(error) );
02146
02147
02148 MonitorResults.SumErrorCaseMonitor += RelError[i];
02149 ++MonitorResults.nSumErrorCaseMonitor;
02150 }
02151 fprintf(ioMONITOR,"\n");
02152 }
02153 else
02154 TotalInsanity();
02155 }
02156
02157 else if( strcmp( chAssertType[i] , "D ") == 0 )
02158 {
02159 long ipZ , n;
02160
02161
02162 ipZ = (long)wavelength[i];
02163 if( !dense.lgElmtOn[ipZ] )
02164 {
02165 fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ);
02166 PredQuan[i] = 0;
02167 RelError[i] = 0.;
02168 }
02169 else
02170 {
02171 RelError[i] = 0.;
02172 PredQuan[i] = 0.;
02173 long numPrintLevels = iso.numLevels_local[ipH_LIKE][ipZ] - (long)AssertQuantity2[i];
02174 for( n=(long)AssertQuantity2[i]; n<numPrintLevels+(long)AssertQuantity2[i]; ++n )
02175 {
02176 PredQuan[i] += iso.DepartCoef[ipH_LIKE][ipZ][n];
02177
02178 RelError[i] = MAX2( RelError[i] , iso.DepartCoef[ipH_LIKE][ipZ][n] -1.);
02179 RelError[i] = MAX2( RelError[i] , 1. - iso.DepartCoef[ipH_LIKE][ipZ][n] );
02180 }
02181 PredQuan[i] /= (double)(numPrintLevels);
02182 }
02183 }
02184
02185
02186 else if( strcmp( chAssertType[i] , "De") == 0 )
02187 {
02188 long ipZ , n;
02189
02190
02191 ipZ = (long)wavelength[i];
02192 if( !dense.lgElmtOn[ipZ] )
02193 {
02194 fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ);
02195 PredQuan[i] = 0.;
02196 RelError[i] = 0.;
02197 }
02198 else
02199 {
02200 RelError[i] = 0.;
02201 PredQuan[i] = 0.;
02202 long numPrintLevels = iso.numLevels_local[ipHE_LIKE][ipZ] - (long)AssertQuantity2[i];
02203 for( n=(long)AssertQuantity2[i]; n<numPrintLevels+(long)AssertQuantity2[i]; ++n )
02204 {
02205 double relerror;
02206 relerror = 0.;
02207 PredQuan[i] += iso.DepartCoef[ipHE_LIKE][ipZ][n];
02208
02209 relerror = iso.DepartCoef[ipHE_LIKE][ipZ][n] -1.;
02210 relerror = MAX2( relerror , 1. - iso.DepartCoef[ipHE_LIKE][ipZ][n] );
02211 RelError[i] = MAX2( RelError[i] , relerror / AssertQuantity[i] );
02212 }
02213 PredQuan[i] /= (double)(numPrintLevels);
02214 }
02215 }
02216
02217
02218 else if( strcmp( chAssertType[i] , "d ") == 0 )
02219 {
02220 double BigError , StdDev;
02221
02222 AssertFeIIDep( &PredQuan[i] , &BigError , &StdDev );
02223
02224
02225
02226 RelError[i] = StdDev;
02227 }
02228
02229
02230 else if( strcmp( chAssertType[i] , "d-") == 0 )
02231 {
02232 PredQuan[i] = hmi.hmidep;
02233 RelError[i] = fabs( hmi.hmidep - 1.);
02234 }
02235
02236
02237 else if( (strcmp( chAssertType[i] , "f ") == 0) ||
02238 (strcmp( chAssertType[i] , "F ") == 0) )
02239 {
02240 char chWeight[7];
02241 if( strcmp( chAssertType[i] , "F ") == 0 )
02242 {
02243 strcpy( chWeight , "VOLUME" );
02244 }
02245 else
02246 {
02247
02248 strcpy( chWeight , "RADIUS" );
02249 }
02250
02251 if( cdIonFrac(
02252
02253 chAssertLineLabel[i],
02254
02255 (long)wavelength[i],
02256
02257 &relint,
02258
02259 chWeight ,
02260
02261 false) )
02262 {
02263 fprintf( ioMONITOR,
02264 " monitor error: lgCheckMonitors could not find a line with label %s %f \n",
02265 chAssertLineLabel[i] , wavelength[i] );
02266
02267 lg1OK[i] = false;
02268 RelError[i] = 0;
02269 PredQuan[i] = 0;
02270 lgFound[i] = false;
02271 continue;
02272 }
02273
02274 PredQuan[i] = relint;
02275 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
02276 }
02277
02278
02279 else if( strcmp( chAssertType[i] , "cd") == 0)
02280 {
02281
02282 if( (strcmp( chAssertLineLabel[i], "H2 " ) == 0) && (Param[i][0] >= 0.) )
02283 {
02284
02285
02286 if( (relint = cdH2_colden( (long)Param[i][0] , (long)Param[i][1] ) ) < 0. )
02287 {
02288 fprintf(ioQQQ," PROBLEM lgCheckMonitors did not find v=%li, J=%li for H2 column density.\n",
02289 (long)Param[i][0] , (long)Param[i][1] );
02290 lg1OK[i] = false;
02291 RelError[i] = 0;
02292 PredQuan[i] = 0;
02293 lgFound[i] = false;
02294 continue;
02295 }
02296 }
02297 else
02298 {
02299
02300 if( cdColm(
02301
02302 chAssertLineLabel[i],
02303
02304
02305 (long)wavelength[i],
02306
02307 &relint) )
02308 {
02309 fprintf( ioMONITOR,
02310 " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
02311 chAssertLineLabel[i] , wavelength[i] );
02312
02313 lg1OK[i] = false;
02314 RelError[i] = 0;
02315 PredQuan[i] = 0;
02316 lgFound[i] = false;
02317 continue;
02318 }
02319 }
02320
02321 PredQuan[i] = relint;
02322 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
02323 }
02324
02325
02326 else if( (strcmp( chAssertType[i] , "MF") == 0) || (strcmp( chAssertType[i] , "mf") == 0) )
02327 {
02328
02329 relint = 0.;
02330 if( strcmp( chAssertLineLabel[i], "H2 ") ==0)
02331 {
02332
02333 if( cdColm("H2 " , 0,
02334
02335 &relint) )
02336 TotalInsanity();
02337
02338 relint = relint / colden.colden[ipCOL_HTOT];
02339 }
02340 else
02341 {
02342 fprintf( ioMONITOR,
02343 " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
02344 chAssertLineLabel[i] , wavelength[i] );
02345
02346 lg1OK[i] = false;
02347 RelError[i] = 0;
02348 PredQuan[i] = 0;
02349 continue;
02350 }
02351
02352 PredQuan[i] = relint;
02353 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
02354 }
02355
02356
02357 else if( strcmp( chAssertType[i] , "mh") == 0 || strcmp( chAssertType[i] , "mc") == 0)
02358 {
02359
02360
02361 if( hcmap.nMapAlloc == 0 )
02362 {
02363
02364 fprintf( ioMONITOR,
02365 " monitor error: lgCheckMonitors cannot check map since map not done.\n");
02366
02367 lg1OK[i] = false;
02368 RelError[i] = 0;
02369 PredQuan[i] = 0;
02370 continue;
02371 }
02372
02373 if( AssertQuantity2[i]< hcmap.temap[0] || AssertQuantity2[i]> hcmap.temap[hcmap.nmap-1] )
02374 {
02375 fprintf( ioMONITOR,
02376 " monitor error: lgCheckMonitors cannot check map since temperature not within range.\n");
02377
02378 lg1OK[i] = false;
02379 RelError[i] = 0;
02380 PredQuan[i] = 0;
02381 continue;
02382 }
02383
02384
02385 j = 0;
02386 while( AssertQuantity2[i]>hcmap.temap[j]*1.001 && j < hcmap.nmap )
02387 {
02388 ++j;
02389 }
02390
02391
02392
02393 if( strcmp( chAssertType[i] , "mh") == 0 )
02394 {
02395
02396 PredQuan[i] = hcmap.hmap[j];
02397 strcpy(chAssertLineLabel[i],"MapH" );
02398 }
02399 else if( strcmp( chAssertType[i] , "mc") == 0)
02400 {
02401
02402 PredQuan[i] = hcmap.cmap[j];
02403 strcpy(chAssertLineLabel[i],"MapC" );
02404 }
02405 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
02406 }
02407
02408
02409 else if( (strcmp( chAssertType[i] , "t ") == 0) ||
02410 (strcmp( chAssertType[i] , "T ") == 0) )
02411 {
02412 char chWeight[7];
02413 if( strcmp( chAssertType[i] , "T ") == 0 )
02414 {
02415 strcpy( chWeight , "VOLUME" );
02416 }
02417 else
02418 {
02419
02420 strcpy( chWeight , "RADIUS" );
02421 }
02422
02423
02424 if( strcmp( chAssertLineLabel[i], "GTem" ) == 0 )
02425 {
02426 size_t nd;
02427
02428
02429 nd = (size_t)wavelength[i]-1;
02430 if( nd >= gv.bin.size() )
02431 {
02432 fprintf( ioQQQ, "Illegal grain number found: %f\n" , wavelength[i] );
02433 fprintf( ioQQQ, "Use 1 for first grain that is turned on, " );
02434 fprintf( ioQQQ, "2 for second, etc....\n" );
02435 fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" );
02436 cdEXIT(EXIT_FAILURE);
02437 }
02438 relint = gv.bin[nd]->avdust/radius.depth_x_fillfac;
02439 }
02440 else if( strcmp( chAssertLineLabel[i], "face" ) == 0 )
02441 {
02442
02443 relint = struc.testr[0];
02444 }
02445 else
02446 {
02447
02448 if( cdTemp(
02449
02450 chAssertLineLabel[i],
02451
02452 (long)wavelength[i],
02453
02454 &relint,
02455
02456 chWeight ) )
02457 {
02458 fprintf( ioMONITOR,
02459 " monitor error: lgCheckMonitors could not find an ion with label %s ion %li \n",
02460 chAssertLineLabel[i] , (long)wavelength[i] );
02461
02462 lg1OK[i] = false;
02463 RelError[i] = 0;
02464 PredQuan[i] = 0;
02465 lgFound[i] = false;
02466 continue;
02467 }
02468 }
02469
02470 PredQuan[i] = relint;
02471 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i];
02472 }
02473
02474
02475 else if( strcmp( chAssertType[i], "gp" ) == 0 )
02476 {
02477
02478
02479 size_t nd = (size_t)wavelength[i]-1;
02480 if( nd >= gv.bin.size() )
02481 {
02482 fprintf( ioQQQ, "Illegal grain number found: %g\n" , wavelength[i] );
02483 fprintf( ioQQQ, "Use 1 for first grain that is turned on, " );
02484 fprintf( ioQQQ, "2 for second, etc....\n" );
02485 fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" );
02486 cdEXIT(EXIT_FAILURE);
02487 }
02488
02489
02490 PredQuan[i] = gv.bin[nd]->avdpot/radius.depth_x_fillfac;
02491
02492 RelError[i] = AssertQuantity[i] - PredQuan[i];
02493 }
02494
02495 else
02496 {
02497 fprintf( ioMONITOR,
02498 " monitor error: lgCheckMonitors received an insane chAssertType=%s, impossible\n",
02499 chAssertType[i] );
02500 ShowMe();
02501 cdEXIT(EXIT_FAILURE);
02502 }
02503
02504 if( chAssertLimit[i] == '=' )
02505 {
02506
02507 if( fabs(RelError[i]) > AssertError[i] )
02508 {
02509 lg1OK[i] = false;
02510 lgMonitorsOK = false;
02511 }
02512 }
02513 else if( chAssertLimit[i] == '<' )
02514 {
02515
02516
02517
02518
02519
02520 if( RelError[i] <= 0. )
02521 {
02522 lg1OK[i] = false;
02523 lgMonitorsOK = false;
02524 }
02525 }
02526 else if( chAssertLimit[i] == '>' )
02527 {
02528
02529
02530 if( RelError[i] >= 0. )
02531 {
02532 lg1OK[i] = false;
02533 lgMonitorsOK = false;
02534 }
02535 }
02536 }
02537
02538
02539 if( called.lgTalk && nAsserts>0 )
02540 {
02541 time_t now;
02542
02543
02544 if( lgDisambiguate )
02545 {
02546
02547 long sigfigsav = LineSave.sig_figs;
02548 double relint1, relint2, absint1;
02549
02550 LineSave.sig_figs = 6;
02551
02552 fprintf( ioMONITOR, "=============Line Disambiguation=======================================================\n" );
02553 fprintf( ioMONITOR, " Wavelengths || Intensities \n" );
02554 fprintf( ioMONITOR, "Label line match1 match2 match3 || asserted match1 match2 match3\n" );
02555
02556 for( i=0; i<nAsserts; ++i )
02557 {
02558 if( ipDisambiguate[i][1] > 0 )
02559 {
02560 fprintf( ioMONITOR , "%4s ", chAssertLineLabel[i] );
02561 prt_wl( ioMONITOR , wavelength[i] );
02562 fprintf( ioMONITOR , " " );
02563 prt_wl( ioMONITOR , LineSv[ipDisambiguate[i][0]].wavelength );
02564 fprintf( ioMONITOR , " " );
02565 prt_wl( ioMONITOR , LineSv[ipDisambiguate[i][1]].wavelength );
02566 fprintf( ioMONITOR , " " );
02567 if( ipDisambiguate[i][2] > 0 )
02568 {
02569 prt_wl( ioMONITOR , LineSv[ipDisambiguate[i][2]].wavelength );
02570 cdLine_ip( ipDisambiguate[i][2], &relint2, &absint1 );
02571 }
02572 else
02573 {
02574 fprintf( ioMONITOR , "--------" );
02575 relint2 = 0.0;
02576 }
02577 fprintf( ioMONITOR , " ||" );
02578
02579 cdLine_ip( ipDisambiguate[i][1], &relint1, &absint1 );
02580
02581 if( lgPrtSciNot )
02582 {
02583 fprintf( ioMONITOR , " %10.3e %10.3e %10.3e %10.3e\n",
02584 AssertQuantity[i],
02585 PredQuan[i] ,
02586 relint1,
02587 relint2 );
02588 }
02589 else
02590 {
02591 fprintf( ioMONITOR , " %10.4f %10.4f %10.4f %10.4f\n",
02592 AssertQuantity[i],
02593 PredQuan[i] ,
02594 relint1,
02595 relint2 );
02596 }
02597 }
02598 }
02599 fprintf( ioMONITOR, "\n" );
02600
02601
02602 LineSave.sig_figs = sigfigsav;
02603 }
02604
02605
02606 fprintf( ioMONITOR, "=============Results of monitors: Cloudy %s ", t_version::Inst().chVersion );
02607
02608
02609
02610 if( prt.lgPrintTime )
02611 {
02612
02613 now = time(NULL);
02614
02615 fprintf(ioMONITOR,"%s", ctime(&now) );
02616 }
02617 else
02618 {
02619 fprintf(ioMONITOR,"\n" );
02620 }
02621
02622 if( lgMonitorsOK )
02623 {
02624 fprintf( ioMONITOR, " No errors were found. Summary follows.\n");
02625 }
02626 else
02627 {
02628 fprintf( ioMONITOR, " Errors were found. Summary follows.\n");
02629 }
02630
02631 fprintf( ioMONITOR,
02632 " Label line computed asserted Rel Err Set err\n");
02633
02634 for( i=0; i<nAsserts; ++i )
02635 {
02636 double prtPredQuan, prtAssertQuantity;
02637
02638
02639 if( lgQuantityLog[i] )
02640 {
02641 prtPredQuan = log10( MAX2( SMALLDOUBLE , PredQuan[i] ) );
02642 prtAssertQuantity = log10( MAX2( SMALLDOUBLE , AssertQuantity[i] ) );
02643 }
02644 else
02645 {
02646 prtPredQuan = PredQuan[i];
02647 prtAssertQuantity = AssertQuantity[i];
02648 }
02649
02650 if( lg1OK[i] )
02651 {
02652
02653
02654 double relative = fabs(RelError[i] / SDIV( AssertError[i]));
02655
02656 if( relative < 0.25 || chAssertLimit[i] != '=' )
02657 {
02658 fprintf( ioMONITOR, " ChkMonitor ");
02659 }
02660 else if( relative < 0.50 )
02661 {
02662 fprintf( ioMONITOR, " ChkMonitor - ");
02663 }
02664 else if( relative < 0.75 )
02665 {
02666 fprintf( ioMONITOR, " ChkMonitor -- ");
02667 }
02668 else if( relative < 0.90 )
02669 {
02670 fprintf( ioMONITOR, " ChkMonitor --- ");
02671 }
02672 else if( relative < 0.95 )
02673 {
02674 fprintf( ioMONITOR, " ChkMonitor ---- ");
02675 }
02676 else if( relative < 0.98 )
02677 {
02678 fprintf( ioMONITOR, " ChkMonitor ----- ");
02679 }
02680 else
02681 {
02682 fprintf( ioMONITOR, " ChkMonitor ------ ");
02683 }
02684
02685 }
02686 else
02687 {
02688 fprintf( ioMONITOR, " ChkMonitor botch>>");
02689 }
02690 if( strcmp( chAssertType[i] , "Ll")==0 || ( strcmp( chAssertType[i] , "Lr")==0 ) )
02691 {
02692
02693 fprintf( ioMONITOR , "%4s ",
02694 chAssertLineLabel[i] );
02695
02696 prt_wl( ioMONITOR , wavelength[i] );
02697
02698 if( lgPrtSciNot )
02699 {
02700 fprintf( ioMONITOR , " %10.3e %c %10.3e %7.3f %7.3f ",
02701 prtPredQuan ,
02702 chAssertLimit[i] ,
02703 prtAssertQuantity ,
02704 RelError[i] ,
02705 AssertError[i]);
02706 }
02707 else
02708 {
02709 fprintf( ioMONITOR , " %10.4f %c %10.4f %7.3f %7.3f ",
02710 prtPredQuan ,
02711 chAssertLimit[i] ,
02712 prtAssertQuantity ,
02713 RelError[i] ,
02714 AssertError[i]);
02715 }
02716
02717 {
02718 enum {DEBUG_LOC=false};
02719
02720 if( DEBUG_LOC )
02721 {
02722 long ipHi, ipLo;
02723 errorwave = WavlenErrorGet( wavelength[i] );
02724
02725 for( ipLo = ipHe1s1S; ipLo <= iso.numLevels_local[ipHE_LIKE][ipHELIUM] -
02726 iso.nCollapsed_local[ipHE_LIKE][ipHELIUM]; ipLo ++ )
02727 {
02728 for( ipHi = ipLo+1; ipHi < iso.numLevels_local[ipHE_LIKE][ipHELIUM] -
02729 iso.nCollapsed_local[ipHE_LIKE][ipHELIUM]; ipHi++ )
02730 {
02731 if( fabs(Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng - wavelength[i])
02732 < errorwave && ipLo!=ipHe2p3P0 && ipLo!=ipHe2p3P1 )
02733 {
02734 fprintf( ioQQQ, "\t%li %li %li %li %li %li",
02735 StatesElemNEW[ipHELIUM][ipHELIUM-ipHE_LIKE][ipHi].n,
02736 StatesElemNEW[ipHELIUM][ipHELIUM-ipHE_LIKE][ipHi].l,
02737 StatesElemNEW[ipHELIUM][ipHELIUM-ipHE_LIKE][ipHi].S,
02738 StatesElemNEW[ipHELIUM][ipHELIUM-ipHE_LIKE][ipLo].n,
02739 StatesElemNEW[ipHELIUM][ipHELIUM-ipHE_LIKE][ipLo].l,
02740 StatesElemNEW[ipHELIUM][ipHELIUM-ipHE_LIKE][ipLo].S );
02741 break;
02742 }
02743 }
02744 }
02745 }
02746 }
02747
02748 }
02749 else
02750 {
02751 fprintf( ioMONITOR , "%4s %6i %10.4f %c %10.4f %7.3f %7.3f ",
02752 chAssertLineLabel[i] ,
02753 (int)wavelength[i] ,
02754 prtPredQuan ,
02755 chAssertLimit[i] ,
02756 prtAssertQuantity ,
02757 RelError[i] ,
02758 AssertError[i]);
02759 }
02760
02761
02762
02763
02764
02765
02766 if( !lg1OK[i] && (fabs(RelError[i]) > 3.*AssertError[i]) && lgFound[i] )
02767 {
02768 fprintf( ioMONITOR , " <<BIG BOTCH!!\n");
02769 lgBigBotch = true;
02770 }
02771 else
02772 {
02773 fprintf( ioMONITOR , "\n");
02774 }
02775 }
02776 fprintf( ioMONITOR , " \n");
02777
02778
02779
02780 if( !lgMonitorsOK && lgBigBotch )
02781 {
02782
02783 fprintf( ioMONITOR, " BIG BOTCHED MONITORS!!! Big Botched Monitors!!! \n");
02784 }
02785 else if( !lgMonitorsOK )
02786 {
02787 fprintf( ioMONITOR, " BOTCHED MONITORS!!! Botched Monitors!!! \n");
02788 }
02789
02790 if( MonitorResults.nSumErrorCaseMonitor>0 )
02791 {
02792 fprintf(ioMONITOR,"\n The mean of the %li monitor Case A, B relative "
02793 "residuals is %.2f\n\n" ,
02794 MonitorResults.nSumErrorCaseMonitor,
02795 MonitorResults.SumErrorCaseMonitor /MonitorResults.nSumErrorCaseMonitor );
02796 }
02797
02798
02799 if( prt.lgPrintTime )
02800 {
02801 fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo );
02802 }
02803 }
02804 return lgMonitorsOK;
02805 }