92 long n, logu, loggamma2;
116 for( nelem=0; nelem<
LIMELM; ++nelem )
119 if(
dense.lgElmtOn[nelem] )
126 fprintf(
ioQQQ,
" SanityCheck found insane H-like As.\n");
140 for( logu=-4; logu<=4; logu++)
143 for(loggamma2=-4; loggamma2<=4; loggamma2++)
145 double SutherlandGff[9][9]=
146 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
147 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
148 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
149 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
150 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
151 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
152 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
153 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
154 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};
157 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
159 if( error>0.11 || ( loggamma2<2 && error>0.05 ) )
161 fprintf(
ioQQQ,
" SanityCheck found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
162 logu, loggamma2, error );
174 for( nelem=1; nelem<
LIMELM; ++nelem )
179 0. ,9.47e+006 ,3.44e+008 ,1.74e+009 ,5.51e+009 ,1.34e+010 ,
180 2.79e+010 ,5.32E+010 ,8.81e+010 ,1.46E+011 ,2.15e+011 ,3.15e+011 ,
181 4.46e+011 ,6.39E+011 ,8.26e+011 ,1.09e+012 ,1.41e+012 ,1.86E+012 ,
182 2.26e+012 ,2.80e+012 ,3.44e+012 ,4.18e+012 ,5.04e+012 ,6.02e+012 ,
183 7.14e+012 ,8.40e+012 ,9.83e+012 ,1.14e+013 ,1.32e+013 ,1.52e+013
189 if( fabs( as[nelem] -
iso_sp[
ipHE_LIKE][nelem].trans(9,1).Emis().Aul() ) /as[nelem] > 0.025 )
192 " SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n",
205 for( i = 0; i <=110; i++ )
207 double DrakeTotalAuls[111] = {
208 -1.0000E+00, -1.0000E+00, -1.0000E+00, 1.02160E+07,
209 1.02160E+07, 1.02160E+07, 1.80090E+09, 2.78530E+07,
210 1.82990E+07, 1.05480E+07, 7.07210E+07, 6.37210E+07,
211 5.79960E+08, 1.60330E+07, 1.13640E+07, 7.21900E+06,
212 3.11920E+07, 2.69830E+07, 1.38380E+07, 1.38330E+07,
213 2.52270E+08, 9.20720E+06, 6.82220E+06, 4.56010E+06,
214 1.64120E+07, 1.39290E+07, 7.16030E+06, 7.15560E+06,
215 4.25840E+06, 4.25830E+06, 1.31150E+08, 5.62960E+06,
216 4.29430E+06, 2.95570E+06, 9.66980E+06, 8.12340E+06,
217 4.19010E+06, 4.18650E+06, 2.48120E+06, 2.48120E+06,
218 1.64590E+06, 1.64590E+06, 7.65750E+07, 3.65330E+06,
219 2.84420E+06, 1.99470E+06, 6.16640E+06, 5.14950E+06,
220 2.66460E+06, 2.66200E+06, 1.57560E+06, 1.57560E+06,
221 1.04170E+06, 1.04170E+06, 7.41210E+05, 7.41210E+05,
222 4.84990E+07, 2.49130E+06, 1.96890E+06, 1.39900E+06,
223 4.16900E+06, 3.46850E+06, 1.79980E+06, 1.79790E+06,
224 1.06410E+06, 1.06410E+06, 7.02480E+05, 7.02480E+05,
225 4.98460E+05, 4.98460E+05, -1.0000E+00, -1.0000E+00,
226 3.26190E+07, 1.76920E+06, 1.41440E+06, 1.01460E+06,
227 2.94830E+06, 2.44680E+06, 1.27280E+06, 1.27140E+06,
228 7.52800E+05, 7.52790E+05, 4.96740E+05, 4.96740E+05,
229 3.51970E+05, 3.51970E+05, -1.0000E+00, -1.0000E+00,
230 -1.0000E+00, -1.0000E+00, 2.29740E+07, 1.29900E+06,
231 1.04800E+06, 7.57160E+05, 2.16090E+06, 1.79030E+06,
232 9.33210E+05, 9.32120E+05, 5.52310E+05, 5.52310E+05,
233 3.64460E+05, 3.64460E+05, 2.58070E+05, 2.58070E+05,
234 -1.0000E+00, -1.0000E+00, -1.0000E+00, -1.0000E+00,
235 -1.0000E+00, -1.0000E+00, 1.67840E+07};
237 if( DrakeTotalAuls[i] > 0. &&
243 " SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n",
261 const int NHE1CS = 20;
262 double he1cs[NHE1CS] =
264 5.480E-18 , 9.253E-18 , 1.598E-17 , 1.598E-17 , 1.598E-17 , 1.348E-17 ,
265 8.025E-18 , 1.449E-17 , 2.852E-17 , 1.848E-17 , 1.813E-17 , 2.699E-17 ,
266 1.077E-17 , 2.038E-17 , 4.159E-17 , 3.670E-17 , 3.575E-17 , 1.900E-17 ,
267 1.900E-17 , 4.175E-17
284 if( fabs( he1cs[n-1] -
opac.OpacStack[i - 1] ) /he1cs[n-1] > 0.15 )
287 " SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n",
290 opac.OpacStack[i - 1]);
292 " n=%li, l=%li, s=%li\n",
306 if( !
iso_ctrl.lgNoRecombInterp[ipISO] )
313 " SanityCheck found insane1 %s %s recom coef: expected, n=%i error: %.2e \n",
326 " SanityCheck found insane2 %s %s recom coef: expected, n=%i error: %.2e \n",
342 const int NSORT = 100 ;
346 ipvector = (
long *)
MALLOC((NSORT)*
sizeof(
long int) );
350 for( i=0; i<NSORT; ++i )
365 if( lgFlag ) lgOK =
false;
367 for( i=1; i<NSORT; ++i )
371 if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] )
373 fprintf(
ioQQQ,
" SanityCheck found insane sort\n");
384 if( fabs(ttemp -
phycon.sqrte )/ttemp > 1e-5 )
386 fprintf(
ioQQQ ,
"SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n",
407 for( i=1; i<
rfield.nupper-1; ++i )
412 fprintf(
ioQQQ,
" SanityCheck found insane widflx anu[i+1]=%e anu[i]=%e widflx=%e delta=%e rel err %e\n",
422 fprintf(
ioQQQ,
" SanityCheck found insane overlap1 widflx %e %e %e %e %e %e\n",
429 fprintf(
ioQQQ,
" big error was %e\n", x);
441 for( nelem=0; nelem<2; ++nelem )
444 if(
dense.lgElmtOn[nelem] )
448 Z4 = (double)(nelem+1);
454 if( fabs(ans1-ans2)/ans2 > 1e-3 )
456 fprintf(
ioQQQ ,
"SanityCheck finds insane A for H Lya %g %g nelem=%li\n",
457 ans1 , ans2 , nelem );
464 for( nelem=0; nelem<
LIMELM; ++nelem )
466 if(
dense.lgElmtOn[nelem] )
472 for( ipLo=0; ipLo<ipHi; ++ipLo )
476 if( fabs(sum-1.)>0.01 )
479 "SanityCheck H branching ratio sum not unity for nelem=%li upper level=%i sum=%.3e\n",
488 for( nelem=0; nelem<
LIMELM; ++nelem )
490 if(
dense.lgElmtOn[nelem] )
495 double inverse_sum = 0.;
503 for( ipLo=0; ipLo<ipHi; ++ipLo )
508 inverse_sum = 1./sum;
511 double percent_error = (1. - inverse_sum/lifetime)*100;
516 if( fabs(percent_error) > 0.5 )
519 "SanityCheck hydrogenic lifetime not consistent with closed form for nelem=%2li upper level=%2i inverse-sum= %.3e lifetime= %.3e error= %.2f%%\n",
520 nelem, ipHi, inverse_sum, lifetime, percent_error );
531 for(
long n=1; n <=
iso_sp[ipISO][nelem].n_HighestResolved_max; ++n )
544 double rel_photon_energy = 1. + FLT_EPSILON*2.;
545 for(
long l=0; l < n; l++ )
548 long index =
iso_sp[ipISO][nelem].QuantumNumbers2Index[n][l][2];
550 cs =
H_photo_cs( rel_photon_energy, n, l, nelem+1 );
552 error = fabs(cs -
opac.OpacStack[
iso_sp[ipISO][nelem].fb[index].ipOpac-1] )/
553 ( (cs +
opac.OpacStack[
iso_sp[ipISO][nelem].fb[index].ipOpac-1] )/2.);
557 fprintf(
ioQQQ ,
"SanityCheck finds insane H photo cs, n,l = %li, %li, expected, found = %e, %e, error = %e\n",
558 n, l,
opac.OpacStack[
iso_sp[ipISO][nelem].fb[index].ipOpac-1], cs, error );
589 ans2 =
expn( 1 , x );
590 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
592 fprintf(
ioQQQ ,
"SanityCheck finds insane E1 %g %g %g\n",
599 ans2 =
expn( 2 , x );
600 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
602 fprintf(
ioQQQ ,
"SanityCheck finds insane E2 %g %g %g\n",
643 A = (
double*)
MALLOC(
sizeof(
double)*NDIM*NDIM );
646 for( i=0; i < 3; ++i )
648 for( j=0; j < 3; ++j )
661 fprintf(
ioQQQ,
" SanityCheck DGETRF error\n" );
675 fprintf(
ioQQQ,
" SanityCheck DGETRS error\n" );
690 rcond =
MAX2( rcond, x );
694 if( rcond>DBL_EPSILON)
697 "SanityCheck found too large a deviation in matrix solver = %g \n",
707 for( nelem=0; nelem<
LIMELM; ++nelem )
713 if(
dense.lgElmtOn[nelem] )
737 for( ion=0; ion<=nelem; ++ion )
740 for( nshells=0; nshells<
Heavy.nsShells[nelem][ion]; ++nshells )
746 if(
opac.ipElement[nelem][ion][nshells][j] <=0 )
750 "SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n",
751 nelem , ion , nshells, j );
753 "value was %li \n",
opac.ipElement[nelem][ion][nshells][j] );
762 vector<realnum> saveion;
763 saveion.resize(nelem+2);
765 for( j=0; j <= (nelem + 1); j++ )
767 saveion[j] =
dense.xIonDense[nelem][j];
768 dense.xIonDense[nelem][j] = 0.;
771 dense.xIonDense[nelem][ion] = 1.;
774 opac.lgRedoStatic =
true;
783 if(
opac.opacity_abs[j]+
opac.OpacStatic[j] < FLT_MIN )
787 "SanityCheck found non-positive photo cs for nelem=%li ion=%li \n",
790 "value was %.2e + %.2e nelem %li ion %li at energy %.2e\n",
791 opac.opacity_abs[j] ,
802 for( j=0; j <= (nelem + 1); j++ )
804 dense.xIonDense[nelem][j] = saveion[j];
834 fprintf(
ioQQQ,
" SanityCheck returns OK\n");
842 fprintf(
ioQQQ ,
"SanityCheck finds insanity so exiting\n");