cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sanity_check.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 /* NB - this routine must not change any global variables - any that are changed as part of
24  * a test must be reset, so that the code retains state */
25 
26 STATIC void SanityCheckBegin(void );
27 STATIC void SanityCheckFinal(void );
28 /* chJob is either "begin" or "final"
29  * "begin is before code starts up
30  * "final is after model is complete */
31 void SanityCheck( const char *chJob )
32 {
33  DEBUG_ENTRY( "SanityCheck()" );
34 
35  if( strcmp(chJob,"begin") == 0 )
36  {
38  }
39  else if( strcmp(chJob,"final") == 0 )
40  {
42  }
43  else
44  {
45  fprintf(ioQQQ,"SanityCheck called with insane argument.\n");
46  cdEXIT(EXIT_FAILURE);
47  }
48 }
49 
51 {
52  /* PrtComment also has some ending checks on sanity */
53 }
54 
56 {
57  bool lgOK=true;
58  int lgFlag;// error return for spsort, 0 success, >=1 for errors
59  int32 ner, ipiv[3];
60  long i ,
61  j ,
62  nelem ,
63  ion ,
64  nshells;
65  double *A;
66 
67  /* this will be charge to the 4th power */
68  double Aul ,
69  error,
70  Z4, gaunt;
71 
72  long n, logu, loggamma2;
73 
74  const int NDIM = 10;
75  double x , ans1 , ans2 , xMatrix[NDIM][NDIM] , yVector[NDIM] ,
76  rcond;
77  realnum *fvector;
78  long int *ipvector;
79 
80  DEBUG_ENTRY( "SanityCheck()" );
81 
82  /*********************************************************
83  * *
84  * confirm that various part of cloudy still work *
85  * *
86  *********************************************************/
87 
88  /* if this is no longer true at end, we have a problem */
89  lgOK = true;
90 
91  /*********************************************************
92  * *
93  * check that all the Lyas As are ok *
94  * *
95  *********************************************************/
96  for( nelem=0; nelem<LIMELM; ++nelem )
97  {
98  /* this element may be turned off */
99  if( dense.lgElmtOn[nelem] )
100  {
101  /* H_Einstein_A( n, l, np, lp, iz ) - all are on physics scale */
102  Aul = H_Einstein_A( 2, 1, 1, 0, nelem+1 );
103  /*fprintf(ioQQQ,"%li\t%.4e\n", nelem+1, Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul );*/
104  if( fabs(Aul - Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul ) /Aul > 0.01 )
105  {
106  fprintf(ioQQQ," SanityCheck found insane H-like As.\n");
107  lgOK = false;
108  }
109  }
110  }
111 
112  /*********************************************************
113  * *
114  * check that gaunt factors are good *
115  * *
116  *********************************************************/
117  /* Uncommenting each of the four print statements here
118  * will produce a nice table comparable to Sutherland 98, Table 2. */
119  /* fprintf(ioQQQ,"u\t-4\t-3\t-2\t-1\t0\t1\t2\t3\t4\n");*/
120  for( logu=-4; logu<=4; logu++)
121  {
122  /*fprintf(ioQQQ,"%li\t", logu);*/
123  for(loggamma2=-4; loggamma2<=4; loggamma2++)
124  {
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}};
135 
136  gaunt = cont_gaunt_calc( TE1RYD/pow(10.,(double)loggamma2), 1., pow(10.,(double)(logu-loggamma2)) );
137  error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
138  /*fprintf(ioQQQ,"%1.3f\t", gaunt);*/
139  if( error>0.11 || ( loggamma2<2 && error>0.05 ) )
140  {
141  fprintf(ioQQQ," SanityCheck found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
142  logu, loggamma2, error );
143  lgOK = false;
144  }
145  }
146  /*fprintf(ioQQQ,"\n");*/
147  }
148 
149  /*********************************************************
150  * *
151  * check some transition probabililties for he-like ions *
152  * *
153  *********************************************************/
154  for( nelem=1; nelem<LIMELM; ++nelem )
155  {
156  /* the helike 9-1 transition, A(3^3P to 2^3S) */
157  double as[]={
158  /* updated with Johnson values */
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
164  };
165 
166  if( iso.numLevels_max[ipHE_LIKE][nelem] > 8 && dense.lgElmtOn[nelem])
167  {
168  /* following used to print current values of As */
169  /*fprintf(ioQQQ,"%.2e\n", Transitions[ipHE_LIKE][nelem][9][1].Emis->Aul );*/
170  if( fabs( as[nelem] - Transitions[ipHE_LIKE][nelem][9][1].Emis->Aul ) /as[nelem] > 0.025 )
171  {
172  fprintf(ioQQQ,
173  " SanityCheck found insane He-like As: expected, nelem=%li found: %.2e %.2e.\n",
174  nelem,
175  as[nelem] ,
176  Transitions[ipHE_LIKE][nelem][9][1].Emis->Aul );
177  lgOK = false;
178  }
179  }
180  }
181 
182  /* only do this test if case b is not in effect */
183  if( !opac.lgCaseB )
184  {
185 
186  for( i = 0; i <=110; i++ )
187  {
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};
217 
218  if( DrakeTotalAuls[i] > 0. &&
220  {
221  if( fabs( DrakeTotalAuls[i] - (1./StatesElem[ipHE_LIKE][ipHELIUM][i].lifetime) ) /DrakeTotalAuls[i] > 0.001 )
222  {
223  fprintf(ioQQQ,
224  " SanityCheck found helium lifetime outside 0.1 pct of Drake values: index, expected, found: %li %.4e %.4e\n",
225  i,
226  DrakeTotalAuls[i],
227  (1./StatesElem[ipHE_LIKE][ipHELIUM][i].lifetime) );
228  lgOK = false;
229  }
230  }
231  }
232  }
233 
234  /*********************************************************
235  * *
236  * check the threshold photoionization cs for He I *
237  * *
238  *********************************************************/
239  if( dense.lgElmtOn[ipHELIUM] )
240  {
241  /* HeI photoionization cross sections at threshold for lowest 20 levels */
242  const int NHE1CS = 20;
243  double he1cs[NHE1CS] =
244  {
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
249  };
250 
251  /* loop over levels and check on photo cross section */
253  for( n=1; n<j; ++n )
254  {
255  /* above list of levels does not include the ground */
256  i = iso.ipOpac[ipHE_LIKE][ipHELIUM][n];
257  ASSERT( i>0 );
258  /*fprintf(ioQQQ,"%li\t%lin", n , i );*/
259  /* >>chng 02 apr 10, from 0.01 to 0.02, values stored
260  * where taken from calc at low contin resolution, when continuum
261  * resolution changed this changes too */
262  /*fprintf(ioQQQ,"%li %.2e\n", n,( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] );*/
263  /* >>chng 02 jul 16, limt from 0.02 to 0.04, so that "set resolution 4" will work */
264  /* >>chng 04 may 18, levels 10 and 11 are about 12% off - because of energy binning, chng from 0.08 to 0.15 */
265  if( fabs( he1cs[n-1] - opac.OpacStack[i - 1] ) /he1cs[n-1] > 0.15 )
266  {
267  fprintf(ioQQQ,
268  " SanityCheck found insane HeI photo cs: expected, n=%li found: %.3e %.3e.\n",
269  n,
270  he1cs[n-1] ,
271  opac.OpacStack[i - 1]);
272  fprintf(ioQQQ,
273  " n=%li, l=%li, s=%li\n",
274  StatesElem[ipHE_LIKE][ipHELIUM][n].n ,
275  StatesElem[ipHE_LIKE][ipHELIUM][n].l ,
276  StatesElem[ipHE_LIKE][ipHELIUM][n].S);
277  lgOK = false;
278  }
279  }
280  }
281 
282  for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
283  {
284  long nelem = ipISO;
285  /* Check for agreement between on-the-fly and interpolation calculations
286  * of recombination, but only if interpolation is turned on. */
287  if( !iso.lgNoRecombInterp[ipISO] )
288  {
289  /* check the recombination coefficients for ground state */
290  error = fabs( iso_recomb_check( ipISO, nelem , 0 , 7500. ) );
291  if( error > 0.01 )
292  {
293  fprintf(ioQQQ,
294  " SanityCheck found insane1 %s %s recom coef: expected, n=%i error: %.2e \n",
295  iso.chISO[ipISO],
296  elementnames.chElementSym[nelem],
297  0,
298  error );
299  lgOK = false;
300  }
301 
302  /* check the recombination coefficients for ground state of the root of each iso sequence */
303  error = fabs( iso_recomb_check( ipISO, nelem , 1 , 12500. ) );
304  if( error > 0.01 )
305  {
306  fprintf(ioQQQ,
307  " SanityCheck found insane2 %s %s recom coef: expected, n=%i error: %.2e \n",
308  iso.chISO[ipISO],
309  elementnames.chElementSym[nelem],
310  1,
311  error );
312  lgOK = false;
313  }
314  }
315  }
316 
317  /*********************************************************
318  * *
319  * check out the sorting routine *
320  * *
321  *********************************************************/
322 
323  const int NSORT = 100 ;
324 
325  fvector = (realnum *)MALLOC((NSORT)*sizeof(realnum) );
326 
327  ipvector = (long *)MALLOC((NSORT)*sizeof(long int) );
328 
329  nelem = 1;
330  /* make up some unsorted values */
331  for( i=0; i<NSORT; ++i )
332  {
333  nelem *= -1;
334  fvector[i] = (realnum)nelem * ((realnum)NSORT-i);
335  }
336 
337  /*spsort netlib routine to sort array returning sorted indices */
338  spsort(fvector,
339  NSORT,
340  ipvector,
341  /* flag saying what to do - 1 sorts into increasing order, not changing
342  * the original routine */
343  1,
344  &lgFlag);
345 
346  if( lgFlag ) lgOK = false;
347 
348  for( i=1; i<NSORT; ++i )
349  {
350  /*fprintf(ioQQQ," %li %li %.0f\n",
351  i, ipvector[i],fvector[ipvector[i]] );*/
352  if( fvector[ipvector[i]] <= fvector[ipvector[i-1]] )
353  {
354  fprintf(ioQQQ," SanityCheck found insane sort\n");
355  lgOK = false;
356  }
357  }
358 
359  free( fvector );
360  free( ipvector);
361 
362 # if 0
363  ttemp = (realnum)sqrt(phycon.te);
364  /* check that the temperatures make sense */
365  if( fabs(ttemp - phycon.sqrte )/ttemp > 1e-5 )
366  {
367  fprintf(ioQQQ , "SanityCheck finds insane te %e sqrt te %e sqrte %e dif %e\n",
368  phycon.te ,
369  sqrt(phycon.te) ,
370  phycon.sqrte ,
371  fabs(sqrt(phycon.te) - phycon.sqrte ) );
372  lgOK = false;
373  }
374 # endif
375 
376  /*********************************************************
377  * *
378  * confirm that widflx and anu arrays correspond *
379  * to one another *
380  * *
381  *********************************************************/
382 
383 # if 0
384  /* this check on widflx can't be used since some sharpling curved continua, like laser,
385  * totally fail due to non-linear nature of widflx and anu relationship */
386 # if !defined(NDEBUG)
387  x = 0.;
388  for( i=1; i<rfield.nupper-1; ++i )
389  {
390  if( fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.) > 0.02 )
391  {
392  ans1 = fabs( ((rfield.anu[i+1]-rfield.anu[i]) + (rfield.anu[i]-rfield.anu[i-1])) /rfield.widflx[i] /2.-1.);
393  fprintf(ioQQQ," SanityCheck found insane widflx anu[i+1]=%e anu[i]=%e widflx=%e delta=%e rel err %e\n",
394  rfield.anu[i+1] , rfield.anu[i] , rfield.widflx[i] , rfield.anu[i+1] -rfield.anu[i] , ans1 );
395  lgOK = false;
396  x = MAX2( ans1 , x);
397  }
398  /* problems when at energy where resolution of grid changes dramatically */
399  /* this is resolution at current energy */
400  ans1 = rfield.widflx[i] / rfield.anu[i];
401  if( (rfield.anu[i]+rfield.widflx[i]/2.)*(1.-ans1/10.) > rfield.anu[i+1] - rfield.widflx[i+1]/2.)
402  {
403  fprintf(ioQQQ," SanityCheck found insane overlap1 widflx %e %e %e %e %e %e\n",
404  rfield.anu[i] , rfield.widflx[i], rfield.anu[i] + rfield.widflx[i]/2. , rfield.anu[i+1],
405  rfield.widflx[i+1], rfield.anu[i+1] -rfield.widflx[i+1]/2. );
406  lgOK = false;
407  }
408  if( !lgOK )
409  {
410  fprintf(ioQQQ," big error was %e\n", x);
411  }
412  }
413 # endif
414 # endif
415 
416 
417  /*********************************************************
418  * *
419  * confirm that hydrogen einstein As are still valid *
420  * *
421  *********************************************************/
422  for( nelem=0; nelem<2; ++nelem )
423  {
424  /* this element may be turned off */
425  if( dense.lgElmtOn[nelem] )
426  {
427  /*Z4 = (double)(POW2(nelem+1)*POW2(nelem+1));*/
428  /* form charge to the 4th power */
429  Z4 = (double)(nelem+1);
430  Z4 *= Z4;
431  Z4 *= Z4;
432  /* H Lya */
433  ans1 = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul;
434  ans2 = 6.265e8*Z4;
435  if( fabs(ans1-ans2)/ans2 > 1e-3 )
436  {
437  fprintf(ioQQQ , "SanityCheck finds insane A for H Lya %g %g nelem=%li\n",
438  ans1 , ans2 , nelem );
439  lgOK = false;
440  }
441 
442 # if 0
443  /*must disable since, at this time, induced is included in Aul */
444  /* H two photon */
445  ans1 = Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Aul;
446  ans2 = 8.226*powi( (1.+nelem) , 6 );
447  if( fabs(ans1-ans2)/ans2 > 1e-3 )
448  {
449  fprintf(ioQQQ , "SanityCheck finds insane A for H 2-phot %g %g nelem=%li\n",
450  ans1 , ans2 , nelem );
451  lgOK = false;
452  }
453 
455  {
456  /* Balmer gamma 5 - 2 + 5-4 + 5-3*/
457  ans1 = Transitions[ipH_LIKE][nelem][ipH5p][ipH2s].Emis->Aul +
458  Transitions[ipH_LIKE][nelem][ipH5s][ipH2p].Emis->Aul +
459  Transitions[ipH_LIKE][nelem][ipH5d][ipH2p].Emis->Aul +
460  Transitions[ipH_LIKE][nelem][ipH5p][ipH3s].Emis->Aul +
461  Transitions[ipH_LIKE][nelem][ipH5s][ipH3p].Emis->Aul +
462  Transitions[ipH_LIKE][nelem][ipH5d][ipH3p].Emis->Aul +
463  Transitions[ipH_LIKE][nelem][ipH5p][ipH3d].Emis->Aul +
464  Transitions[ipH_LIKE][nelem][ipH5f][ipH3d].Emis->Aul +
465  Transitions[ipH_LIKE][nelem][ipH5p][ipH4s].Emis->Aul +
466  Transitions[ipH_LIKE][nelem][ipH5s][ipH4p].Emis->Aul +
467  Transitions[ipH_LIKE][nelem][ipH5d][ipH4p].Emis->Aul +
468  Transitions[ipH_LIKE][nelem][ipH5p][ipH4d].Emis->Aul +
469  Transitions[ipH_LIKE][nelem][ipH5f][ipH4d].Emis->Aul +
470  Transitions[ipH_LIKE][nelem][ipH5d][ipH4f].Emis->Aul +
471  Transitions[ipH_LIKE][nelem][ipH5g][ipH4f].Emis->Aul;
472  ans2 = (2.532e6+2.20e6+2.70e6)*Z4;
473  if( fabs(ans1-ans2)/ans2 > 1e-2 )
474  {
475  fprintf(ioQQQ ,
476  "SanityCheck finds insane A for H5-2 found=%g correct=%g nelem=%li\n",
477  ans1 , ans2 , nelem );
478  lgOK = false;
479  }
480  }
481 # endif
482  }
483  }
484 
485  /* check that hydrogenic branching ratios add up to unity */
486  for( nelem=0; nelem<LIMELM; ++nelem )
487  {
488  if( dense.lgElmtOn[nelem] )
489  {
490  int ipHi, ipLo;
491  for( ipHi=4; ipHi< iso.numLevels_max[ipH_LIKE][nelem]-iso.nCollapsed_max[ipH_LIKE][nelem]; ++ipHi )
492  {
493  double sum = 0.;
494  for( ipLo=0; ipLo<ipHi; ++ipLo )
495  {
496  sum += iso.BranchRatio[ipH_LIKE][nelem][ipHi][ipLo];
497  }
498  if( fabs(sum-1.)>0.01 )
499  {
500  fprintf(ioQQQ ,
501  "SanityCheck H branching ratio sum not unity for nelem=%li upper n=%i sum=%.3e\n",
502  nelem, ipHi, sum );
503  lgOK = false;
504  }
505  }
506  }
507  }
508 
509  fixit();/* make this work */
510 #if 0
511  /* check photo cross sections for H */
512  long ipISO = ipH_LIKE;
513  nelem = 0;
514  /* loop starts at 3, the first level with n = n and full l */
515  for( n=3; n < MIN2(100, iso.n_HighestResolved_max[ipISO][nelem]); ++n )
516  {
517  realnum anu[1]={1.} , cs[1];
518  double energy;
519  long index = iso.QuantumNumbers2Index[ipISO][nelem][n][0][2];
520 
521  /* photon energy where cross section will be evaluated,
522  * this is in Ryd */
523  energy = iso.xIsoLevNIonRyd[ipISO][nelem][index];
524  anu[0] = (realnum)(energy*1.01);
525 
526  H_photo_cs( photon , N_(n), L_(n), nelem+1 );
527  hypho(
528  /* Z-Nelec, the residual charge, 1 for hydrogen, 2 for helium */
529  (double)(nelem+1),
530  /* principal quantum number */
531  n,
532  /* lowest angular momentum */
533  0,
534  /* highest angular momentum */
535  n-1,
536  /* scaled lowest photon energy,
537  * => incorrect?? in units of zed^2/n^2,
538  * at which the cs will be done */
539  energy,
540  /* number of points to do */
541  1,
542  /* array of energies (in units given above) */
543  anu ,
544  /* calculated photoionization cross section in cm^-2 */
545  cs);
546 
547  error = fabs(cs[0] - opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/
548  ( (cs[0] + opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/2.);
549  /*fprintf(ioQQQ,"z=%ld n=%ld error %g old %e new %e\n",nelem,n, error,
550  opac.OpacStack[iso.ipOpac[ipISO][nelem][n]-1] ,cs[0] );*/
551  if( error > 0.05 )
552  {
553  fprintf(ioQQQ , "SanityCheck finds insane H photo cs\n");
554  lgOK = false;
555  }
556  }
557 #endif
558 
559  /*********************************************************
560  * *
561  * confirm that exponential integral routines still work *
562  * *
563  *********************************************************/
564 
565  /* check that first and second exponential integrals are ok,
566  * step through range of values, beginning with following */
567  x = 1e-3;
568  do
569  {
570  /* check that fast e1 routine is ok */
571  ans1 = ee1(x);
572  ans2 = expn( 1 , x );
573  if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
574  {
575  fprintf(ioQQQ , "SanityCheck finds insane E1 %g %g %g\n",
576  x , ans1 , ans2 );
577  lgOK = false;
578  }
579 
580  /* check that e2 is ok */
581  ans1 = e2(x);
582  ans2 = expn( 2 , x );
583  if( fabs(ans1-ans2)/(ans1+ans2) > 1e-6 )
584  {
585  fprintf(ioQQQ , "SanityCheck finds insane E2 %g %g %g\n",
586  x , ans1 , ans2 );
587  lgOK = false;
588  }
589 
590  /* now increment x */
591  x *= 2.;
592  /* following limit set by sexp returning zero, used in ee1 */
593  } while( x < 64. );
594 
595  /*********************************************************
596  * *
597  * confirm that matrix inversion routine still works *
598  * *
599  *********************************************************/
600 
601  /* these are the answer, chosen to get xvec 1,2,3 */
602  yVector[0] = 1.;
603  yVector[1] = 3.;
604  yVector[2] = 3.;
605 
606  /* zero out the main matrix */
607  for(i=0;i<3;++i)
608  {
609  for( j=0;j<3;++j )
610  {
611  xMatrix[i][j] = 0.;
612  }
613  }
614 
615  /* remember that order is column, row, alphabetical order, rc */
616  xMatrix[0][0] = 1.;
617  xMatrix[0][1] = 1.;
618  xMatrix[1][1] = 1.;
619  xMatrix[2][2] = 1.;
620 
621  /* this is the default matrix solver */
622  /* this test is the 1-d matrix with 2-d macro simulation */
623  /* LDA is right dimension of matrix */
624 
625  /* MALLOC space for the 1-d array */
626  A = (double*)MALLOC( sizeof(double)*NDIM*NDIM );
627 
628  /* copy over the main matrix */
629  for( i=0; i < 3; ++i )
630  {
631  for( j=0; j < 3; ++j )
632  {
633  A[i*NDIM+j] = xMatrix[i][j];
634  }
635  }
636 
637  ner = 0;
638 
639  /*void DGETRF(long,long,double*,long,long[],long*);*/
640  /*void DGETRF(int,int,double*,int,int[],int*);*/
641  getrf_wrapper(3, 3, A, NDIM, ipiv, &ner);
642  if( ner != 0 )
643  {
644  fprintf( ioQQQ, " SanityCheck DGETRF error\n" );
645  cdEXIT(EXIT_FAILURE);
646  }
647 
648  /* usage DGETRS, 'N' = no transpose
649  * order of matrix,
650  * number of cols in bvec, =1
651  * array
652  * leading dim of array */
653  /*void DGETRS(char,int,int,double*,int,int[],double*,int,int*);*/
654  getrs_wrapper('N', 3, 1, A, NDIM, ipiv, yVector, 3, &ner);
655 
656  if( ner != 0 )
657  {
658  fprintf( ioQQQ, " SanityCheck DGETRS error\n" );
659  cdEXIT(EXIT_FAILURE);
660  }
661  /* release the vector */
662  free( A );
663 
664  /* now check on validity of the solution, demand that this
665  * simple problem have come within a few epsilons of the
666  * correct answer */
667 
668  /* find largest deviation */
669  rcond = 0.;
670  for(i=0;i<3;++i)
671  {
672  x = fabs( yVector[i]-i-1.);
673  rcond = MAX2( rcond, x );
674  /*printf(" %g ", yVector[i]);*/
675  }
676 
677  if( rcond>DBL_EPSILON)
678  {
679  fprintf(ioQQQ,
680  "SanityCheck found too large a deviation in matrix solver = %g \n",
681  rcond);
682  /* set flag saying that things are not ok */
683  lgOK = false;
684  }
685  /* end matrix inversion check */
686 
687 
688  /* these pointers were set to INT_MIN in ContCreatePointers,
689  * then set to valid numbers in ipShells and OpacityCreate1Element
690  * this checks that all values have been properly filled */
691  for( nelem=0; nelem<LIMELM; ++nelem )
692  {
693  /* must reset state of code after tests performed, remember state here */
694  realnum xIonF[NISO][LIMELM];
695  double hbn[NISO][LIMELM], hn[NISO][LIMELM];
696 
697  if( dense.lgElmtOn[nelem] )
698  {
699  /* set these abundances so that opacities can be checked below */
700  hbn[ipH_LIKE][nelem] = iso.DepartCoef[ipH_LIKE][nelem][0];
701  hn[ipH_LIKE][nelem] = StatesElem[ipH_LIKE][nelem][0].Pop;
702  xIonF[ipH_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipH_LIKE];
703 
704  iso.DepartCoef[ipH_LIKE][nelem][0] = 0.;
705  StatesElem[ipH_LIKE][nelem][0].Pop = 1.;
706  dense.xIonDense[nelem][nelem+1-ipH_LIKE] = 1.;
707 
708  if( nelem > ipHYDROGEN )
709  {
710 
711  hbn[ipHE_LIKE][nelem] = iso.DepartCoef[ipHE_LIKE][nelem][0];
712  hn[ipHE_LIKE][nelem] = StatesElem[ipHE_LIKE][nelem][0].Pop;
713  xIonF[ipHE_LIKE][nelem] = dense.xIonDense[nelem][nelem+1-ipHE_LIKE];
714 
715  /* this does not exist for hydrogen itself */
716  iso.DepartCoef[ipHE_LIKE][nelem][0] = 0.;
717  StatesElem[ipHE_LIKE][nelem][0].Pop = 1.;
718  dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = 1.;
719  }
720 
721  for( ion=0; ion<=nelem; ++ion )
722  {
723  /* loop over all shells that are defined */
724  for( nshells=0; nshells<Heavy.nsShells[nelem][ion]; ++nshells )
725  {
726  for( j=0; j<3; ++j )
727  {
728  /* >>chng 00 apr 05, array index is on fortran scale so must be
729  * >= 1. This test had been <0, correct for C. Caught by Peter van Hoof */
730  if( opac.ipElement[nelem][ion][nshells][j] <=0 )
731  {
732  /* this is not possible */
733  fprintf(ioQQQ,
734  "SanityCheck found insane ipElement for nelem=%li ion=%li nshells=%li j=%li \n",
735  nelem , ion , nshells, j );
736  fprintf(ioQQQ,
737  "value was %li \n", opac.ipElement[nelem][ion][nshells][j] );
738  /* set flag saying that things are not ok */
739  lgOK = false;
740  }
741  }
742  }
743 
744  if( nelem > 1 )
745  {
746  realnum saveion[LIMELM+3];
747  /* check that photoionization cross sections are ok */
748  for( j=1; j <= (nelem + 2); j++ )
749  {
750  saveion[j] = dense.xIonDense[nelem][j-1];
751  dense.xIonDense[nelem][j-1] = 0.;
752  }
753 
754  dense.xIonDense[nelem][ion] = 1.;
755 
756  OpacityZero();
757  opac.lgRedoStatic = true;
758 
759  /* generate opacity with standard routine - this is the one
760  * called in OpacityAddTotal to make opacities in usual calculations */
761  OpacityAdd1Element(nelem);
762 
763  /* this starts one beyond energy of threshold since cs may be zero there */
764  for( j=Heavy.ipHeavy[nelem][ion]; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
765  {
766  if( opac.opacity_abs[j]+opac.OpacStatic[j] < FLT_MIN )
767  {
768  /* this is not possible */
769  fprintf(ioQQQ,
770  "SanityCheck found non-positive photo cs for nelem=%li ion=%li \n",
771  nelem , ion );
772  fprintf(ioQQQ,
773  "value was %.2e + %.2e nelem %li ion %li at energy %.2e\n",
774  opac.opacity_abs[j] ,
775  opac.OpacStatic[j] ,
776  nelem ,
777  ion ,
778  rfield.anu[j]);
779  /* set flag saying that things are not ok */
780  lgOK = false;
781  break;
782  }
783  }
784  /* reset the ionization distribution */
785  for( j=1; j <= (nelem + 2); j++ )
786  {
787  dense.xIonDense[nelem][j-1] = saveion[j];
788  }
789 
790  }
791  }
792  iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = hbn[ipH_LIKE][nelem];
793  StatesElem[ipH_LIKE][nelem][ipH1s].Pop = hn[ipH_LIKE][nelem];
794  dense.xIonDense[nelem][nelem+1-ipH_LIKE] = xIonF[ipH_LIKE][nelem];
795 
796  if( nelem > ipHYDROGEN )
797  {
798  iso.DepartCoef[ipHE_LIKE][nelem][ipHe1s1S] = hbn[ipHE_LIKE][nelem];
799  StatesElem[ipHE_LIKE][nelem][ipHe1s1S].Pop = hn[ipHE_LIKE][nelem];
800  dense.xIonDense[nelem][nelem+1-ipHE_LIKE] = xIonF[ipHE_LIKE][nelem];
801  }
802  }
803  }
804 
805 
806  /*********************************************************
807  * *
808  * everything is done, all checks make, did we pass them?*
809  * *
810  *********************************************************/
811 
812  if( lgOK )
813  {
814  /*return if ok */
815  if( trace.lgTrace )
816  {
817  fprintf( ioQQQ, " SanityCheck returns OK\n");
818  }
819  return;
820  }
821 
822  else
823  {
824  /* stop since problem encountered, lgEOF set false */
825  fprintf(ioQQQ , "SanityCheck finds insanity so exiting\n");
826  ShowMe();
827  cdEXIT(EXIT_FAILURE);
828  }
829 }

Generated for cloudy by doxygen 1.8.1.1