00001
00002
00003
00004
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 namespace {
00024 class my_Integrand: public std::unary_function<double, double>
00025 {
00026 public:
00027 double operator() (double x)
00028 {
00029 return sin(x);
00030 }
00031
00032 my_Integrand() {}
00033 };
00034
00035
00036
00037 double mySin(double x)
00038 {
00039 return sin(x);
00040 }
00041 }
00042
00043
00044
00045
00046 STATIC void SanityCheckBegin(void );
00047 STATIC void SanityCheckFinal(void );
00048
00049
00050
00051 void SanityCheck( const char *chJob )
00052 {
00053 DEBUG_ENTRY( "SanityCheck()" );
00054
00055 if( strcmp(chJob,"begin") == 0 )
00056 {
00057 SanityCheckBegin();
00058 }
00059 else if( strcmp(chJob,"final") == 0 )
00060 {
00061 SanityCheckFinal();
00062 }
00063 else
00064 {
00065 fprintf(ioQQQ,"SanityCheck called with insane argument.\n");
00066 cdEXIT(EXIT_FAILURE);
00067 }
00068 }
00069
00070 STATIC void SanityCheckFinal(void )
00071 {
00072
00073 }
00074
00075 STATIC void SanityCheckBegin(void )
00076 {
00077 bool lgOK=true;
00078 int lgFlag;
00079 int32 ner, ipiv[3];
00080 long i ,
00081 j ,
00082 nelem ,
00083 ion ,
00084 nshells;
00085 double *A;
00086
00087
00088 double Aul ,
00089 error,
00090 Z4, gaunt;
00091
00092 long n, logu, loggamma2;
00093
00094 const int NDIM = 10;
00095 double x , ans1 , ans2 , xMatrix[NDIM][NDIM] , yVector[NDIM] ,
00096 rcond;
00097 realnum *fvector;
00098 long int *ipvector;
00099
00100 DEBUG_ENTRY( "SanityCheck()" );
00101
00102
00103
00104
00105
00106
00107
00108
00109 lgOK = true;
00110
00111
00112
00113
00114
00115
00116 for( nelem=0; nelem<LIMELM; ++nelem )
00117 {
00118
00119 if( dense.lgElmtOn[nelem] )
00120 {
00121
00122 Aul = H_Einstein_A( 2, 1, 1, 0, nelem+1 );
00123
00124 if( fabs(Aul - iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul() ) /Aul > 0.01 )
00125 {
00126 fprintf(ioQQQ," SanityCheck found insane H-like As.\n");
00127 lgOK = false;
00128 }
00129 }
00130 }
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 for( logu=-4; logu<=4; logu++)
00141 {
00142
00143 for(loggamma2=-4; loggamma2<=4; loggamma2++)
00144 {
00145 double SutherlandGff[9][9]=
00146 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
00147 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
00148 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
00149 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
00150 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
00151 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
00152 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
00153 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
00154 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};
00155
00156 gaunt = cont_gaunt_calc( TE1RYD/pow(10.,(double)loggamma2), 1., pow(10.,(double)(logu-loggamma2)) );
00157 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
00158
00159 if( error>0.11 || ( loggamma2<2 && error>0.05 ) )
00160 {
00161 fprintf(ioQQQ," SanityCheck found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
00162 logu, loggamma2, error );
00163 lgOK = false;
00164 }
00165 }
00166
00167 }
00168
00169
00170
00171
00172
00173
00174 for( nelem=1; nelem<LIMELM; ++nelem )
00175 {
00176
00177 double as[]={
00178
00179 0. ,9.47e+006 ,3.44e+008 ,1.74e+009 ,5.51e+009 ,1.34e+010 ,
00180 2.79e+010 ,5.32E+010 ,8.81e+010 ,1.46E+011 ,2.15e+011 ,3.15e+011 ,
00181 4.46e+011 ,6.39E+011 ,8.26e+011 ,1.09e+012 ,1.41e+012 ,1.86E+012 ,
00182 2.26e+012 ,2.80e+012 ,3.44e+012 ,4.18e+012 ,5.04e+012 ,6.02e+012 ,
00183 7.14e+012 ,8.40e+012 ,9.83e+012 ,1.14e+013 ,1.32e+013 ,1.52e+013
00184 };
00185
00186 if( iso_sp[ipHE_LIKE][nelem].numLevels_max > 8 && dense.lgElmtOn[nelem])
00187 {
00188
00189 if( fabs( as[nelem] - iso_sp[ipHE_LIKE][nelem].trans(9,1).Emis().Aul() ) /as[nelem] > 0.025 )
00190 {
00191 fprintf(ioQQQ,
00192 " SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n",
00193 nelem,
00194 as[nelem] ,
00195 iso_sp[ipHE_LIKE][nelem].trans(9,1).Emis().Aul() );
00196 lgOK = false;
00197 }
00198 }
00199 }
00200
00201
00202 if( !opac.lgCaseB )
00203 {
00204
00205 for( i = 0; i <=110; i++ )
00206 {
00207 double DrakeTotalAuls[111] = {
00208 -1.0000E+00, -1.0000E+00, -1.0000E+00, 1.02160E+07,
00209 1.02160E+07, 1.02160E+07, 1.80090E+09, 2.78530E+07,
00210 1.82990E+07, 1.05480E+07, 7.07210E+07, 6.37210E+07,
00211 5.79960E+08, 1.60330E+07, 1.13640E+07, 7.21900E+06,
00212 3.11920E+07, 2.69830E+07, 1.38380E+07, 1.38330E+07,
00213 2.52270E+08, 9.20720E+06, 6.82220E+06, 4.56010E+06,
00214 1.64120E+07, 1.39290E+07, 7.16030E+06, 7.15560E+06,
00215 4.25840E+06, 4.25830E+06, 1.31150E+08, 5.62960E+06,
00216 4.29430E+06, 2.95570E+06, 9.66980E+06, 8.12340E+06,
00217 4.19010E+06, 4.18650E+06, 2.48120E+06, 2.48120E+06,
00218 1.64590E+06, 1.64590E+06, 7.65750E+07, 3.65330E+06,
00219 2.84420E+06, 1.99470E+06, 6.16640E+06, 5.14950E+06,
00220 2.66460E+06, 2.66200E+06, 1.57560E+06, 1.57560E+06,
00221 1.04170E+06, 1.04170E+06, 7.41210E+05, 7.41210E+05,
00222 4.84990E+07, 2.49130E+06, 1.96890E+06, 1.39900E+06,
00223 4.16900E+06, 3.46850E+06, 1.79980E+06, 1.79790E+06,
00224 1.06410E+06, 1.06410E+06, 7.02480E+05, 7.02480E+05,
00225 4.98460E+05, 4.98460E+05, -1.0000E+00, -1.0000E+00,
00226 3.26190E+07, 1.76920E+06, 1.41440E+06, 1.01460E+06,
00227 2.94830E+06, 2.44680E+06, 1.27280E+06, 1.27140E+06,
00228 7.52800E+05, 7.52790E+05, 4.96740E+05, 4.96740E+05,
00229 3.51970E+05, 3.51970E+05, -1.0000E+00, -1.0000E+00,
00230 -1.0000E+00, -1.0000E+00, 2.29740E+07, 1.29900E+06,
00231 1.04800E+06, 7.57160E+05, 2.16090E+06, 1.79030E+06,
00232 9.33210E+05, 9.32120E+05, 5.52310E+05, 5.52310E+05,
00233 3.64460E+05, 3.64460E+05, 2.58070E+05, 2.58070E+05,
00234 -1.0000E+00, -1.0000E+00, -1.0000E+00, -1.0000E+00,
00235 -1.0000E+00, -1.0000E+00, 1.67840E+07};
00236
00237 if( DrakeTotalAuls[i] > 0. &&
00238 i < iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max)
00239 {
00240 if( fabs( DrakeTotalAuls[i] - (1./iso_sp[ipHE_LIKE][ipHELIUM].st[i].lifetime()) ) /DrakeTotalAuls[i] > 0.001 )
00241 {
00242 fprintf(ioQQQ,
00243 " SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n",
00244 i,
00245 DrakeTotalAuls[i],
00246 (1./iso_sp[ipHE_LIKE][ipHELIUM].st[i].lifetime()) );
00247 lgOK = false;
00248 }
00249 }
00250 }
00251 }
00252
00253
00254
00255
00256
00257
00258 if( dense.lgElmtOn[ipHELIUM] )
00259 {
00260
00261 const int NHE1CS = 20;
00262 double he1cs[NHE1CS] =
00263 {
00264 5.480E-18 , 9.253E-18 , 1.598E-17 , 1.598E-17 , 1.598E-17 , 1.348E-17 ,
00265 8.025E-18 , 1.449E-17 , 2.852E-17 , 1.848E-17 , 1.813E-17 , 2.699E-17 ,
00266 1.077E-17 , 2.038E-17 , 4.159E-17 , 3.670E-17 , 3.575E-17 , 1.900E-17 ,
00267 1.900E-17 , 4.175E-17
00268 };
00269
00270
00271 j = MIN2( NHE1CS+1 , iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max -iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max );
00272 for( n=1; n<j; ++n )
00273 {
00274
00275 i = iso_sp[ipHE_LIKE][ipHELIUM].fb[n].ipOpac;
00276 ASSERT( i>0 );
00277
00278
00279
00280
00281
00282
00283
00284 if( fabs( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] > 0.15 )
00285 {
00286 fprintf(ioQQQ,
00287 " SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n",
00288 n,
00289 he1cs[n-1] ,
00290 opac.OpacStack[i - 1]);
00291 fprintf(ioQQQ,
00292 " n=%li, l=%li, s=%li\n",
00293 iso_sp[ipHE_LIKE][ipHELIUM].st[n].n() ,
00294 iso_sp[ipHE_LIKE][ipHELIUM].st[n].l() ,
00295 iso_sp[ipHE_LIKE][ipHELIUM].st[n].S());
00296 lgOK = false;
00297 }
00298 }
00299 }
00300
00301 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00302 {
00303 long nelem = ipISO;
00304
00305
00306 if( !iso_ctrl.lgNoRecombInterp[ipISO] )
00307 {
00308
00309 error = fabs( iso_recomb_check( ipISO, nelem , 0 , 7500. ) );
00310 if( error > 0.01 )
00311 {
00312 fprintf(ioQQQ,
00313 " SanityCheck found insane1 %s %s recom coef: expected, n=%i error: %.2e \n",
00314 iso_ctrl.chISO[ipISO],
00315 elementnames.chElementSym[nelem],
00316 0,
00317 error );
00318 lgOK = false;
00319 }
00320
00321
00322 error = fabs( iso_recomb_check( ipISO, nelem , 1 , 12500. ) );
00323 if( error > 0.01 )
00324 {
00325 fprintf(ioQQQ,
00326 " SanityCheck found insane2 %s %s recom coef: expected, n=%i error: %.2e \n",
00327 iso_ctrl.chISO[ipISO],
00328 elementnames.chElementSym[nelem],
00329 1,
00330 error );
00331 lgOK = false;
00332 }
00333 }
00334 }
00335
00336
00337
00338
00339
00340
00341
00342 const int NSORT = 100 ;
00343
00344 fvector = (realnum *)MALLOC((NSORT)*sizeof(realnum) );
00345
00346 ipvector = (long *)MALLOC((NSORT)*sizeof(long int) );
00347
00348 nelem = 1;
00349
00350 for( i=0; i<NSORT; ++i )
00351 {
00352 nelem *= -1;
00353 fvector[i] = (realnum)nelem * ((realnum)NSORT-i);
00354 }
00355
00356
00357 spsort(fvector,
00358 NSORT,
00359 ipvector,
00360
00361
00362 1,
00363 &lgFlag);
00364
00365 if( lgFlag ) lgOK = false;
00366
00367 for( i=1; i<NSORT; ++i )
00368 {
00369
00370
00371 if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] )
00372 {
00373 fprintf(ioQQQ," SanityCheck found insane sort\n");
00374 lgOK = false;
00375 }
00376 }
00377
00378 free( fvector );
00379 free( ipvector);
00380
00381 # if 0
00382 ttemp = (realnum)sqrt(phycon.te);
00383
00384 if( fabs(ttemp - phycon.sqrte )/ttemp > 1e-5 )
00385 {
00386 fprintf(ioQQQ , "SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n",
00387 phycon.te ,
00388 sqrt(phycon.te) ,
00389 phycon.sqrte ,
00390 fabs(sqrt(phycon.te) - phycon.sqrte ) );
00391 lgOK = false;
00392 }
00393 # endif
00394
00395
00396
00397
00398
00399
00400
00401
00402 # if 0
00403
00404
00405 # if !defined(NDEBUG)
00406 x = 0.;
00407 for( i=1; i<rfield.nupper-1; ++i )
00408 {
00409 if( fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.) > 0.02 )
00410 {
00411 ans1 = fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.);
00412 fprintf(ioQQQ," SanityCheck found insane widflx anu[i+1]=%e anu[i]=%e widflx=%e delta=%e rel err %e\n",
00413 rfield.anu[i+1] , rfield.anu[i] , rfield.widflx[i] , rfield.anu[i+1] -rfield.anu[i] , ans1 );
00414 lgOK = false;
00415 x = MAX2( ans1 , x);
00416 }
00417
00418
00419 ans1 = rfield.widflx[i] / rfield.anu[i];
00420 if( (rfield.anu[i]+rfield.widflx[i]/2.)*(1.-ans1/10.) > rfield.anu[i+1] - rfield.widflx[i+1]/2.)
00421 {
00422 fprintf(ioQQQ," SanityCheck found insane overlap1 widflx %e %e %e %e %e %e\n",
00423 rfield.anu[i] , rfield.widflx[i], rfield.anu[i] + rfield.widflx[i]/2. , rfield.anu[i+1],
00424 rfield.widflx[i+1], rfield.anu[i+1] -rfield.widflx[i+1]/2. );
00425 lgOK = false;
00426 }
00427 if( !lgOK )
00428 {
00429 fprintf(ioQQQ," big error was %e\n", x);
00430 }
00431 }
00432 # endif
00433 # endif
00434
00435
00436
00437
00438
00439
00440
00441 for( nelem=0; nelem<2; ++nelem )
00442 {
00443
00444 if( dense.lgElmtOn[nelem] )
00445 {
00446
00447
00448 Z4 = (double)(nelem+1);
00449 Z4 *= Z4;
00450 Z4 *= Z4;
00451
00452 ans1 = iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul();
00453 ans2 = 6.265e8*Z4;
00454 if( fabs(ans1-ans2)/ans2 > 1e-3 )
00455 {
00456 fprintf(ioQQQ , "SanityCheck finds insane A for H Lya %g %g nelem=%li\n",
00457 ans1 , ans2 , nelem );
00458 lgOK = false;
00459 }
00460 }
00461 }
00462
00463
00464 for( nelem=0; nelem<LIMELM; ++nelem )
00465 {
00466 if( dense.lgElmtOn[nelem] )
00467 {
00468 int ipHi, ipLo;
00469 for( ipHi=4; ipHi< iso_sp[ipH_LIKE][nelem].numLevels_max-iso_sp[ipH_LIKE][nelem].nCollapsed_max; ++ipHi )
00470 {
00471 double sum = 0.;
00472 for( ipLo=0; ipLo<ipHi; ++ipLo )
00473 {
00474 sum += iso_sp[ipH_LIKE][nelem].BranchRatio[ipHi][ipLo];
00475 }
00476 if( fabs(sum-1.)>0.01 )
00477 {
00478 fprintf(ioQQQ ,
00479 "SanityCheck H branching ratio sum not unity for nelem=%li upper level=%i sum=%.3e\n",
00480 nelem, ipHi, sum );
00481 lgOK = false;
00482 }
00483 }
00484 }
00485 }
00486
00487
00488 for( nelem=0; nelem<LIMELM; ++nelem )
00489 {
00490 if( dense.lgElmtOn[nelem] )
00491 {
00492 int ipHi, ipLo;
00493 for( ipHi=1; ipHi< iso_sp[ipH_LIKE][nelem].numLevels_max-iso_sp[ipH_LIKE][nelem].nCollapsed_max; ++ipHi )
00494 {
00495 double inverse_sum = 0.;
00496 double sum = 0.;
00497 long ipISO = ipH_LIKE;
00498
00499
00500 if( L_(ipHi)==0 )
00501 continue;
00502
00503 for( ipLo=0; ipLo<ipHi; ++ipLo )
00504 {
00505 sum += iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul();
00506 }
00507
00508 inverse_sum = 1./sum;
00509
00510 double lifetime = iso_state_lifetime( ipH_LIKE, nelem, N_(ipHi), L_(ipHi) );
00511 double percent_error = (1. - inverse_sum/lifetime)*100;
00512
00513
00514
00515
00516 if( fabs(percent_error) > 0.5 )
00517 {
00518 fprintf(ioQQQ ,
00519 "SanityCheck hydrogenic lifetime not consistent with closed form for nelem=%2li upper level=%2i inverse-sum= %.3e lifetime= %.3e error= %.2f%%\n",
00520 nelem, ipHi, inverse_sum, lifetime, percent_error );
00521 lgOK = false;
00522 }
00523 }
00524 }
00525 }
00526
00527
00528
00529 long ipISO = ipH_LIKE;
00530 nelem = 0;
00531 for( long n=1; n <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++n )
00532 {
00533
00534
00535
00536
00537
00538
00539
00540
00541 if( n > 60 )
00542 break;
00543
00544 double rel_photon_energy = 1. + FLT_EPSILON*2.;
00545 for( long l=0; l < n; l++ )
00546 {
00547 double cs;
00548 long index = iso_sp[ipISO][nelem].QuantumNumbers2Index[n][l][2];
00549
00550 cs = H_photo_cs( rel_photon_energy, n, l, nelem+1 );
00551
00552 error = fabs(cs - opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1] )/
00553 ( (cs + opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1] )/2.);
00554
00555 if( error > 0.05 )
00556 {
00557 fprintf(ioQQQ , "SanityCheck finds insane H photo cs, n,l = %li, %li, expected, found = %e, %e, error = %e\n",
00558 n, l, opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1], cs, error );
00559 lgOK = false;
00560 }
00561 }
00562 }
00563
00564
00565
00566
00567
00568
00569 ASSERT( fp_equal( qg32( 0., PI, mySin ) , 2. ) );
00570
00571
00572 my_Integrand func;
00573 Integrator<my_Integrand,Gaussian32> mySine;
00574 ASSERT( fp_equal( mySine.sum( 0., PI, func ), 2. ) );
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584 x = 1e-3;
00585 do
00586 {
00587
00588 ans1 = ee1(x);
00589 ans2 = expn( 1 , x );
00590 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
00591 {
00592 fprintf(ioQQQ , "SanityCheck finds insane E1 %g %g %g\n",
00593 x , ans1 , ans2 );
00594 lgOK = false;
00595 }
00596
00597
00598 ans1 = e2(x);
00599 ans2 = expn( 2 , x );
00600 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
00601 {
00602 fprintf(ioQQQ , "SanityCheck finds insane E2 %g %g %g\n",
00603 x , ans1 , ans2 );
00604 lgOK = false;
00605 }
00606
00607
00608 x *= 2.;
00609
00610 } while( x < 64. );
00611
00612
00613
00614
00615
00616
00617
00618
00619 yVector[0] = 1.;
00620 yVector[1] = 3.;
00621 yVector[2] = 3.;
00622
00623
00624 for(i=0;i<3;++i)
00625 {
00626 for( j=0;j<3;++j )
00627 {
00628 xMatrix[i][j] = 0.;
00629 }
00630 }
00631
00632
00633 xMatrix[0][0] = 1.;
00634 xMatrix[0][1] = 1.;
00635 xMatrix[1][1] = 1.;
00636 xMatrix[2][2] = 1.;
00637
00638
00639
00640
00641
00642
00643 A = (double*)MALLOC( sizeof(double)*NDIM*NDIM );
00644
00645
00646 for( i=0; i < 3; ++i )
00647 {
00648 for( j=0; j < 3; ++j )
00649 {
00650 A[i*NDIM+j] = xMatrix[i][j];
00651 }
00652 }
00653
00654 ner = 0;
00655
00656
00657
00658 getrf_wrapper(3, 3, A, NDIM, ipiv, &ner);
00659 if( ner != 0 )
00660 {
00661 fprintf( ioQQQ, " SanityCheck DGETRF error\n" );
00662 cdEXIT(EXIT_FAILURE);
00663 }
00664
00665
00666
00667
00668
00669
00670
00671 getrs_wrapper('N', 3, 1, A, NDIM, ipiv, yVector, 3, &ner);
00672
00673 if( ner != 0 )
00674 {
00675 fprintf( ioQQQ, " SanityCheck DGETRS error\n" );
00676 cdEXIT(EXIT_FAILURE);
00677 }
00678
00679 free( A );
00680
00681
00682
00683
00684
00685
00686 rcond = 0.;
00687 for(i=0;i<3;++i)
00688 {
00689 x = fabs( yVector[i]-i-1.);
00690 rcond = MAX2( rcond, x );
00691
00692 }
00693
00694 if( rcond>DBL_EPSILON)
00695 {
00696 fprintf(ioQQQ,
00697 "SanityCheck found too large a deviation in matrix solver = %g \n",
00698 rcond);
00699
00700 lgOK = false;
00701 }
00702
00703
00704
00705
00706
00707 for( nelem=0; nelem<LIMELM; ++nelem )
00708 {
00709
00710 realnum xIonF[NISO][LIMELM];
00711 double hbn[NISO][LIMELM], hn[NISO][LIMELM];
00712
00713 if( dense.lgElmtOn[nelem] )
00714 {
00715
00716 hbn[ipH_LIKE][nelem] = iso_sp[ipH_LIKE][nelem].fb[0].DepartCoef;
00717 hn[ipH_LIKE][nelem] = iso_sp[ipH_LIKE][nelem].st[0].Pop();
00718 xIonF[ipH_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipH_LIKE];
00719
00720 iso_sp[ipH_LIKE][nelem].fb[0].DepartCoef = 0.;
00721 iso_sp[ipH_LIKE][nelem].st[0].Pop() = 1.;
00722 dense.xIonDense[nelem][nelem+1-ipH_LIKE] = 1.;
00723
00724 if( nelem > ipHYDROGEN )
00725 {
00726
00727 hbn[ipHE_LIKE][nelem] = iso_sp[ipHE_LIKE][nelem].fb[0].DepartCoef;
00728 hn[ipHE_LIKE][nelem] = iso_sp[ipHE_LIKE][nelem].st[0].Pop();
00729 xIonF[ipHE_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipHE_LIKE];
00730
00731
00732 iso_sp[ipHE_LIKE][nelem].fb[0].DepartCoef = 0.;
00733 iso_sp[ipHE_LIKE][nelem].st[0].Pop() = 1.;
00734 dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = 1.;
00735 }
00736
00737 for( ion=0; ion<=nelem; ++ion )
00738 {
00739
00740 for( nshells=0; nshells<Heavy.nsShells[nelem][ion]; ++nshells )
00741 {
00742 for( j=0; j<3; ++j )
00743 {
00744
00745
00746 if( opac.ipElement[nelem][ion][nshells][j] <=0 )
00747 {
00748
00749 fprintf(ioQQQ,
00750 "SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n",
00751 nelem , ion , nshells, j );
00752 fprintf(ioQQQ,
00753 "value was %li \n", opac.ipElement[nelem][ion][nshells][j] );
00754
00755 lgOK = false;
00756 }
00757 }
00758 }
00759
00760 if( nelem > 1 )
00761 {
00762 vector<realnum> saveion;
00763 saveion.resize(nelem+2);
00764
00765 for( j=0; j <= (nelem + 1); j++ )
00766 {
00767 saveion[j] = dense.xIonDense[nelem][j];
00768 dense.xIonDense[nelem][j] = 0.;
00769 }
00770
00771 dense.xIonDense[nelem][ion] = 1.;
00772
00773 OpacityZero();
00774 opac.lgRedoStatic = true;
00775
00776
00777
00778 OpacityAdd1Element(nelem);
00779
00780
00781 for( j=Heavy.ipHeavy[nelem][ion]; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
00782 {
00783 if( opac.opacity_abs[j]+opac.OpacStatic[j] < FLT_MIN )
00784 {
00785
00786 fprintf(ioQQQ,
00787 "SanityCheck found non-positive photo cs for nelem=%li ion=%li \n",
00788 nelem , ion );
00789 fprintf(ioQQQ,
00790 "value was %.2e + %.2e nelem %li ion %li at energy %.2e\n",
00791 opac.opacity_abs[j] ,
00792 opac.OpacStatic[j] ,
00793 nelem ,
00794 ion ,
00795 rfield.anu[j]);
00796
00797 lgOK = false;
00798 break;
00799 }
00800 }
00801
00802 for( j=0; j <= (nelem + 1); j++ )
00803 {
00804 dense.xIonDense[nelem][j] = saveion[j];
00805 }
00806
00807 }
00808 }
00809 iso_sp[ipH_LIKE][nelem].fb[ipH1s].DepartCoef = hbn[ipH_LIKE][nelem];
00810 iso_sp[ipH_LIKE][nelem].st[ipH1s].Pop() = hn[ipH_LIKE][nelem];
00811 dense.xIonDense[nelem][nelem+1-ipH_LIKE] = xIonF[ipH_LIKE][nelem];
00812
00813 if( nelem > ipHYDROGEN )
00814 {
00815 iso_sp[ipHE_LIKE][nelem].fb[ipHe1s1S].DepartCoef = hbn[ipHE_LIKE][nelem];
00816 iso_sp[ipHE_LIKE][nelem].st[ipHe1s1S].Pop() = hn[ipHE_LIKE][nelem];
00817 dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = xIonF[ipHE_LIKE][nelem];
00818 }
00819 }
00820 }
00821
00822
00823
00824
00825
00826
00827
00828
00829 if( lgOK )
00830 {
00831
00832 if( trace.lgTrace )
00833 {
00834 fprintf( ioQQQ, " SanityCheck returns OK\n");
00835 }
00836 return;
00837 }
00838
00839 else
00840 {
00841
00842 fprintf(ioQQQ , "SanityCheck finds insanity so exiting\n");
00843 ShowMe();
00844 cdEXIT(EXIT_FAILURE);
00845 }
00846 }