35 if( strcmp(chJob,
"begin") == 0 )
39 else if( strcmp(chJob,
"final") == 0 )
45 fprintf(
ioQQQ,
"SanityCheck called with insane argument.\n");
72 long n, logu, loggamma2;
96 for( nelem=0; nelem<
LIMELM; ++nelem )
106 fprintf(
ioQQQ,
" SanityCheck found insane H-like As.\n");
120 for( logu=-4; logu<=4; logu++)
123 for(loggamma2=-4; loggamma2<=4; loggamma2++)
125 double SutherlandGff[9][9]=
126 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
127 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
128 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
129 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
130 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
131 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
132 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
133 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
134 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};
137 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
139 if( error>0.11 || ( loggamma2<2 && error>0.05 ) )
141 fprintf(
ioQQQ,
" SanityCheck found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
142 logu, loggamma2, error );
154 for( nelem=1; nelem<
LIMELM; ++nelem )
159 0. ,9.47e+006 ,3.44e+008 ,1.74e+009 ,5.51e+009 ,1.34e+010 ,
160 2.79e+010 ,5.32E+010 ,8.81e+010 ,1.46E+011 ,2.15e+011 ,3.15e+011 ,
161 4.46e+011 ,6.39E+011 ,8.26e+011 ,1.09e+012 ,1.41e+012 ,1.86E+012 ,
162 2.26e+012 ,2.80e+012 ,3.44e+012 ,4.18e+012 ,5.04e+012 ,6.02e+012 ,
163 7.14e+012 ,8.40e+012 ,9.83e+012 ,1.14e+013 ,1.32e+013 ,1.52e+013
173 " SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n",
186 for( i = 0; i <=110; i++ )
188 double DrakeTotalAuls[111] = {
189 -1.0000E+00, -1.0000E+00, -1.0000E+00, 1.02160E+07,
190 1.02160E+07, 1.02160E+07, 1.80090E+09, 2.78530E+07,
191 1.82990E+07, 1.05480E+07, 7.07210E+07, 6.37210E+07,
192 5.79960E+08, 1.60330E+07, 1.13640E+07, 7.21900E+06,
193 3.11920E+07, 2.69830E+07, 1.38380E+07, 1.38330E+07,
194 2.52270E+08, 9.20720E+06, 6.82220E+06, 4.56010E+06,
195 1.64120E+07, 1.39290E+07, 7.16030E+06, 7.15560E+06,
196 4.25840E+06, 4.25830E+06, 1.31150E+08, 5.62960E+06,
197 4.29430E+06, 2.95570E+06, 9.66980E+06, 8.12340E+06,
198 4.19010E+06, 4.18650E+06, 2.48120E+06, 2.48120E+06,
199 1.64590E+06, 1.64590E+06, 7.65750E+07, 3.65330E+06,
200 2.84420E+06, 1.99470E+06, 6.16640E+06, 5.14950E+06,
201 2.66460E+06, 2.66200E+06, 1.57560E+06, 1.57560E+06,
202 1.04170E+06, 1.04170E+06, 7.41210E+05, 7.41210E+05,
203 4.84990E+07, 2.49130E+06, 1.96890E+06, 1.39900E+06,
204 4.16900E+06, 3.46850E+06, 1.79980E+06, 1.79790E+06,
205 1.06410E+06, 1.06410E+06, 7.02480E+05, 7.02480E+05,
206 4.98460E+05, 4.98460E+05, -1.0000E+00, -1.0000E+00,
207 3.26190E+07, 1.76920E+06, 1.41440E+06, 1.01460E+06,
208 2.94830E+06, 2.44680E+06, 1.27280E+06, 1.27140E+06,
209 7.52800E+05, 7.52790E+05, 4.96740E+05, 4.96740E+05,
210 3.51970E+05, 3.51970E+05, -1.0000E+00, -1.0000E+00,
211 -1.0000E+00, -1.0000E+00, 2.29740E+07, 1.29900E+06,
212 1.04800E+06, 7.57160E+05, 2.16090E+06, 1.79030E+06,
213 9.33210E+05, 9.32120E+05, 5.52310E+05, 5.52310E+05,
214 3.64460E+05, 3.64460E+05, 2.58070E+05, 2.58070E+05,
215 -1.0000E+00, -1.0000E+00, -1.0000E+00, -1.0000E+00,
216 -1.0000E+00, -1.0000E+00, 1.67840E+07};
218 if( DrakeTotalAuls[i] > 0. &&
221 if( fabs( DrakeTotalAuls[i] - (1./
StatesElem[
ipHE_LIKE][ipHELIUM][i].lifetime) ) /DrakeTotalAuls[i] > 0.001 )
224 " SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n",
242 const int NHE1CS = 20;
243 double he1cs[NHE1CS] =
245 5.480E-18 , 9.253E-18 , 1.598E-17 , 1.598E-17 , 1.598E-17 , 1.348E-17 ,
246 8.025E-18 , 1.449E-17 , 2.852E-17 , 1.848E-17 , 1.813E-17 , 2.699E-17 ,
247 1.077E-17 , 2.038E-17 , 4.159E-17 , 3.670E-17 , 3.575E-17 , 1.900E-17 ,
248 1.900E-17 , 4.175E-17
265 if( fabs( he1cs[n-1] -
opac.
OpacStack[i - 1] ) /he1cs[n-1] > 0.15 )
268 " SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n",
273 " n=%li, l=%li, s=%li\n",
294 " SanityCheck found insane1 %s %s recom coef: expected, n=%i error: %.2e \n",
307 " SanityCheck found insane2 %s %s recom coef: expected, n=%i error: %.2e \n",
323 const int NSORT = 100 ;
327 ipvector = (
long *)
MALLOC((NSORT)*
sizeof(
long int) );
331 for( i=0; i<NSORT; ++i )
346 if( lgFlag ) lgOK =
false;
348 for( i=1; i<NSORT; ++i )
352 if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] )
354 fprintf(
ioQQQ,
" SanityCheck found insane sort\n");
367 fprintf(
ioQQQ ,
"SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n",
386 # if !defined(NDEBUG)
393 fprintf(
ioQQQ,
" SanityCheck found insane widflx anu[i+1]=%e anu[i]=%e widflx=%e delta=%e rel err %e\n",
403 fprintf(
ioQQQ,
" SanityCheck found insane overlap1 widflx %e %e %e %e %e %e\n",
410 fprintf(
ioQQQ,
" big error was %e\n", x);
422 for( nelem=0; nelem<2; ++nelem )
429 Z4 = (double)(nelem+1);
435 if( fabs(ans1-ans2)/ans2 > 1e-3 )
437 fprintf(
ioQQQ ,
"SanityCheck finds insane A for H Lya %g %g nelem=%li\n",
438 ans1 , ans2 , nelem );
446 ans2 = 8.226*
powi( (1.+nelem) , 6 );
447 if( fabs(ans1-ans2)/ans2 > 1e-3 )
449 fprintf(
ioQQQ ,
"SanityCheck finds insane A for H 2-phot %g %g nelem=%li\n",
450 ans1 , ans2 , nelem );
472 ans2 = (2.532e6+2.20e6+2.70e6)*Z4;
473 if( fabs(ans1-ans2)/ans2 > 1e-2 )
476 "SanityCheck finds insane A for H5-2 found=%g correct=%g nelem=%li\n",
477 ans1 , ans2 , nelem );
486 for( nelem=0; nelem<
LIMELM; ++nelem )
494 for( ipLo=0; ipLo<ipHi; ++ipLo )
498 if( fabs(sum-1.)>0.01 )
501 "SanityCheck H branching ratio sum not unity for nelem=%li upper n=%i sum=%.3e\n",
524 anu[0] = (
realnum)(energy*1.01);
553 fprintf(
ioQQQ ,
"SanityCheck finds insane H photo cs\n");
572 ans2 =
expn( 1 , x );
573 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
575 fprintf(
ioQQQ ,
"SanityCheck finds insane E1 %g %g %g\n",
582 ans2 =
expn( 2 , x );
583 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
585 fprintf(
ioQQQ ,
"SanityCheck finds insane E2 %g %g %g\n",
626 A = (
double*)
MALLOC(
sizeof(
double)*NDIM*NDIM );
629 for( i=0; i < 3; ++i )
631 for( j=0; j < 3; ++j )
633 A[i*NDIM+j] = xMatrix[i][j];
644 fprintf(
ioQQQ,
" SanityCheck DGETRF error\n" );
658 fprintf(
ioQQQ,
" SanityCheck DGETRS error\n" );
672 x = fabs( yVector[i]-i-1.);
673 rcond =
MAX2( rcond, x );
677 if( rcond>DBL_EPSILON)
680 "SanityCheck found too large a deviation in matrix solver = %g \n",
691 for( nelem=0; nelem<
LIMELM; ++nelem )
721 for( ion=0; ion<=nelem; ++ion )
724 for( nshells=0; nshells<
Heavy.
nsShells[nelem][ion]; ++nshells )
734 "SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n",
735 nelem , ion , nshells, j );
737 "value was %li \n",
opac.
ipElement[nelem][ion][nshells][j] );
748 for( j=1; j <= (nelem + 2); j++ )
770 "SanityCheck found non-positive photo cs for nelem=%li ion=%li \n",
773 "value was %.2e + %.2e nelem %li ion %li at energy %.2e\n",
785 for( j=1; j <= (nelem + 2); j++ )
817 fprintf(
ioQQQ,
" SanityCheck returns OK\n");
825 fprintf(
ioQQQ ,
"SanityCheck finds insanity so exiting\n");