cloudy  trunk
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 
23 namespace {
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 
46 STATIC void SanityCheckBegin(void );
47 STATIC void SanityCheckFinal(void );
48 /* chJob is either "begin" or "final"
49  * "begin is before code starts up
50  * "final is after model is complete */
51 void SanityCheck( const char *chJob )
52 {
53  DEBUG_ENTRY( "SanityCheck()" );
54 
55  if( strcmp(chJob,"begin") == 0 )
56  {
57  SanityCheckBegin();
58  }
59  else if( strcmp(chJob,"final") == 0 )
60  {
61  SanityCheckFinal();
62  }
63  else
64  {
65  fprintf(ioQQQ,"SanityCheck called with insane argument.\n");
66  cdEXIT(EXIT_FAILURE);
67  }
68 }
69 
70 STATIC void SanityCheckFinal(void )
71 {
72  /* PrtComment also has some ending checks on sanity */
73 }
74 
75 STATIC void SanityCheckBegin(void )
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;
573  Integrator<my_Integrand,Gaussian32> mySine;
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" );
662  cdEXIT(EXIT_FAILURE);
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" );
676  cdEXIT(EXIT_FAILURE);
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();
844  cdEXIT(EXIT_FAILURE);
845  }
846 }
H_Einstein_A
double H_Einstein_A(long int n, long int l, long int np, long int lp, long int iz)
Definition: hydro_bauman.cpp:2324
t_continuum::KshellLimit
long int KshellLimit
Definition: continuum.h:122
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
cont_gaunt_calc
double cont_gaunt_calc(double temp, double z, double photon)
Definition: cont_gaunt.cpp:26
Integrator
Definition: cddefines.h:1504
dense
t_dense dense
Definition: dense.cpp:24
t_opac::lgCaseB
bool lgCaseB
Definition: opacity.h:161
elementnames.h
cloudy.h
t_iso_sp::n_HighestResolved_max
long int n_HighestResolved_max
Definition: iso.h:505
rfield
t_rfield rfield
Definition: rfield.cpp:8
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
t_isoCTRL::chISO
const char * chISO[NISO]
Definition: iso.h:330
OpacityZero
void OpacityZero(void)
Definition: opacity_zero.cpp:8
realnum
float realnum
Definition: cddefines.h:103
rfield.h
SanityCheckBegin
STATIC void SanityCheckBegin(void)
Definition: sanity_check.cpp:75
STATIC
#define STATIC
Definition: cddefines.h:97
hydro_bauman.h
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
spsort
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1100
L_
#define L_(A_)
Definition: iso.h:21
thirdparty.h
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
getrs_wrapper
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
Definition: thirdparty_lapack.cpp:69
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
opac
t_opac opac
Definition: opacity.cpp:5
MIN2
#define MIN2
Definition: cddefines.h:761
ee1
double ee1(double x)
Definition: service.cpp:312
t_isoCTRL::lgNoRecombInterp
bool lgNoRecombInterp[NISO]
Definition: iso.h:385
PI
const UNUSED double PI
Definition: physconst.h:29
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_opac::lgRedoStatic
bool lgRedoStatic
Definition: opacity.h:147
expn
double expn(int n, double x)
Definition: thirdparty.cpp:2121
dense.h
trace
t_trace trace
Definition: trace.cpp:5
cddefines.h
Integrator::sum
double sum(double min, double max, Integrand func)
Definition: cddefines.h:1531
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
e2
double e2(double t)
Definition: service.cpp:299
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_opac::OpacStack
double * OpacStack
Definition: opacity.h:151
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
t_rfield::nflux
long int nflux
Definition: rfield.h:43
heavy.h
iso_recomb_check
double iso_recomb_check(long ipISO, long nelem, long level, double temperature)
Definition: iso_radiative_recomb.cpp:734
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:95
hypho.h
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
SanityCheck
void SanityCheck(const char *chJob)
Definition: sanity_check.cpp:51
yVector
static double * yVector
Definition: atom_feii.cpp:167
MAX2
#define MAX2
Definition: cddefines.h:782
LIMELM
const int LIMELM
Definition: cddefines.h:258
getrf_wrapper
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
Definition: thirdparty_lapack.cpp:47
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:487
t_rfield::nupper
long int nupper
Definition: rfield.h:46
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
hydrogenic.h
ipH2p
const int ipH2p
Definition: iso.h:29
t_opac::OpacStatic
double * OpacStatic
Definition: opacity.h:114
t_rfield::anu
double * anu
Definition: rfield.h:58
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
physconst.h
t_iso_sp::BranchRatio
multi_arr< double, 2 > BranchRatio
Definition: iso.h:451
t_rfield::widflx
realnum * widflx
Definition: rfield.h:65
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
iso_state_lifetime
double iso_state_lifetime(long ipISO, long nelem, long n, long l)
Definition: iso_create.cpp:1057
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
ShowMe
void ShowMe(void)
Definition: service.cpp:181
ipH1s
const int ipH1s
Definition: iso.h:27
t_phycon::sqrte
double sqrte
Definition: phycon.h:48
continuum
t_continuum continuum
Definition: continuum.cpp:5
TE1RYD
const UNUSED double TE1RYD
Definition: physconst.h:183
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
qg32
double qg32(double, double, double(*)(double))
Definition: service.cpp:1053
continuum.h
SanityCheckFinal
STATIC void SanityCheckFinal(void)
Definition: sanity_check.cpp:70
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
t_phycon::te
double te
Definition: phycon.h:11
NISO
const int NISO
Definition: cddefines.h:261
t_opac::ipElement
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:269
t_Heavy::ipHeavy
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
OpacityAdd1Element
void OpacityAdd1Element(long int ipZ)
Definition: opacity_add1element.cpp:12
xMatrix
static double ** xMatrix
Definition: atom_feii.cpp:170
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
H_photo_cs
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
Definition: hydro_bauman.cpp:564
N_
#define N_(A_)
Definition: iso.h:20
helike_recom.h
Generated by   doxygen 1.8.17