cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*SanityCheck, check that various parts of the code still work, called by Cloudy after continuum 00004 * and optical depth arrays are set up, but before initial temperature and ionization */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "thirdparty.h" 00008 #include "dense.h" 00009 #include "elementnames.h" 00010 #include "continuum.h" 00011 #include "helike_recom.h" 00012 #include "rfield.h" 00013 #include "taulines.h" 00014 #include "hypho.h" 00015 #include "iso.h" 00016 #include "opacity.h" 00017 #include "hydro_bauman.h" 00018 #include "hydrogenic.h" 00019 #include "heavy.h" 00020 #include "trace.h" 00021 #include "cloudy.h" 00022 00023 /* NB - this routine must not change any global variables - any that are changed as part of 00024 * a test must be reset, so that the code retains state */ 00025 00026 STATIC void SanityCheckBegin(void ); 00027 STATIC void SanityCheckFinal(void ); 00028 /* chJob is either "begin" or "final" 00029 * "begin is before code starts up 00030 * "final is after model is complete */ 00031 void SanityCheck( const char *chJob ) 00032 { 00033 DEBUG_ENTRY( "SanityCheck()" ); 00034 00035 if( strcmp(chJob,"begin") == 0 ) 00036 { 00037 SanityCheckBegin(); 00038 } 00039 else if( strcmp(chJob,"final") == 0 ) 00040 { 00041 SanityCheckFinal(); 00042 } 00043 else 00044 { 00045 fprintf(ioQQQ,"SanityCheck called with insane argument.\n"); 00046 cdEXIT(EXIT_FAILURE); 00047 } 00048 } 00049 00050 STATIC void SanityCheckFinal(void ) 00051 { 00052 /* PrtComment also has some ending checks on sanity */ 00053 } 00054 00055 STATIC void SanityCheckBegin(void ) 00056 { 00057 bool lgOK=true; 00058 int lgFlag;// error return for spsort, 0 success, >=1 for errors 00059 int32 ner, ipiv[3]; 00060 long i , 00061 j , 00062 nelem , 00063 ion , 00064 nshells; 00065 double *A; 00066 00067 /* this will be charge to the 4th power */ 00068 double Aul , 00069 error, 00070 Z4, gaunt; 00071 00072 long n, logu, loggamma2; 00073 00074 const int NDIM = 10; 00075 double x , ans1 , ans2 , xMatrix[NDIM][NDIM] , yVector[NDIM] , 00076 rcond; 00077 realnum *fvector; 00078 long int *ipvector; 00079 00080 DEBUG_ENTRY( "SanityCheck()" ); 00081 00082 /********************************************************* 00083 * * 00084 * confirm that various part of cloudy still work * 00085 * * 00086 *********************************************************/ 00087 00088 /* if this is no longer true at end, we have a problem */ 00089 lgOK = true; 00090 00091 /********************************************************* 00092 * * 00093 * check that all the Lyas As are ok * 00094 * * 00095 *********************************************************/ 00096 for( nelem=0; nelem<LIMELM; ++nelem ) 00097 { 00098 /* this element may be turned off */ 00099 if( dense.lgElmtOn[nelem] ) 00100 { 00101 /* H_Einstein_A( n, l, np, lp, iz ) - all are on physics scale */ 00102 Aul = H_Einstein_A( 2, 1, 1, 0, nelem+1 ); 00103 /*fprintf(ioQQQ,"%li\t%.4e\n", nelem+1, Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul );*/ 00104 if( fabs(Aul - Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul ) /Aul > 0.01 ) 00105 { 00106 fprintf(ioQQQ," SanityCheck found insane H-like As.\n"); 00107 lgOK = false; 00108 } 00109 } 00110 } 00111 00112 /********************************************************* 00113 * * 00114 * check that gaunt factors are good * 00115 * * 00116 *********************************************************/ 00117 /* Uncommenting each of the four print statements here 00118 * will produce a nice table comparable to Sutherland 98, Table 2. */ 00119 /* fprintf(ioQQQ,"u\t-4\t-3\t-2\t-1\t0\t1\t2\t3\t4\n");*/ 00120 for( logu=-4; logu<=4; logu++) 00121 { 00122 /*fprintf(ioQQQ,"%li\t", logu);*/ 00123 for(loggamma2=-4; loggamma2<=4; loggamma2++) 00124 { 00125 double SutherlandGff[9][9]= 00126 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008}, 00127 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041}, 00128 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771}, 00129 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747}, 00130 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237}, 00131 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202}, 00132 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355}, 00133 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680}, 00134 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}}; 00135 00136 gaunt = cont_gaunt_calc( TE1RYD/pow(10.,(double)loggamma2), 1., pow(10.,(double)(logu-loggamma2)) ); 00137 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt; 00138 /*fprintf(ioQQQ,"%1.3f\t", gaunt);*/ 00139 if( error>0.11 || ( loggamma2<2 && error>0.05 ) ) 00140 { 00141 fprintf(ioQQQ," SanityCheck found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n", 00142 logu, loggamma2, error ); 00143 lgOK = false; 00144 } 00145 } 00146 /*fprintf(ioQQQ,"\n");*/ 00147 } 00148 00149 /********************************************************* 00150 * * 00151 * check some transition probabililties for he-like ions * 00152 * * 00153 *********************************************************/ 00154 for( nelem=1; nelem<LIMELM; ++nelem ) 00155 { 00156 /* the helike 9-1 transition, A(3^3P to 2^3S) */ 00157 double as[]={ 00158 /* updated with Johnson values */ 00159 0. ,9.47e+006 ,3.44e+008 ,1.74e+009 ,5.51e+009 ,1.34e+010 , 00160 2.79e+010 ,5.32E+010 ,8.81e+010 ,1.46E+011 ,2.15e+011 ,3.15e+011 , 00161 4.46e+011 ,6.39E+011 ,8.26e+011 ,1.09e+012 ,1.41e+012 ,1.86E+012 , 00162 2.26e+012 ,2.80e+012 ,3.44e+012 ,4.18e+012 ,5.04e+012 ,6.02e+012 , 00163 7.14e+012 ,8.40e+012 ,9.83e+012 ,1.14e+013 ,1.32e+013 ,1.52e+013 00164 }; 00165 00166 if( iso.numLevels_max[ipHE_LIKE][nelem] > 8 && dense.lgElmtOn[nelem]) 00167 { 00168 /* following used to print current values of As */ 00169 /*fprintf(ioQQQ,"%.2e\n", Transitions[ipHE_LIKE][nelem][9][1].Emis->Aul );*/ 00170 if( fabs( as[nelem] - Transitions[ipHE_LIKE][nelem][9][1].Emis->Aul ) /as[nelem] > 0.025 ) 00171 { 00172 fprintf(ioQQQ, 00173 " SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n", 00174 nelem, 00175 as[nelem] , 00176 Transitions[ipHE_LIKE][nelem][9][1].Emis->Aul ); 00177 lgOK = false; 00178 } 00179 } 00180 } 00181 00182 /* only do this test if case b is not in effect */ 00183 if( !opac.lgCaseB ) 00184 { 00185 00186 for( i = 0; i <=110; i++ ) 00187 { 00188 double DrakeTotalAuls[111] = { 00189 -1.0000E+00, -1.0000E+00, -1.0000E+00, 1.02160E+07, 00190 1.02160E+07, 1.02160E+07, 1.80090E+09, 2.78530E+07, 00191 1.82990E+07, 1.05480E+07, 7.07210E+07, 6.37210E+07, 00192 5.79960E+08, 1.60330E+07, 1.13640E+07, 7.21900E+06, 00193 3.11920E+07, 2.69830E+07, 1.38380E+07, 1.38330E+07, 00194 2.52270E+08, 9.20720E+06, 6.82220E+06, 4.56010E+06, 00195 1.64120E+07, 1.39290E+07, 7.16030E+06, 7.15560E+06, 00196 4.25840E+06, 4.25830E+06, 1.31150E+08, 5.62960E+06, 00197 4.29430E+06, 2.95570E+06, 9.66980E+06, 8.12340E+06, 00198 4.19010E+06, 4.18650E+06, 2.48120E+06, 2.48120E+06, 00199 1.64590E+06, 1.64590E+06, 7.65750E+07, 3.65330E+06, 00200 2.84420E+06, 1.99470E+06, 6.16640E+06, 5.14950E+06, 00201 2.66460E+06, 2.66200E+06, 1.57560E+06, 1.57560E+06, 00202 1.04170E+06, 1.04170E+06, 7.41210E+05, 7.41210E+05, 00203 4.84990E+07, 2.49130E+06, 1.96890E+06, 1.39900E+06, 00204 4.16900E+06, 3.46850E+06, 1.79980E+06, 1.79790E+06, 00205 1.06410E+06, 1.06410E+06, 7.02480E+05, 7.02480E+05, 00206 4.98460E+05, 4.98460E+05, -1.0000E+00, -1.0000E+00, 00207 3.26190E+07, 1.76920E+06, 1.41440E+06, 1.01460E+06, 00208 2.94830E+06, 2.44680E+06, 1.27280E+06, 1.27140E+06, 00209 7.52800E+05, 7.52790E+05, 4.96740E+05, 4.96740E+05, 00210 3.51970E+05, 3.51970E+05, -1.0000E+00, -1.0000E+00, 00211 -1.0000E+00, -1.0000E+00, 2.29740E+07, 1.29900E+06, 00212 1.04800E+06, 7.57160E+05, 2.16090E+06, 1.79030E+06, 00213 9.33210E+05, 9.32120E+05, 5.52310E+05, 5.52310E+05, 00214 3.64460E+05, 3.64460E+05, 2.58070E+05, 2.58070E+05, 00215 -1.0000E+00, -1.0000E+00, -1.0000E+00, -1.0000E+00, 00216 -1.0000E+00, -1.0000E+00, 1.67840E+07}; 00217 00218 if( DrakeTotalAuls[i] > 0. && 00219 i < iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM]) 00220 { 00221 if( fabs( DrakeTotalAuls[i] - (1./StatesElem[ipHE_LIKE][ipHELIUM][i].lifetime) ) /DrakeTotalAuls[i] > 0.001 ) 00222 { 00223 fprintf(ioQQQ, 00224 " SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n", 00225 i, 00226 DrakeTotalAuls[i], 00227 (1./StatesElem[ipHE_LIKE][ipHELIUM][i].lifetime) ); 00228 lgOK = false; 00229 } 00230 } 00231 } 00232 } 00233 00234 /********************************************************* 00235 * * 00236 * check the threshold photoionization cs for He I * 00237 * * 00238 *********************************************************/ 00239 if( dense.lgElmtOn[ipHELIUM] ) 00240 { 00241 /* HeI photoionization cross sections at threshold for lowest 20 levels */ 00242 const int NHE1CS = 20; 00243 double he1cs[NHE1CS] = 00244 { 00245 5.480E-18 , 9.253E-18 , 1.598E-17 , 1.598E-17 , 1.598E-17 , 1.348E-17 , 00246 8.025E-18 , 1.449E-17 , 2.852E-17 , 1.848E-17 , 1.813E-17 , 2.699E-17 , 00247 1.077E-17 , 2.038E-17 , 4.159E-17 , 3.670E-17 , 3.575E-17 , 1.900E-17 , 00248 1.900E-17 , 4.175E-17 00249 }; 00250 00251 /* loop over levels and check on photo cross section */ 00252 j = MIN2( NHE1CS+1 , iso.numLevels_max[ipHE_LIKE][ipHELIUM] -iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] ); 00253 for( n=1; n<j; ++n ) 00254 { 00255 /* above list of levels does not include the ground */ 00256 i = iso.ipOpac[ipHE_LIKE][ipHELIUM][n]; 00257 ASSERT( i>0 ); 00258 /*fprintf(ioQQQ,"%li\t%lin", n , i );*/ 00259 /* >>chng 02 apr 10, from 0.01 to 0.02, values stored 00260 * where taken from calc at low contin resolution, when continuum 00261 * resolution changed this changes too */ 00262 /*fprintf(ioQQQ,"%li %.2e\n", n,( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] );*/ 00263 /* >>chng 02 jul 16, limt from 0.02 to 0.04, so that "set resolution 4" will work */ 00264 /* >>chng 04 may 18, levels 10 and 11 are about 12% off - because of energy binning, chng from 0.08 to 0.15 */ 00265 if( fabs( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] > 0.15 ) 00266 { 00267 fprintf(ioQQQ, 00268 " SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n", 00269 n, 00270 he1cs[n-1] , 00271 opac.OpacStack[i - 1]); 00272 fprintf(ioQQQ, 00273 " n=%li, l=%li, s=%li\n", 00274 StatesElem[ipHE_LIKE][ipHELIUM][n].n , 00275 StatesElem[ipHE_LIKE][ipHELIUM][n].l , 00276 StatesElem[ipHE_LIKE][ipHELIUM][n].S); 00277 lgOK = false; 00278 } 00279 } 00280 } 00281 00282 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ ) 00283 { 00284 long nelem = ipISO; 00285 /* Check for agreement between on-the-fly and interpolation calculations 00286 * of recombination, but only if interpolation is turned on. */ 00287 if( !iso.lgNoRecombInterp[ipISO] ) 00288 { 00289 /* check the recombination coefficients for ground state */ 00290 error = fabs( iso_recomb_check( ipISO, nelem , 0 , 7500. ) ); 00291 if( error > 0.01 ) 00292 { 00293 fprintf(ioQQQ, 00294 " SanityCheck found insane1 %s %s recom coef: expected, n=%i error: %.2e \n", 00295 iso.chISO[ipISO], 00296 elementnames.chElementSym[nelem], 00297 0, 00298 error ); 00299 lgOK = false; 00300 } 00301 00302 /* check the recombination coefficients for ground state of the root of each iso sequence */ 00303 error = fabs( iso_recomb_check( ipISO, nelem , 1 , 12500. ) ); 00304 if( error > 0.01 ) 00305 { 00306 fprintf(ioQQQ, 00307 " SanityCheck found insane2 %s %s recom coef: expected, n=%i error: %.2e \n", 00308 iso.chISO[ipISO], 00309 elementnames.chElementSym[nelem], 00310 1, 00311 error ); 00312 lgOK = false; 00313 } 00314 } 00315 } 00316 00317 /********************************************************* 00318 * * 00319 * check out the sorting routine * 00320 * * 00321 *********************************************************/ 00322 00323 const int NSORT = 100 ; 00324 00325 fvector = (realnum *)MALLOC((NSORT)*sizeof(realnum) ); 00326 00327 ipvector = (long *)MALLOC((NSORT)*sizeof(long int) ); 00328 00329 nelem = 1; 00330 /* make up some unsorted values */ 00331 for( i=0; i<NSORT; ++i ) 00332 { 00333 nelem *= -1; 00334 fvector[i] = (realnum)nelem * ((realnum)NSORT-i); 00335 } 00336 00337 /*spsort netlib routine to sort array returning sorted indices */ 00338 spsort(fvector, 00339 NSORT, 00340 ipvector, 00341 /* flag saying what to do - 1 sorts into increasing order, not changing 00342 * the original routine */ 00343 1, 00344 &lgFlag); 00345 00346 if( lgFlag ) lgOK = false; 00347 00348 for( i=1; i<NSORT; ++i ) 00349 { 00350 /*fprintf(ioQQQ," %li %li %.0f\n", 00351 i, ipvector[i],fvector[ipvector[i]] );*/ 00352 if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] ) 00353 { 00354 fprintf(ioQQQ," SanityCheck found insane sort\n"); 00355 lgOK = false; 00356 } 00357 } 00358 00359 free( fvector ); 00360 free( ipvector); 00361 00362 # if 0 00363 ttemp = (realnum)sqrt(phycon.te); 00364 /* check that the temperatures make sense */ 00365 if( fabs(ttemp - phycon.sqrte )/ttemp > 1e-5 ) 00366 { 00367 fprintf(ioQQQ , "SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n", 00368 phycon.te , 00369 sqrt(phycon.te) , 00370 phycon.sqrte , 00371 fabs(sqrt(phycon.te) - phycon.sqrte ) ); 00372 lgOK = false; 00373 } 00374 # endif 00375 00376 /********************************************************* 00377 * * 00378 * confirm that widflx and anu arrays correspond * 00379 * to one another * 00380 * * 00381 *********************************************************/ 00382 00383 # if 0 00384 /* this check on widflx can't be used since some sharpling curved continua, like laser, 00385 * totally fail due to non-linear nature of widflx and anu relationship */ 00386 # if !defined(NDEBUG) 00387 x = 0.; 00388 for( i=1; i<rfield.nupper-1; ++i ) 00389 { 00390 if( fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.) > 0.02 ) 00391 { 00392 ans1 = fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.); 00393 fprintf(ioQQQ," SanityCheck found insane widflx anu[i+1]=%e anu[i]=%e widflx=%e delta=%e rel err %e\n", 00394 rfield.anu[i+1] , rfield.anu[i] , rfield.widflx[i] , rfield.anu[i+1] -rfield.anu[i] , ans1 ); 00395 lgOK = false; 00396 x = MAX2( ans1 , x); 00397 } 00398 /* problems when at energy where resolution of grid changes dramatically */ 00399 /* this is resolution at current energy */ 00400 ans1 = rfield.widflx[i] / rfield.anu[i]; 00401 if( (rfield.anu[i]+rfield.widflx[i]/2.)*(1.-ans1/10.) > rfield.anu[i+1] - rfield.widflx[i+1]/2.) 00402 { 00403 fprintf(ioQQQ," SanityCheck found insane overlap1 widflx %e %e %e %e %e %e\n", 00404 rfield.anu[i] , rfield.widflx[i], rfield.anu[i] + rfield.widflx[i]/2. , rfield.anu[i+1], 00405 rfield.widflx[i+1], rfield.anu[i+1] -rfield.widflx[i+1]/2. ); 00406 lgOK = false; 00407 } 00408 if( !lgOK ) 00409 { 00410 fprintf(ioQQQ," big error was %e\n", x); 00411 } 00412 } 00413 # endif 00414 # endif 00415 00416 00417 /********************************************************* 00418 * * 00419 * confirm that hydrogen einstein As are still valid * 00420 * * 00421 *********************************************************/ 00422 for( nelem=0; nelem<2; ++nelem ) 00423 { 00424 /* this element may be turned off */ 00425 if( dense.lgElmtOn[nelem] ) 00426 { 00427 /*Z4 = (double)(POW2(nelem+1)*POW2(nelem+1));*/ 00428 /* form charge to the 4th power */ 00429 Z4 = (double)(nelem+1); 00430 Z4 *= Z4; 00431 Z4 *= Z4; 00432 /* H Lya */ 00433 ans1 = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul; 00434 ans2 = 6.265e8*Z4; 00435 if( fabs(ans1-ans2)/ans2 > 1e-3 ) 00436 { 00437 fprintf(ioQQQ , "SanityCheck finds insane A for H Lya %g %g nelem=%li\n", 00438 ans1 , ans2 , nelem ); 00439 lgOK = false; 00440 } 00441 00442 # if 0 00443 /*must disable since, at this time, induced is included in Aul */ 00444 /* H two photon */ 00445 ans1 = Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Aul; 00446 ans2 = 8.226*powi( (1.+nelem) , 6 ); 00447 if( fabs(ans1-ans2)/ans2 > 1e-3 ) 00448 { 00449 fprintf(ioQQQ , "SanityCheck finds insane A for H 2-phot %g %g nelem=%li\n", 00450 ans1 , ans2 , nelem ); 00451 lgOK = false; 00452 } 00453 00454 if( iso.numLevels_max[ipH_LIKE][ipHYDROGEN] > ipH5g ) 00455 { 00456 /* Balmer gamma 5 - 2 + 5-4 + 5-3*/ 00457 ans1 = Transitions[ipH_LIKE][nelem][ipH5p][ipH2s].Emis->Aul + 00458 Transitions[ipH_LIKE][nelem][ipH5s][ipH2p].Emis->Aul + 00459 Transitions[ipH_LIKE][nelem][ipH5d][ipH2p].Emis->Aul + 00460 Transitions[ipH_LIKE][nelem][ipH5p][ipH3s].Emis->Aul + 00461 Transitions[ipH_LIKE][nelem][ipH5s][ipH3p].Emis->Aul + 00462 Transitions[ipH_LIKE][nelem][ipH5d][ipH3p].Emis->Aul + 00463 Transitions[ipH_LIKE][nelem][ipH5p][ipH3d].Emis->Aul + 00464 Transitions[ipH_LIKE][nelem][ipH5f][ipH3d].Emis->Aul + 00465 Transitions[ipH_LIKE][nelem][ipH5p][ipH4s].Emis->Aul + 00466 Transitions[ipH_LIKE][nelem][ipH5s][ipH4p].Emis->Aul + 00467 Transitions[ipH_LIKE][nelem][ipH5d][ipH4p].Emis->Aul + 00468 Transitions[ipH_LIKE][nelem][ipH5p][ipH4d].Emis->Aul + 00469 Transitions[ipH_LIKE][nelem][ipH5f][ipH4d].Emis->Aul + 00470 Transitions[ipH_LIKE][nelem][ipH5d][ipH4f].Emis->Aul + 00471 Transitions[ipH_LIKE][nelem][ipH5g][ipH4f].Emis->Aul; 00472 ans2 = (2.532e6+2.20e6+2.70e6)*Z4; 00473 if( fabs(ans1-ans2)/ans2 > 1e-2 ) 00474 { 00475 fprintf(ioQQQ , 00476 "SanityCheck finds insane A for H5-2 found=%g correct=%g nelem=%li\n", 00477 ans1 , ans2 , nelem ); 00478 lgOK = false; 00479 } 00480 } 00481 # endif 00482 } 00483 } 00484 00485 /* check that hydrogenic branching ratios add up to unity */ 00486 for( nelem=0; nelem<LIMELM; ++nelem ) 00487 { 00488 if( dense.lgElmtOn[nelem] ) 00489 { 00490 int ipHi, ipLo; 00491 for( ipHi=4; ipHi< iso.numLevels_max[ipH_LIKE][nelem]-iso.nCollapsed_max[ipH_LIKE][nelem]; ++ipHi ) 00492 { 00493 double sum = 0.; 00494 for( ipLo=0; ipLo<ipHi; ++ipLo ) 00495 { 00496 sum += iso.BranchRatio[ipH_LIKE][nelem][ipHi][ipLo]; 00497 } 00498 if( fabs(sum-1.)>0.01 ) 00499 { 00500 fprintf(ioQQQ , 00501 "SanityCheck H branching ratio sum not unity for nelem=%li upper n=%i sum=%.3e\n", 00502 nelem, ipHi, sum ); 00503 lgOK = false; 00504 } 00505 } 00506 } 00507 } 00508 00509 fixit();/* make this work */ 00510 #if 0 00511 /* check photo cross sections for H */ 00512 long ipISO = ipH_LIKE; 00513 nelem = 0; 00514 /* loop starts at 3, the first level with n = n and full l */ 00515 for( n=3; n < MIN2(100, iso.n_HighestResolved_max[ipISO][nelem]); ++n ) 00516 { 00517 realnum anu[1]={1.} , cs[1]; 00518 double energy; 00519 long index = iso.QuantumNumbers2Index[ipISO][nelem][n][0][2]; 00520 00521 /* photon energy where cross section will be evaluated, 00522 * this is in Ryd */ 00523 energy = iso.xIsoLevNIonRyd[ipISO][nelem][index]; 00524 anu[0] = (realnum)(energy*1.01); 00525 00526 H_photo_cs( photon , N_(n), L_(n), nelem+1 ); 00527 hypho( 00528 /* Z-Nelec, the residual charge, 1 for hydrogen, 2 for helium */ 00529 (double)(nelem+1), 00530 /* principal quantum number */ 00531 n, 00532 /* lowest angular momentum */ 00533 0, 00534 /* highest angular momentum */ 00535 n-1, 00536 /* scaled lowest photon energy, 00537 * => incorrect?? in units of zed^2/n^2, 00538 * at which the cs will be done */ 00539 energy, 00540 /* number of points to do */ 00541 1, 00542 /* array of energies (in units given above) */ 00543 anu , 00544 /* calculated photoionization cross section in cm^-2 */ 00545 cs); 00546 00547 error = fabs(cs[0] - opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/ 00548 ( (cs[0] + opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/2.); 00549 /*fprintf(ioQQQ,"z=%ld n=%ld error %g old %e new %e\n",nelem,n, error, 00550 opac.OpacStack[iso.ipOpac[ipISO][nelem][n]-1] ,cs[0] );*/ 00551 if( error > 0.05 ) 00552 { 00553 fprintf(ioQQQ , "SanityCheck finds insane H photo cs\n"); 00554 lgOK = false; 00555 } 00556 } 00557 #endif 00558 00559 /********************************************************* 00560 * * 00561 * confirm that exponential integral routines still work * 00562 * * 00563 *********************************************************/ 00564 00565 /* check that first and second exponential integrals are ok, 00566 * step through range of values, beginning with following */ 00567 x = 1e-3; 00568 do 00569 { 00570 /* check that fast e1 routine is ok */ 00571 ans1 = ee1(x); 00572 ans2 = expn( 1 , x ); 00573 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 ) 00574 { 00575 fprintf(ioQQQ , "SanityCheck finds insane E1 %g %g %g\n", 00576 x , ans1 , ans2 ); 00577 lgOK = false; 00578 } 00579 00580 /* check that e2 is ok */ 00581 ans1 = e2(x); 00582 ans2 = expn( 2 , x ); 00583 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 ) 00584 { 00585 fprintf(ioQQQ , "SanityCheck finds insane E2 %g %g %g\n", 00586 x , ans1 , ans2 ); 00587 lgOK = false; 00588 } 00589 00590 /* now increment x */ 00591 x *= 2.; 00592 /* following limit set by sexp returning zero, used in ee1 */ 00593 } while( x < 64. ); 00594 00595 /********************************************************* 00596 * * 00597 * confirm that matrix inversion routine still works * 00598 * * 00599 *********************************************************/ 00600 00601 /* these are the answer, chosen to get xvec 1,2,3 */ 00602 yVector[0] = 1.; 00603 yVector[1] = 3.; 00604 yVector[2] = 3.; 00605 00606 /* zero out the main matrix */ 00607 for(i=0;i<3;++i) 00608 { 00609 for( j=0;j<3;++j ) 00610 { 00611 xMatrix[i][j] = 0.; 00612 } 00613 } 00614 00615 /* remember that order is column, row, alphabetical order, rc */ 00616 xMatrix[0][0] = 1.; 00617 xMatrix[0][1] = 1.; 00618 xMatrix[1][1] = 1.; 00619 xMatrix[2][2] = 1.; 00620 00621 /* this is the default matrix solver */ 00622 /* this test is the 1-d matrix with 2-d macro simulation */ 00623 /* LDA is right dimension of matrix */ 00624 00625 /* MALLOC space for the 1-d array */ 00626 A = (double*)MALLOC( sizeof(double)*NDIM*NDIM ); 00627 00628 /* copy over the main matrix */ 00629 for( i=0; i < 3; ++i ) 00630 { 00631 for( j=0; j < 3; ++j ) 00632 { 00633 A[i*NDIM+j] = xMatrix[i][j]; 00634 } 00635 } 00636 00637 ner = 0; 00638 00639 /*void DGETRF(long,long,double*,long,long[],long*);*/ 00640 /*void DGETRF(int,int,double*,int,int[],int*);*/ 00641 getrf_wrapper(3, 3, A, NDIM, ipiv, &ner); 00642 if( ner != 0 ) 00643 { 00644 fprintf( ioQQQ, " SanityCheck DGETRF error\n" ); 00645 cdEXIT(EXIT_FAILURE); 00646 } 00647 00648 /* usage DGETRS, 'N' = no transpose 00649 * order of matrix, 00650 * number of cols in bvec, =1 00651 * array 00652 * leading dim of array */ 00653 /*void DGETRS(char,int,int,double*,int,int[],double*,int,int*);*/ 00654 getrs_wrapper('N', 3, 1, A, NDIM, ipiv, yVector, 3, &ner); 00655 00656 if( ner != 0 ) 00657 { 00658 fprintf( ioQQQ, " SanityCheck DGETRS error\n" ); 00659 cdEXIT(EXIT_FAILURE); 00660 } 00661 /* release the vector */ 00662 free( A ); 00663 00664 /* now check on validity of the solution, demand that this 00665 * simple problem have come within a few epsilons of the 00666 * correct answer */ 00667 00668 /* find largest deviation */ 00669 rcond = 0.; 00670 for(i=0;i<3;++i) 00671 { 00672 x = fabs( yVector[i]-i-1.); 00673 rcond = MAX2( rcond, x ); 00674 /*printf(" %g ", yVector[i]);*/ 00675 } 00676 00677 if( rcond>DBL_EPSILON) 00678 { 00679 fprintf(ioQQQ, 00680 "SanityCheck found too large a deviation in matrix solver = %g \n", 00681 rcond); 00682 /* set flag saying that things are not ok */ 00683 lgOK = false; 00684 } 00685 /* end matrix inversion check */ 00686 00687 00688 /* these pointers were set to INT_MIN in ContCreatePointers, 00689 * then set to valid numbers in ipShells and OpacityCreate1Element 00690 * this checks that all values have been properly filled */ 00691 for( nelem=0; nelem<LIMELM; ++nelem ) 00692 { 00693 /* must reset state of code after tests performed, remember state here */ 00694 realnum xIonF[NISO][LIMELM]; 00695 double hbn[NISO][LIMELM], hn[NISO][LIMELM]; 00696 00697 if( dense.lgElmtOn[nelem] ) 00698 { 00699 /* set these abundances so that opacities can be checked below */ 00700 hbn[ipH_LIKE][nelem] = iso.DepartCoef[ipH_LIKE][nelem][0]; 00701 hn[ipH_LIKE][nelem] = StatesElem[ipH_LIKE][nelem][0].Pop; 00702 xIonF[ipH_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipH_LIKE]; 00703 00704 iso.DepartCoef[ipH_LIKE][nelem][0] = 0.; 00705 StatesElem[ipH_LIKE][nelem][0].Pop = 1.; 00706 dense.xIonDense[nelem][nelem+1-ipH_LIKE] = 1.; 00707 00708 if( nelem > ipHYDROGEN ) 00709 { 00710 00711 hbn[ipHE_LIKE][nelem] = iso.DepartCoef[ipHE_LIKE][nelem][0]; 00712 hn[ipHE_LIKE][nelem] = StatesElem[ipHE_LIKE][nelem][0].Pop; 00713 xIonF[ipHE_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipHE_LIKE]; 00714 00715 /* this does not exist for hydrogen itself */ 00716 iso.DepartCoef[ipHE_LIKE][nelem][0] = 0.; 00717 StatesElem[ipHE_LIKE][nelem][0].Pop = 1.; 00718 dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = 1.; 00719 } 00720 00721 for( ion=0; ion<=nelem; ++ion ) 00722 { 00723 /* loop over all shells that are defined */ 00724 for( nshells=0; nshells<Heavy.nsShells[nelem][ion]; ++nshells ) 00725 { 00726 for( j=0; j<3; ++j ) 00727 { 00728 /* >>chng 00 apr 05, array index is on fortran scale so must be 00729 * >= 1. This test had been <0, correct for C. Caught by Peter van Hoof */ 00730 if( opac.ipElement[nelem][ion][nshells][j] <=0 ) 00731 { 00732 /* this is not possible */ 00733 fprintf(ioQQQ, 00734 "SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n", 00735 nelem , ion , nshells, j ); 00736 fprintf(ioQQQ, 00737 "value was %li \n", opac.ipElement[nelem][ion][nshells][j] ); 00738 /* set flag saying that things are not ok */ 00739 lgOK = false; 00740 } 00741 } 00742 } 00743 00744 if( nelem > 1 ) 00745 { 00746 realnum saveion[LIMELM+3]; 00747 /* check that photoionization cross sections are ok */ 00748 for( j=1; j <= (nelem + 2); j++ ) 00749 { 00750 saveion[j] = dense.xIonDense[nelem][j-1]; 00751 dense.xIonDense[nelem][j-1] = 0.; 00752 } 00753 00754 dense.xIonDense[nelem][ion] = 1.; 00755 00756 OpacityZero(); 00757 opac.lgRedoStatic = true; 00758 00759 /* generate opacity with standard routine - this is the one 00760 * called in OpacityAddTotal to make opacities in usual calculations */ 00761 OpacityAdd1Element(nelem); 00762 00763 /* this starts one beyond energy of threshold since cs may be zero there */ 00764 for( j=Heavy.ipHeavy[nelem][ion]; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ ) 00765 { 00766 if( opac.opacity_abs[j]+opac.OpacStatic[j] < FLT_MIN ) 00767 { 00768 /* this is not possible */ 00769 fprintf(ioQQQ, 00770 "SanityCheck found non-positive photo cs for nelem=%li ion=%li \n", 00771 nelem , ion ); 00772 fprintf(ioQQQ, 00773 "value was %.2e + %.2e nelem %li ion %li at energy %.2e\n", 00774 opac.opacity_abs[j] , 00775 opac.OpacStatic[j] , 00776 nelem , 00777 ion , 00778 rfield.anu[j]); 00779 /* set flag saying that things are not ok */ 00780 lgOK = false; 00781 break; 00782 } 00783 } 00784 /* reset the ionization distribution */ 00785 for( j=1; j <= (nelem + 2); j++ ) 00786 { 00787 dense.xIonDense[nelem][j-1] = saveion[j]; 00788 } 00789 00790 } 00791 } 00792 iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = hbn[ipH_LIKE][nelem]; 00793 StatesElem[ipH_LIKE][nelem][ipH1s].Pop = hn[ipH_LIKE][nelem]; 00794 dense.xIonDense[nelem][nelem+1-ipH_LIKE] = xIonF[ipH_LIKE][nelem]; 00795 00796 if( nelem > ipHYDROGEN ) 00797 { 00798 iso.DepartCoef[ipHE_LIKE][nelem][ipHe1s1S] = hbn[ipHE_LIKE][nelem]; 00799 StatesElem[ipHE_LIKE][nelem][ipHe1s1S].Pop = hn[ipHE_LIKE][nelem]; 00800 dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = xIonF[ipHE_LIKE][nelem]; 00801 } 00802 } 00803 } 00804 00805 00806 /********************************************************* 00807 * * 00808 * everything is done, all checks make, did we pass them?* 00809 * * 00810 *********************************************************/ 00811 00812 if( lgOK ) 00813 { 00814 /*return if ok */ 00815 if( trace.lgTrace ) 00816 { 00817 fprintf( ioQQQ, " SanityCheck returns OK\n"); 00818 } 00819 return; 00820 } 00821 00822 else 00823 { 00824 /* stop since problem encountered, lgEOF set false */ 00825 fprintf(ioQQQ , "SanityCheck finds insanity so exiting\n"); 00826 ShowMe(); 00827 cdEXIT(EXIT_FAILURE); 00828 } 00829 }