cloudy trunk
Loading...
Searching...
No Matches
sanity_check.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*SanityCheck, check that various parts of the code still work, called by Cloudy after continuum
4 * and optical depth arrays are set up, but before initial temperature and ionization */
5#include "cddefines.h"
6#include "physconst.h"
7#include "thirdparty.h"
8#include "dense.h"
9#include "elementnames.h"
10#include "continuum.h"
11#include "helike_recom.h"
12#include "rfield.h"
13#include "taulines.h"
14#include "hypho.h"
15#include "iso.h"
16#include "opacity.h"
17#include "hydro_bauman.h"
18#include "hydrogenic.h"
19#include "heavy.h"
20#include "trace.h"
21#include "cloudy.h"
22
23namespace {
24 class my_Integrand: public std::unary_function<double, double>
25 {
26 public:
27 double operator() (double x)
28 {
29 return sin(x);
30 }
31 // the constructor is needed to avoid warnings by the Sun Studio compiler
32 my_Integrand() {}
33 };
34
35 // create a version of sin with C++ linkage
36 // this is done because we need a C++ pointer to this routine below
37 double mySin(double x)
38 {
39 return sin(x);
40 }
41}
42
43/* NB - this routine must not change any global variables - any that are changed as part of
44 * a test must be reset, so that the code retains state */
45
46STATIC void SanityCheckBegin(void );
47STATIC void SanityCheckFinal(void );
48/* chJob is either "begin" or "final"
49 * "begin is before code starts up
50 * "final is after model is complete */
51void SanityCheck( const char *chJob )
52{
53 DEBUG_ENTRY( "SanityCheck()" );
54
55 if( strcmp(chJob,"begin") == 0 )
56 {
58 }
59 else if( strcmp(chJob,"final") == 0 )
60 {
62 }
63 else
64 {
65 fprintf(ioQQQ,"SanityCheck called with insane argument.\n");
67 }
68}
69
71{
72 /* PrtComment also has some ending checks on sanity */
73}
74
76{
77 bool lgOK=true;
78 int lgFlag;// error return for spsort, 0 success, >=1 for errors
79 int32 ner, ipiv[3];
80 long i ,
81 j ,
82 nelem ,
83 ion ,
84 nshells;
85 double *A;
86
87 /* this will be charge to the 4th power */
88 double Aul ,
89 error,
90 Z4, gaunt;
91
92 long n, logu, loggamma2;
93
94 const int NDIM = 10;
95 double x , ans1 , ans2 , xMatrix[NDIM][NDIM] , yVector[NDIM] ,
96 rcond;
97 realnum *fvector;
98 long int *ipvector;
99
100 DEBUG_ENTRY( "SanityCheck()" );
101
102 /*********************************************************
103 * *
104 * confirm that various part of cloudy still work *
105 * *
106 *********************************************************/
107
108 /* if this is no longer true at end, we have a problem */
109 lgOK = true;
110
111 /*********************************************************
112 * *
113 * check that all the Lyas As are ok *
114 * *
115 *********************************************************/
116 for( nelem=0; nelem<LIMELM; ++nelem )
117 {
118 /* this element may be turned off */
119 if( dense.lgElmtOn[nelem] )
120 {
121 /* H_Einstein_A( n, l, np, lp, iz ) - all are on physics scale */
122 Aul = H_Einstein_A( 2, 1, 1, 0, nelem+1 );
123 /*fprintf(ioQQQ,"%li\t%.4e\n", nelem+1, iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul() );*/
124 if( fabs(Aul - iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul() ) /Aul > 0.01 )
125 {
126 fprintf(ioQQQ," SanityCheck found insane H-like As.\n");
127 lgOK = false;
128 }
129 }
130 }
131
132 /*********************************************************
133 * *
134 * check that gaunt factors are good *
135 * *
136 *********************************************************/
137 /* Uncommenting each of the four print statements here
138 * will produce a nice table comparable to Sutherland 98, Table 2. */
139 /* fprintf(ioQQQ,"u\t-4\t-3\t-2\t-1\t0\t1\t2\t3\t4\n");*/
140 for( logu=-4; logu<=4; logu++)
141 {
142 /*fprintf(ioQQQ,"%li\t", logu);*/
143 for(loggamma2=-4; loggamma2<=4; loggamma2++)
144 {
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}};
155
156 gaunt = cont_gaunt_calc( TE1RYD/pow(10.,(double)loggamma2), 1., pow(10.,(double)(logu-loggamma2)) );
157 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
158 /*fprintf(ioQQQ,"%1.3f\t", gaunt);*/
159 if( error>0.11 || ( loggamma2<2 && error>0.05 ) )
160 {
161 fprintf(ioQQQ," SanityCheck found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
162 logu, loggamma2, error );
163 lgOK = false;
164 }
165 }
166 /*fprintf(ioQQQ,"\n");*/
167 }
168
169 /*********************************************************
170 * *
171 * check some transition probabililties for he-like ions *
172 * *
173 *********************************************************/
174 for( nelem=1; nelem<LIMELM; ++nelem )
175 {
176 /* the helike 9-1 transition, A(3^3P to 2^3S) */
177 double as[]={
178 /* updated with Johnson values */
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
184 };
185
186 if( iso_sp[ipHE_LIKE][nelem].numLevels_max > 8 && dense.lgElmtOn[nelem])
187 {
188 /* following used to print current values of As */
189 if( fabs( as[nelem] - iso_sp[ipHE_LIKE][nelem].trans(9,1).Emis().Aul() ) /as[nelem] > 0.025 )
190 {
191 fprintf(ioQQQ,
192 " SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n",
193 nelem,
194 as[nelem] ,
195 iso_sp[ipHE_LIKE][nelem].trans(9,1).Emis().Aul() );
196 lgOK = false;
197 }
198 }
199 }
200
201 /* only do this test if case b is not in effect */
202 if( !opac.lgCaseB )
203 {
204
205 for( i = 0; i <=110; i++ )
206 {
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};
236
237 if( DrakeTotalAuls[i] > 0. &&
238 i < iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max)
239 {
240 if( fabs( DrakeTotalAuls[i] - (1./iso_sp[ipHE_LIKE][ipHELIUM].st[i].lifetime()) ) /DrakeTotalAuls[i] > 0.001 )
241 {
242 fprintf(ioQQQ,
243 " SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n",
244 i,
245 DrakeTotalAuls[i],
246 (1./iso_sp[ipHE_LIKE][ipHELIUM].st[i].lifetime()) );
247 lgOK = false;
248 }
249 }
250 }
251 }
252
253 /*********************************************************
254 * *
255 * check the threshold photoionization cs for He I *
256 * *
257 *********************************************************/
258 if( dense.lgElmtOn[ipHELIUM] )
259 {
260 /* HeI photoionization cross sections at threshold for lowest 20 levels */
261 const int NHE1CS = 20;
262 double he1cs[NHE1CS] =
263 {
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
268 };
269
270 /* loop over levels and check on photo cross section */
271 j = MIN2( NHE1CS+1 , iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max -iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_max );
272 for( n=1; n<j; ++n )
273 {
274 /* above list of levels does not include the ground */
275 i = iso_sp[ipHE_LIKE][ipHELIUM].fb[n].ipOpac;
276 ASSERT( i>0 );
277 /*fprintf(ioQQQ,"%li\t%lin", n , i );*/
278 /* >>chng 02 apr 10, from 0.01 to 0.02, values stored
279 * where taken from calc at low contin resolution, when continuum
280 * resolution changed this changes too */
281 /*fprintf(ioQQQ,"%li %.2e\n", n,( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] );*/
282 /* >>chng 02 jul 16, limt from 0.02 to 0.04, so that "set resolution 4" will work */
283 /* >>chng 04 may 18, levels 10 and 11 are about 12% off - because of energy binning, chng from 0.08 to 0.15 */
284 if( fabs( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] > 0.15 )
285 {
286 fprintf(ioQQQ,
287 " SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n",
288 n,
289 he1cs[n-1] ,
290 opac.OpacStack[i - 1]);
291 fprintf(ioQQQ,
292 " n=%li, l=%li, s=%li\n",
293 iso_sp[ipHE_LIKE][ipHELIUM].st[n].n() ,
294 iso_sp[ipHE_LIKE][ipHELIUM].st[n].l() ,
295 iso_sp[ipHE_LIKE][ipHELIUM].st[n].S());
296 lgOK = false;
297 }
298 }
299 }
300
301 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
302 {
303 long nelem = ipISO;
304 /* Check for agreement between on-the-fly and interpolation calculations
305 * of recombination, but only if interpolation is turned on. */
306 if( !iso_ctrl.lgNoRecombInterp[ipISO] )
307 {
308 /* check the recombination coefficients for ground state */
309 error = fabs( iso_recomb_check( ipISO, nelem , 0 , 7500. ) );
310 if( error > 0.01 )
311 {
312 fprintf(ioQQQ,
313 " SanityCheck found insane1 %s %s recom coef: expected, n=%i error: %.2e \n",
314 iso_ctrl.chISO[ipISO],
315 elementnames.chElementSym[nelem],
316 0,
317 error );
318 lgOK = false;
319 }
320
321 /* check the recombination coefficients for ground state of the root of each iso sequence */
322 error = fabs( iso_recomb_check( ipISO, nelem , 1 , 12500. ) );
323 if( error > 0.01 )
324 {
325 fprintf(ioQQQ,
326 " SanityCheck found insane2 %s %s recom coef: expected, n=%i error: %.2e \n",
327 iso_ctrl.chISO[ipISO],
328 elementnames.chElementSym[nelem],
329 1,
330 error );
331 lgOK = false;
332 }
333 }
334 }
335
336 /*********************************************************
337 * *
338 * check out the sorting routine *
339 * *
340 *********************************************************/
341
342 const int NSORT = 100 ;
343
344 fvector = (realnum *)MALLOC((NSORT)*sizeof(realnum) );
345
346 ipvector = (long *)MALLOC((NSORT)*sizeof(long int) );
347
348 nelem = 1;
349 /* make up some unsorted values */
350 for( i=0; i<NSORT; ++i )
351 {
352 nelem *= -1;
353 fvector[i] = (realnum)nelem * ((realnum)NSORT-i);
354 }
355
356 /*spsort netlib routine to sort array returning sorted indices */
357 spsort(fvector,
358 NSORT,
359 ipvector,
360 /* flag saying what to do - 1 sorts into increasing order, not changing
361 * the original routine */
362 1,
363 &lgFlag);
364
365 if( lgFlag ) lgOK = false;
366
367 for( i=1; i<NSORT; ++i )
368 {
369 /*fprintf(ioQQQ," %li %li %.0f\n",
370 i, ipvector[i],fvector[ipvector[i]] );*/
371 if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] )
372 {
373 fprintf(ioQQQ," SanityCheck found insane sort\n");
374 lgOK = false;
375 }
376 }
377
378 free( fvector );
379 free( ipvector);
380
381# if 0
382 ttemp = (realnum)sqrt(phycon.te);
383 /* check that the temperatures make sense */
384 if( fabs(ttemp - phycon.sqrte )/ttemp > 1e-5 )
385 {
386 fprintf(ioQQQ , "SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n",
387 phycon.te ,
388 sqrt(phycon.te) ,
389 phycon.sqrte ,
390 fabs(sqrt(phycon.te) - phycon.sqrte ) );
391 lgOK = false;
392 }
393# endif
394
395 /*********************************************************
396 * *
397 * confirm that widflx and anu arrays correspond *
398 * to one another *
399 * *
400 *********************************************************/
401
402# if 0
403 /* this check on widflx can't be used since some sharpling curved continua, like laser,
404 * totally fail due to non-linear nature of widflx and anu relationship */
405# if !defined(NDEBUG)
406 x = 0.;
407 for( i=1; i<rfield.nupper-1; ++i )
408 {
409 if( fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.) > 0.02 )
410 {
411 ans1 = fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.);
412 fprintf(ioQQQ," SanityCheck found insane widflx anu[i+1]=%e anu[i]=%e widflx=%e delta=%e rel err %e\n",
413 rfield.anu[i+1] , rfield.anu[i] , rfield.widflx[i] , rfield.anu[i+1] -rfield.anu[i] , ans1 );
414 lgOK = false;
415 x = MAX2( ans1 , x);
416 }
417 /* problems when at energy where resolution of grid changes dramatically */
418 /* this is resolution at current energy */
419 ans1 = rfield.widflx[i] / rfield.anu[i];
420 if( (rfield.anu[i]+rfield.widflx[i]/2.)*(1.-ans1/10.) > rfield.anu[i+1] - rfield.widflx[i+1]/2.)
421 {
422 fprintf(ioQQQ," SanityCheck found insane overlap1 widflx %e %e %e %e %e %e\n",
423 rfield.anu[i] , rfield.widflx[i], rfield.anu[i] + rfield.widflx[i]/2. , rfield.anu[i+1],
424 rfield.widflx[i+1], rfield.anu[i+1] -rfield.widflx[i+1]/2. );
425 lgOK = false;
426 }
427 if( !lgOK )
428 {
429 fprintf(ioQQQ," big error was %e\n", x);
430 }
431 }
432# endif
433# endif
434
435
436 /*********************************************************
437 * *
438 * confirm that hydrogen einstein As are still valid *
439 * *
440 *********************************************************/
441 for( nelem=0; nelem<2; ++nelem )
442 {
443 /* this element may be turned off */
444 if( dense.lgElmtOn[nelem] )
445 {
446 /*Z4 = (double)(POW2(nelem+1)*POW2(nelem+1));*/
447 /* form charge to the 4th power */
448 Z4 = (double)(nelem+1);
449 Z4 *= Z4;
450 Z4 *= Z4;
451 /* H Lya */
452 ans1 = iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul();
453 ans2 = 6.265e8*Z4;
454 if( fabs(ans1-ans2)/ans2 > 1e-3 )
455 {
456 fprintf(ioQQQ , "SanityCheck finds insane A for H Lya %g %g nelem=%li\n",
457 ans1 , ans2 , nelem );
458 lgOK = false;
459 }
460 }
461 }
462
463 /* check that hydrogenic branching ratios add up to unity */
464 for( nelem=0; nelem<LIMELM; ++nelem )
465 {
466 if( dense.lgElmtOn[nelem] )
467 {
468 int ipHi, ipLo;
469 for( ipHi=4; ipHi< iso_sp[ipH_LIKE][nelem].numLevels_max-iso_sp[ipH_LIKE][nelem].nCollapsed_max; ++ipHi )
470 {
471 double sum = 0.;
472 for( ipLo=0; ipLo<ipHi; ++ipLo )
473 {
474 sum += iso_sp[ipH_LIKE][nelem].BranchRatio[ipHi][ipLo];
475 }
476 if( fabs(sum-1.)>0.01 )
477 {
478 fprintf(ioQQQ ,
479 "SanityCheck H branching ratio sum not unity for nelem=%li upper level=%i sum=%.3e\n",
480 nelem, ipHi, sum );
481 lgOK = false;
482 }
483 }
484 }
485 }
486
487 /* check that hydrogenic lifetimes are sane (compare inverse sum of As with closed form of lifetime) */
488 for( nelem=0; nelem<LIMELM; ++nelem )
489 {
490 if( dense.lgElmtOn[nelem] )
491 {
492 int ipHi, ipLo;
493 for( ipHi=1; ipHi< iso_sp[ipH_LIKE][nelem].numLevels_max-iso_sp[ipH_LIKE][nelem].nCollapsed_max; ++ipHi )
494 {
495 double inverse_sum = 0.;
496 double sum = 0.;
497 long ipISO = ipH_LIKE;
498
499 /* we do not have an accurate closed form for l=0 lifetimes. Everything else should be very accurate. */
500 if( L_(ipHi)==0 )
501 continue;
502
503 for( ipLo=0; ipLo<ipHi; ++ipLo )
504 {
505 sum += iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul();
506 }
507
508 inverse_sum = 1./sum;
509
510 double lifetime = iso_state_lifetime( ipH_LIKE, nelem, N_(ipHi), L_(ipHi) );
511 double percent_error = (1. - inverse_sum/lifetime)*100;
512 /* The closed form seems to consistently yield lifetimes roughly 0.2% smaller than the inverse sum
513 * Transition probabilities go like energy cubed. The cube of the finite mass rydberg is about 0.16% less than unity.
514 * is that the difference? */
515 /* for now enforce to better than 0.5% */
516 if( fabs(percent_error) > 0.5 )
517 {
518 fprintf(ioQQQ ,
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 );
521 lgOK = false;
522 }
523 }
524 }
525 }
526
527
528 /* check photo cross sections for H */
529 long ipISO = ipH_LIKE;
530 nelem = 0;
531 for( long n=1; n <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++n )
532 {
533 // This sanity check is of little value, as the cross-section is calculated
534 // with the same routine here as when populating the opacity stack.
535 // However, the cell mid-point energy used in populating OpacStack can be
536 // greater than the energy implied by the relative energy below.
537 // For very Rydberg, very yrare levels, the cross-section can fall off fast
538 // enough over that energy difference to fail this check.
539 // Will reported a sim that failed the check for some levels with n >= 71.
540 // So just don't bother checking above 60.
541 if( n > 60 )
542 break;
543
544 double rel_photon_energy = 1. + FLT_EPSILON*2.;
545 for( long l=0; l < n; l++ )
546 {
547 double cs;
548 long index = iso_sp[ipISO][nelem].QuantumNumbers2Index[n][l][2];
549
550 cs = H_photo_cs( rel_photon_energy, n, l, nelem+1 );
551
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.);
554
555 if( error > 0.05 )
556 {
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 );
559 lgOK = false;
560 }
561 }
562 }
563
564 /*********************************************************
565 * *
566 * confirm that Gaussian integration routines still work *
567 * *
568 *********************************************************/
569 ASSERT( fp_equal( qg32( 0., PI, mySin ) , 2. ) );
570
571 /* And again with the new structure */
572 my_Integrand func;
574 ASSERT( fp_equal( mySine.sum( 0., PI, func ), 2. ) );
575
576 /*********************************************************
577 * *
578 * confirm that exponential integral routines still work *
579 * *
580 *********************************************************/
581
582 /* check that first and second exponential integrals are ok,
583 * step through range of values, beginning with following */
584 x = 1e-3;
585 do
586 {
587 /* check that fast e1 routine is ok */
588 ans1 = ee1(x);
589 ans2 = expn( 1 , x );
590 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
591 {
592 fprintf(ioQQQ , "SanityCheck finds insane E1 %g %g %g\n",
593 x , ans1 , ans2 );
594 lgOK = false;
595 }
596
597 /* check that e2 is ok */
598 ans1 = e2(x);
599 ans2 = expn( 2 , x );
600 if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
601 {
602 fprintf(ioQQQ , "SanityCheck finds insane E2 %g %g %g\n",
603 x , ans1 , ans2 );
604 lgOK = false;
605 }
606
607 /* now increment x */
608 x *= 2.;
609 /* following limit set by sexp returning zero, used in ee1 */
610 } while( x < 64. );
611
612 /*********************************************************
613 * *
614 * confirm that matrix inversion routine still works *
615 * *
616 *********************************************************/
617
618 /* these are the answer, chosen to get xvec 1,2,3 */
619 yVector[0] = 1.;
620 yVector[1] = 3.;
621 yVector[2] = 3.;
622
623 /* zero out the main matrix */
624 for(i=0;i<3;++i)
625 {
626 for( j=0;j<3;++j )
627 {
628 xMatrix[i][j] = 0.;
629 }
630 }
631
632 /* remember that order is column, row, alphabetical order, rc */
633 xMatrix[0][0] = 1.;
634 xMatrix[0][1] = 1.;
635 xMatrix[1][1] = 1.;
636 xMatrix[2][2] = 1.;
637
638 /* this is the default matrix solver */
639 /* this test is the 1-d matrix with 2-d macro simulation */
640 /* LDA is right dimension of matrix */
641
642 /* MALLOC space for the 1-d array */
643 A = (double*)MALLOC( sizeof(double)*NDIM*NDIM );
644
645 /* copy over the main matrix */
646 for( i=0; i < 3; ++i )
647 {
648 for( j=0; j < 3; ++j )
649 {
650 A[i*NDIM+j] = xMatrix[i][j];
651 }
652 }
653
654 ner = 0;
655
656 /*void DGETRF(long,long,double*,long,long[],long*);*/
657 /*void DGETRF(int,int,double*,int,int[],int*);*/
658 getrf_wrapper(3, 3, A, NDIM, ipiv, &ner);
659 if( ner != 0 )
660 {
661 fprintf( ioQQQ, " SanityCheck DGETRF error\n" );
663 }
664
665 /* usage DGETRS, 'N' = no transpose
666 * order of matrix,
667 * number of cols in bvec, =1
668 * array
669 * leading dim of array */
670 /*void DGETRS(char,int,int,double*,int,int[],double*,int,int*);*/
671 getrs_wrapper('N', 3, 1, A, NDIM, ipiv, yVector, 3, &ner);
672
673 if( ner != 0 )
674 {
675 fprintf( ioQQQ, " SanityCheck DGETRS error\n" );
677 }
678 /* release the vector */
679 free( A );
680
681 /* now check on validity of the solution, demand that this
682 * simple problem have come within a few epsilons of the
683 * correct answer */
684
685 /* find largest deviation */
686 rcond = 0.;
687 for(i=0;i<3;++i)
688 {
689 x = fabs( yVector[i]-i-1.);
690 rcond = MAX2( rcond, x );
691 /*printf(" %g ", yVector[i]);*/
692 }
693
694 if( rcond>DBL_EPSILON)
695 {
696 fprintf(ioQQQ,
697 "SanityCheck found too large a deviation in matrix solver = %g \n",
698 rcond);
699 /* set flag saying that things are not ok */
700 lgOK = false;
701 }
702 /* end matrix inversion check */
703
704 /* these pointers were set to INT_MIN in ContCreatePointers,
705 * then set to valid numbers in ipShells and OpacityCreate1Element
706 * this checks that all values have been properly filled */
707 for( nelem=0; nelem<LIMELM; ++nelem )
708 {
709 /* must reset state of code after tests performed, remember state here */
710 realnum xIonF[NISO][LIMELM];
711 double hbn[NISO][LIMELM], hn[NISO][LIMELM];
712
713 if( dense.lgElmtOn[nelem] )
714 {
715 /* set these abundances so that opacities can be checked below */
716 hbn[ipH_LIKE][nelem] = iso_sp[ipH_LIKE][nelem].fb[0].DepartCoef;
717 hn[ipH_LIKE][nelem] = iso_sp[ipH_LIKE][nelem].st[0].Pop();
718 xIonF[ipH_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipH_LIKE];
719
720 iso_sp[ipH_LIKE][nelem].fb[0].DepartCoef = 0.;
721 iso_sp[ipH_LIKE][nelem].st[0].Pop() = 1.;
722 dense.xIonDense[nelem][nelem+1-ipH_LIKE] = 1.;
723
724 if( nelem > ipHYDROGEN )
725 {
726
727 hbn[ipHE_LIKE][nelem] = iso_sp[ipHE_LIKE][nelem].fb[0].DepartCoef;
728 hn[ipHE_LIKE][nelem] = iso_sp[ipHE_LIKE][nelem].st[0].Pop();
729 xIonF[ipHE_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipHE_LIKE];
730
731 /* this does not exist for hydrogen itself */
732 iso_sp[ipHE_LIKE][nelem].fb[0].DepartCoef = 0.;
733 iso_sp[ipHE_LIKE][nelem].st[0].Pop() = 1.;
734 dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = 1.;
735 }
736
737 for( ion=0; ion<=nelem; ++ion )
738 {
739 /* loop over all shells that are defined */
740 for( nshells=0; nshells<Heavy.nsShells[nelem][ion]; ++nshells )
741 {
742 for( j=0; j<3; ++j )
743 {
744 /* >>chng 00 apr 05, array index is on fortran scale so must be
745 * >= 1. This test had been <0, correct for C. Caught by Peter van Hoof */
746 if( opac.ipElement[nelem][ion][nshells][j] <=0 )
747 {
748 /* this is not possible */
749 fprintf(ioQQQ,
750 "SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n",
751 nelem , ion , nshells, j );
752 fprintf(ioQQQ,
753 "value was %li \n", opac.ipElement[nelem][ion][nshells][j] );
754 /* set flag saying that things are not ok */
755 lgOK = false;
756 }
757 }
758 }
759
760 if( nelem > 1 )
761 {
762 vector<realnum> saveion;
763 saveion.resize(nelem+2);
764 /* check that photoionization cross sections are ok */
765 for( j=0; j <= (nelem + 1); j++ )
766 {
767 saveion[j] = dense.xIonDense[nelem][j];
768 dense.xIonDense[nelem][j] = 0.;
769 }
770
771 dense.xIonDense[nelem][ion] = 1.;
772
773 OpacityZero();
774 opac.lgRedoStatic = true;
775
776 /* generate opacity with standard routine - this is the one
777 * called in OpacityAddTotal to make opacities in usual calculations */
778 OpacityAdd1Element(nelem);
779
780 /* this starts one beyond energy of threshold since cs may be zero there */
781 for( j=Heavy.ipHeavy[nelem][ion]; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
782 {
783 if( opac.opacity_abs[j]+opac.OpacStatic[j] < FLT_MIN )
784 {
785 /* this is not possible */
786 fprintf(ioQQQ,
787 "SanityCheck found non-positive photo cs for nelem=%li ion=%li \n",
788 nelem , ion );
789 fprintf(ioQQQ,
790 "value was %.2e + %.2e nelem %li ion %li at energy %.2e\n",
791 opac.opacity_abs[j] ,
792 opac.OpacStatic[j] ,
793 nelem ,
794 ion ,
795 rfield.anu[j]);
796 /* set flag saying that things are not ok */
797 lgOK = false;
798 break;/**/
799 }
800 }
801 /* reset the ionization distribution */
802 for( j=0; j <= (nelem + 1); j++ )
803 {
804 dense.xIonDense[nelem][j] = saveion[j];
805 }
806
807 }
808 }
809 iso_sp[ipH_LIKE][nelem].fb[ipH1s].DepartCoef = hbn[ipH_LIKE][nelem];
810 iso_sp[ipH_LIKE][nelem].st[ipH1s].Pop() = hn[ipH_LIKE][nelem];
811 dense.xIonDense[nelem][nelem+1-ipH_LIKE] = xIonF[ipH_LIKE][nelem];
812
813 if( nelem > ipHYDROGEN )
814 {
815 iso_sp[ipHE_LIKE][nelem].fb[ipHe1s1S].DepartCoef = hbn[ipHE_LIKE][nelem];
816 iso_sp[ipHE_LIKE][nelem].st[ipHe1s1S].Pop() = hn[ipHE_LIKE][nelem];
817 dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = xIonF[ipHE_LIKE][nelem];
818 }
819 }
820 }
821
822
823 /*********************************************************
824 * *
825 * everything is done, all checks make, did we pass them?*
826 * *
827 *********************************************************/
828
829 if( lgOK )
830 {
831 /*return if ok */
832 if( trace.lgTrace )
833 {
834 fprintf( ioQQQ, " SanityCheck returns OK\n");
835 }
836 return;
837 }
838
839 else
840 {
841 /* stop since problem encountered, lgEOF set false */
842 fprintf(ioQQQ , "SanityCheck finds insanity so exiting\n");
843 ShowMe();
845 }
846}
static double ** xMatrix
static double * yVector
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
double qg32(double, double, double(*)(double))
Definition service.cpp:1053
double ee1(double x)
Definition service.cpp:312
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
double e2(double t)
Definition service.cpp:299
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition service.cpp:1100
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
double sum(double min, double max, Integrand func)
Definition cddefines.h:1531
double cont_gaunt_calc(double temp, double z, double photon)
t_continuum continuum
Definition continuum.cpp:5
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_Heavy Heavy
Definition heavy.cpp:5
double H_Einstein_A(long int n, long int l, long int np, long int lp, long int iz)
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
const int ipH1s
Definition iso.h:27
#define N_(A_)
Definition iso.h:20
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
double iso_state_lifetime(long ipISO, long nelem, long n, long l)
double iso_recomb_check(long ipISO, long nelem, long level, double temperature)
const int ipH_LIKE
Definition iso.h:62
#define L_(A_)
Definition iso.h:21
t_opac opac
Definition opacity.cpp:5
void OpacityZero(void)
void OpacityAdd1Element(long int ipZ)
#define S(I_, J_)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double PI
Definition physconst.h:29
UNUSED const double TE1RYD
Definition physconst.h:183
t_rfield rfield
Definition rfield.cpp:8
STATIC void SanityCheckFinal(void)
STATIC void SanityCheckBegin(void)
void SanityCheck(const char *chJob)
double expn(int n, double x)
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
t_trace trace
Definition trace.cpp:5