87 if( chAssertLimit1 ==
'=' )
89 else if( chAssertLimit1 ==
'<' )
91 else if( chAssertLimit1 ==
'>' )
125 fprintf(
ioQQQ,
" ParseMonitorResults called before InitAsserResults\n" );
189 vector<double> AssertVector;
203 fprintf(
ioQQQ,
"PROBLEM the file %s does not have enough values. sorry\n", chLabel.c_str() );
237 fprintf(
ioQQQ,
" The default monitor error is being changed after"
238 " some asserts were entered. \n This will only affect asserts"
239 " that come after this command.\n");
262 " I could not identify an option on this ASSERT SET XXX command.\n");
269 else if( p.
nMatch(
"IONI" ) )
288 if( (nelem = p.
GetElem()) < 0 )
291 " I could not identify an element name on this line.\n");
305 p.
NoNumb(
"ionization stage");
311 " The ionization stage is inappropriate for this element.\n");
326 p.
NoNumb(
"ionization fraction");
347 " The ionization fraction must be less than one.\n");
356 " The ionization ionization fraction is too small, or zero. Check input\n");
387 else if( p.
nMatch(
" CO " ) )
394 " I could not identify CO or H2 on this line.\n");
406 p.
NoNumb(
"molecular fraction");
428 else if( p.
nMatch(
" LINE " ) )
475 if(0 && b !=
blends.end())
479 else if( chLabel.length() >
NCHLAB-1 )
481 fprintf(
ioQQQ,
" The label must be no more than %d char long, between double quotes.\n",
498 p.
NoNumb(
"intensity/luminosity");
507 " The asserted quantity is a log, but is too large or small, value is %e.\n",
509 fprintf(
ioQQQ,
" I would crash if I tried to use it.\n Sorry.\n" );
518 " The relative intensity must be non-negative, and was %e.\n",
AssertQuantity[nAsserts] );
533 else if( p.
nMatch(
"CASE" ) )
547 if( p.
GetParam(
"FAINT",&Param[nAsserts][4]) )
552 fprintf(
ioQQQ,
" The monitor Case B faint option must have a number,"
553 " the relative intensity of the fainest line to monitor.\n");
557 if( Param[nAsserts][4]<=0. )
567 if( p.
GetRange(
"RANG",&Param[nAsserts][2],&Param[nAsserts][3]) )
572 fprintf(
ioQQQ,
" The monitor Case B range option must have two numbers,"
573 " the lower and upper limit to the wavelengths in Angstroms.\n");
574 fprintf(
ioQQQ,
" There must be a total of three numbers on the line,"
575 " the relative error followed by the lower and upper limits to the "
576 "wavelength in Angstroms.\n");
579 if( Param[nAsserts][2]>Param[nAsserts][3])
597 if( (nelem = p.
GetElem()) < 0 )
600 fprintf(
ioQQQ,
"monitor H-like case B did not find an element on this line, sorry\n");
607 fprintf(
ioQQQ,
"monitor H-like cannot do elements heavier than O, sorry\n");
615 if( p.
nMatch(
"CASE A " ) )
620 else if( p.
nMatch(
"HE-L") )
632 " I could not identify an iso-sequence on this Case A/B command.\n");
639 else if( p.
nMatch(
"DEPA") )
649 for (
long i=chLabel.length();i<
NCHLAB-1;++i)
650 chAssertLineLabel[nAsserts] +=
' ';
656 p.
NoNumb(
"average departure coefficient");
670 else if( p.
nMatch(
"H-LI" ) )
680 " I could not identify an element name on this line.\n");
685 Param[nAsserts][1] = 1.;
691 else if( p.
nMatch(
"HE-L" ) )
701 " I could not identify an element name on this line.\n");
706 Param[nAsserts][1] = 1.;
712 else if( p.
nMatch(
"HMIN" ) )
724 " There must be a second key: H-LIke, HMINus, or HE-Like followed by element.\n");
743 else if( p.
nMatch(
" MAP" ) )
752 else if( p.
nMatch(
"COOL" ) )
759 " There must be a second key, HEATing or COOLing.\n");
785 p.
NoNumb(
"heating/cooling");
805 else if( p.
nMatch(
"COLU" ) )
821 for (
long i=chLabel.length();i<
NCHLAB-1;++i)
868 else if( p.
nMatch(
" CO "))
872 else if( p.
nMatch(
"SIO "))
876 else if( p.
nMatch(
" OH "))
880 else if( p.
nMatch(
" CN ") )
884 else if( p.
nMatch(
" CH ") )
888 else if( p.
nMatch(
" CH+") )
895 " I could not identify H2, H3+, H2+, H2g, H2*, H2H+, CO, O2, SiO, OH, C2 or C3 or on this line.\n");
904 fprintf(
ioQQQ,
" LEVEL option only implemented for H2.\n");
931 p.
NoNumb(
"column density");
969 " I did not find a 2, as in H2, as the first number on this line.\n");
971 " The rate should be the second number.\n");
979 p.
NoNumb(
"grain H2 formation");
1014 p.
NoNumb(
"grain potential");
1029 else if( p.
nMatch(
"TEMP") )
1062 else if( p.
nMatch(
"FACE") )
1070 else if( p.
nMatch(
" H2 ") )
1079 else if( (nelem = p.
GetElem()) >= 0 )
1089 " I could not identify an element name on this line.\n");
1115 p.
NoNumb(
"ionization stage");
1121 " The ionization stage is inappropriate for this element.\n");
1157 else if( p.
nMatch(
"HHEI") )
1169 p.
NoNumb(
"ionization correction factor");
1185 else if( p.
nMatch(
" H2 ") )
1198 " I could not identify a second keyword on this line.\n");
1205 if( p.
lgEOL() || n2 != 2 )
1207 p.
NoNumb(
"the 2 in H2 ?!");
1224 else if( p.
nMatch(
" MPI") )
1237 else if( p.
nMatch(
"NZON") )
1270 p.
NoNumb(
"pressure error");
1288 else if( p.
nMatch(
"PRADMAX") )
1321 else if( p.
nMatch(
"CSUP") )
1332 p.
NoNumb(
"secondary ionization rate");
1354 else if( p.
nMatch(
"HTOT") )
1366 p.
NoNumb(
"heating rate");
1388 else if( p.
nMatch(
"CTOT") )
1407 p.
NoNumb(
"cooling rate");
1430 else if( p.
nMatch(
"ITRZ") )
1441 p.
NoNumb(
"iterations per zone");
1453 else if( p.
nMatch(
"EDEN") )
1465 p.
NoNumb(
" electron density of the last zone");
1483 else if( p.
nMatch(
" TU ") )
1492 p.
NoNumb(
"energy density of last zone");
1526 p.
NoNumb(
"thickness or depth of model");
1544 else if( p.
nMatch(
"RADI") )
1555 p.
NoNumb(
"outer radius");
1573 else if( p.
nMatch(
"NITE") )
1586 p.
NoNumb(
"number of iterations");
1601 else if( p.
nMatch(
"VELO") )
1612 p.
NoNumb(
"terminal velocity");
1624 else if( p.
nMatch(
"NOTH") )
1636 " Unrecognized command. The line image was\n");
1639 " The options I know about are: ionization, line, departure coefficient, map, column, "
1640 "temperature, nzone, csupre, htot, itrz, eden, thickness, niter, \n");
1650 " ParseMonitorResults: too many asserts, limit is NASSERT=%d\n",
1659 return 1. -
safe_div(pred,assert,1.);
1673 double relint , absint;
1683 long lgDisambiguate =
false;
1725 ipDisambiguate[i][j] = -1;
1736 if( ipDisambiguate[i][0] <= 0 )
1738 fprintf( ioMONITOR,
" monitor error: lgCheckMonitors could not find line ");
1742 " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1745 RelError[i] = 100000.;
1769 if( j==ipDisambiguate[i][0] )
1784 if( strcmp(chCaps,chFind) == 0 )
1786 double relint1, absint1, current_error;
1796 current_error < 2.*fabs(RelError[i]) )
1798 lgDisambiguate =
true;
1801 if( ipDisambiguate[i][1] > 0 )
1803 ipDisambiguate[i][2] = j;
1808 ipDisambiguate[i][1] = j;
1816 PredQuan[i] = relint;
1827 fprintf( ioMONITOR,
" monitor error: lgCheckMonitors could not find line ");
1831 " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1833 RelError[i] = 10000000.;
1841 double relint_cb = 0.,
1845 fprintf( ioMONITOR,
" monitor error: lgCheckMonitors could not find line ");
1849 " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1851 RelError[i] = 10000000.;
1859 PredQuan[i] = absint / absint_cb;
1876 fprintf( ioMONITOR,
" monitor error: lgCheckMonitors could not find line ");
1880 " monitor error: The \"save line labels\" command is a good way to get a list of line labels.\n\n");
1882 RelError[i] = 10000000.;
1889 PredQuan[i] = absint;
1899 double hfrac , hefrac;
1914 " monitor error: lgCheckMonitors could not find h ionization fraction \n");
1934 " monitor error: lgCheckMonitors could not find h ionization fraction \n");
1942 PredQuan[i] = hefrac-hfrac;
1949 PredQuan[i] =
cpu.
i().
lgMPI() ? 1. : 0.;
1957 PredQuan[i] = (double)
nzone;
2007 double sumx=0., sumx2=0., average;
2010 for( n=0; n<
nzone; n++ )
2018 average = sumx/
nzone;
2020 sumx = sqrt( (sumx2-
POW2(sumx)/nzone)/(nzone-1) );
2022 PredQuan[i] = sumx / average;
2030 RelError[i] = PredQuan[i];
2069 (4.*STEFAN_BOLTZ),1,4);
2169 long int nISOCaseB = (long)Param[i][0];
2170 long int nelemCaseB = (long)Param[i][1];
2171 string chElemLabelCaseB =
chIonLbl( nelemCaseB+1, nelemCaseB+1-nISOCaseB );
2188 fprintf(ioMONITOR,
" Species nHi nLo Wl Computed Asserted error\n");
2193 for(
long int ipLo=1+iCase; ipLo<
MIN2(10,nHighestPrinted-1); ++ipLo )
2195 for(
long int ipHi=ipLo+1; ipHi<
MIN2(25,nHighestPrinted); ++ipHi )
2200 if( wl < Param[i][2] || wl > Param[i][3] )
2203 double relint, absint, CBrelint, CBabsint;
2206 if( CBrelint < Param[i][4] )
2210 if(
cdLine( chElemLabelCaseB.c_str(), wl, &relint, &absint,
iLineType[i] ) > 0 )
2213 error = (CBabsint - absint)/
MAX2(CBabsint , absint);
2215 error = (CBabsint - absint);
2216 double RelativeError = fabs(error /
AssertError[i]);
2218 if( RelativeError < 1. )
2220 if( RelativeError < 0.25 )
2222 fprintf( ioMONITOR,
" ChkMonitor ");
2224 else if( RelativeError < 0.50 )
2226 fprintf( ioMONITOR,
" ChkMonitor - ");
2228 else if( RelativeError < 0.75 )
2230 fprintf( ioMONITOR,
" ChkMonitor -- ");
2232 else if( RelativeError < 0.90 )
2234 fprintf( ioMONITOR,
" ChkMonitor --- ");
2236 else if( RelativeError < 0.95 )
2238 fprintf( ioMONITOR,
" ChkMonitor ---- ");
2240 else if( RelativeError < 0.98 )
2242 fprintf( ioMONITOR,
" ChkMonitor ----- ");
2246 fprintf( ioMONITOR,
" ChkMonitor ------ ");
2252 fprintf( ioMONITOR,
" ChkMonitor botch>>");
2254 fprintf(ioMONITOR,
" %s %3li %3li ",
2255 chElemLabelCaseB.c_str() , ipHi , ipLo );
2257 fprintf(ioMONITOR,
" %.2e %.2e %10.3f",
2258 log10(absint) , log10(CBabsint) , error );
2263 fprintf(ioMONITOR ,
" botch \n");
2269 RelError[i] =
MAX2( RelError[i] , fabs(error) );
2283 fprintf(
ioQQQ,
"PROBLEM monitor case B for a He is requested but He is not "
2285 fprintf(
ioQQQ,
"Do not turn off He if you want to monitor its spectrum.\n");
2290 fprintf(ioMONITOR,
" Wl Computed Asserted error\n");
2296 if( wl < Param[i][2] || wl > Param[i][3] )
2298 double relint , absint,CBrelint , CBabsint;
2300 if( CBrelint < Param[i][4] )
2303 if(
cdLine( chElemLabelCaseB.c_str(), wl, &relint, &absint,
iLineType[i] ) > 0)
2306 error = (CBabsint - absint)/
MAX2(CBabsint , absint);
2308 error = (CBabsint - absint);
2309 double RelativeError = fabs(error /
AssertError[i]);
2311 if( RelativeError < 1. )
2313 if( RelativeError < 0.25 )
2315 fprintf( ioMONITOR,
" ChkMonitor ");
2317 else if( RelativeError < 0.50 )
2319 fprintf( ioMONITOR,
" ChkMonitor - ");
2321 else if( RelativeError < 0.75 )
2323 fprintf( ioMONITOR,
" ChkMonitor -- ");
2325 else if( RelativeError < 0.90 )
2327 fprintf( ioMONITOR,
" ChkMonitor --- ");
2329 else if( RelativeError < 0.95 )
2331 fprintf( ioMONITOR,
" ChkMonitor ---- ");
2333 else if( RelativeError < 0.98 )
2335 fprintf( ioMONITOR,
" ChkMonitor ----- ");
2339 fprintf( ioMONITOR,
" ChkMonitor ------ ");
2345 fprintf( ioMONITOR,
" ChkMonitor botch>>");
2348 fprintf(ioMONITOR,
" %.2e %.2e %10.3f",
2349 absint , CBabsint , error );
2354 fprintf(ioMONITOR ,
" botch \n");
2360 RelError[i] =
MAX2( RelError[i] , fabs(error) );
2376 if( this_species.find(
'[' ) == string::npos )
2380 this_species +=
"[2:]";
2382 this_species +=
"[:]";
2385 vector<long> speciesLevels;
2389 fprintf(
ioQQQ,
"PROBLEM Could not find species between quotes: \"%s\".\n",
2396 fprintf(
ioQQQ,
"WARNING Species '%s' has no internal structure."
2397 " Cannot compute departure coefficient\n",
2403 else if( speciesLevels.size() > 0 )
2409 long numPrintLevels = 0;
2410 for( vector<long>::const_iterator ilvl = speciesLevels.begin(); ilvl != speciesLevels.end(); ilvl++ )
2418 ASSERT( numPrintLevels > 0 );
2419 PredQuan[i] /= (double)(numPrintLevels);
2423 if( 0 &&
fp_equal( Param[i][1], 1. ) && PredQuan[i]==0. )
2434 fprintf(
ioQQQ,
"WARNING Requested level(s) in '%s' do not exist\n",
2446 long ipISO = (long)Param[i][0];
2450 fprintf(
ioQQQ,
"PROBLEM asserted element %ld is not turned on!\n",nelem);
2461 PredQuan[i] +=
iso_sp[ipISO][nelem].
st[n].DepartCoef();
2463 ASSERT( numPrintLevels > 0 );
2464 PredQuan[i] /= (double)(numPrintLevels);
2467 if(
fp_equal( Param[i][1], 1. ) && PredQuan[i]==0. )
2491 strcpy( chWeight ,
"VOLUME" );
2496 strcpy( chWeight ,
"RADIUS" );
2512 " monitor error: lgCheckMonitors could not find a line with label %s %f \n",
2522 PredQuan[i] = relint;
2537 if( (relint =
cdH2_colden( (
long)Param[i][0] , (
long)Param[i][1] ) ) < 0. )
2539 fprintf(
ioQQQ,
" PROBLEM lgCheckMonitors did not find v=%li, J=%li for H2 column density.\n",
2540 (
long)Param[i][0] , (
long)Param[i][1] );
2561 " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
2572 PredQuan[i] = relint;
2597 " monitor error: lgCheckMonitors could not find a molecule with label %s %f \n",
2606 PredQuan[i] = relint;
2622 " monitor error: lgCheckMonitors cannot check map since map not done.\n");
2633 " monitor error: lgCheckMonitors cannot check map since temperature not within range.\n");
2675 strcpy( chWeight ,
"VOLUME" );
2680 strcpy( chWeight ,
"RADIUS" );
2690 if( nd >=
gv.
bin.size() )
2693 fprintf(
ioQQQ,
"Use 1 for first grain that is turned on, " );
2695 fprintf(
ioQQQ,
"Old style grain numbers are not valid anymore !!\n" );
2719 " monitor error: lgCheckMonitors could not find an ion with label %s ion %li \n",
2730 PredQuan[i] = relint;
2743 if( nd >=
gv.
bin.size() )
2746 fprintf(
ioQQQ,
"Use 1 for first grain that is turned on, " );
2748 fprintf(
ioQQQ,
"Old style grain numbers are not valid anymore !!\n" );
2761 " monitor error: lgCheckMonitors received an insane chAssertType=%s, impossible\n",
2783 if( RelError[i] <= 0. )
2793 if( RelError[i] >= 0. )
2807 if( lgDisambiguate )
2811 double relint1, relint2, absint1;
2815 fprintf( ioMONITOR,
"=============Line Disambiguation=======================================================\n" );
2816 fprintf( ioMONITOR,
" Wavelengths || Intensities \n" );
2817 fprintf( ioMONITOR,
"Label line match1 match2 match3 || asserted match1 match2 match3\n" );
2821 if( ipDisambiguate[i][1] > 0 )
2830 if( ipDisambiguate[i][2] > 0 )
2837 fprintf( ioMONITOR ,
"--------" );
2846 fprintf( ioMONITOR ,
" %10.3e %10.3e %10.3e %10.3e\n",
2854 fprintf( ioMONITOR ,
" %10.4f %10.4f %10.4f %10.4f\n",
2878 fprintf(ioMONITOR,
"%s", ctime(&now) );
2887 fprintf( ioMONITOR,
" No errors were found. Summary follows.\n");
2891 fprintf( ioMONITOR,
" Errors were found. Summary follows.\n");
2895 " %-*s%*s computed asserted Rel Err Set err type \n",
2911 fprintf( ioMONITOR,
" BIG BOTCHED MONITORS!!! Big Botched Monitors!!! \n");
2915 fprintf( ioMONITOR,
" BOTCHED MONITORS!!! Botched Monitors!!! \n");
2920 fprintf(ioMONITOR,
"\n The mean of the %li monitor Case A, B relative "
2921 "residuals is %.2f\n\n" ,
2944 fprintf( ioMONITOR,
" intr " );
2947 fprintf( ioMONITOR,
" emer " );
2950 fprintf( ioMONITOR,
" intr cumu" );
2953 fprintf( ioMONITOR,
" emer cumu" );
2956 fprintf( ioMONITOR,
"ERROR: Unrecognized line type: %d\n",
2969 double prtPredQuan, prtAssertQuantity;
2982 prtPredQuan = PredQuan;
2991 double relative = fabs( RelError /
SDIV( fabs(AssertError)));
2993 if( relative < 0.25 || chAssertLimit !=
'=' )
2995 fprintf( ioMONITOR,
" ChkMonitor ");
2997 else if( relative < 0.50 )
2999 fprintf( ioMONITOR,
" ChkMonitor - ");
3001 else if( relative < 0.75 )
3003 fprintf( ioMONITOR,
" ChkMonitor -- ");
3005 else if( relative < 0.90 )
3007 fprintf( ioMONITOR,
" ChkMonitor --- ");
3009 else if( relative < 0.95 )
3011 fprintf( ioMONITOR,
" ChkMonitor ---- ");
3013 else if( relative < 0.98 )
3015 fprintf( ioMONITOR,
" ChkMonitor ----- ");
3019 fprintf( ioMONITOR,
" ChkMonitor ------ ");
3024 fprintf( ioMONITOR,
" ChkMonitor botch>>");
3027 fprintf( ioMONITOR ,
"%-*s ",
NCHLAB-1, chAssertLineLabel.c_str() );
3030 if( strcmp( chAssertType,
"Ll" )==0 || strcmp( chAssertType,
"Lr" )==0 ||
3031 strcmp( chAssertType,
"Lb" )==0 )
3033 prt_wl( ioMONITOR , wavelength );
3040 const char* format =
" %10.4f %c %10.4f %7.3f %7.3f ";
3042 format =
" %10.3e %c %10.3e %7.3f %7.3f ";
3059 if( !lg1OK && (fabs(RelError) > 3.*AssertError) && lgFound )
3061 fprintf( ioMONITOR ,
" <<BIG BOTCH!!\n");
bool nMatch(const char *chKey) const
void cdLine_ip(long int ipLine, double *relint, double *absint)
realnum WaveLengthCaseB[8][25][24]
void prt_wl(FILE *ioOUT, realnum wl)
static bool lgSpaceAllocated
void InitMonitorResults(void)
string chIonLbl(const TransitionProxy &t)
realnum WavlenErrorGet(realnum wavelength, long sig_figs)
double ForcePass(char chAssertLimit1)
static vector< bool > lgQuantityLog
map< std::string, std::vector< TransitionProxy > >::iterator blend_iterator
NORETURN void TotalInsanity(void)
double get_error_ratio(double pred, double assert)
t_monitorresults MonitorResults
long int nSumErrorCaseMonitor
void ParseMonitorResults(Parser &p)
int cdIonFrac(const char *chLabel, long int IonStage, double *fracin, const char *chWeight, bool lgDensity)
long int nTotalIoniz_start
int GetQuote(string &chLabel)
static vector< char > chAssertLimit
bool nMatchErase(const char *chKey)
int cdTemp(const char *chLabel, long int IonStage, double *TeMean, const char *chWeight)
double DepartCoef() const
static realnum ErrorDefaultPerformance
static multi_arr< double, 2 > Param
void cap4(char *chCAP, const char *chLab)
static vector< string > strAssertSpecies
void trimTrailingWhiteSpace(string &str)
NORETURN void StringError() const
static t_version & Inst()
t_elementnames elementnames
t_iso_sp iso_sp[NISO][LIMELM]
void prt_line_err(FILE *ioOUT, const char *label, realnum wvlng)
double xIonDense[LIMELM][LIMELM+1]
bool fp_equal(sys_float x, sys_float y, int n=3)
double SumErrorCaseMonitor
double cdH2_colden(long iVib, long iRot)
long int n_HighestResolved_max
static vector< double > AssertError
static vector< double > AssertQuantity
STATIC void prtLineType(FILE *ioMONITOR, const int iLineType)
const molezone * getLevelsGeneric(const char *chLabel, bool lgValidate, vector< long > &LevelList)
NORETURN void NoNumb(const char *chDesc) const
static realnum ErrorDefault
long int GetElem(void) const
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
static vector< std::vector< TransitionProxy > * > assertBlends
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
static vector< string > chAssertLineLabel
enum level_status status() const
void PrtOneMonitor(FILE *ioMONITOR, const char *chAssertType, const string chAssertLineLabel, const realnum wavelength, const int iLineType, const double PredQuan, const char chAssertLimit, const double AssertQuantity, const double RelError, const double AssertError, const bool lg1OK, const bool lgQuantityLog, const bool lgFound)
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
void reserve(size_type i1)
int cdColm(const char *chLabel, long int ion, double *theocl)
bool GetParam(const char *chKey, double *val)
#define DEBUG_ENTRY(funcname)
double powpq(double x, int p, int q)
static vector< realnum > wavelength
static vector< double > AssertQuantity2
bool lgCheckMonitors(FILE *ioMONITOR)
int fprintf(const Output &stream, const char *format,...)
map< std::string, std::vector< TransitionProxy > > blends
sys_float SDIV(sys_float x)
vector< realnum > CaseBWlHeI
int PrintLine(FILE *fp) const
static const int NASSERTS
static const long sig_figs_max
t_secondaries secondaries
static vector< int > iLineType
bool GetRange(const char *chKey, double *val1, double *val2)
static char ** chAssertType
double rate_h2_form_grains_used_total