cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atom_leveln.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*atom_levelN compute an arbitrary N level atom */
4 #include "cddefines.h"
5 #include "atoms.h"
6 #include "physconst.h"
7 #include "thermal.h"
8 #include "conv.h"
9 #include "phycon.h"
10 #include "dense.h"
11 #include "trace.h"
12 #include "newton_step.h"
13 #include "dynamics.h"
14 #include "vectorize.h"
15 #include "prt.h"
16 
17 // #define EIGEN
18 #ifdef EIGEN
19 extern "C" {
20  void dgeev_(const char*, const char*,int*,double*,int*,double*,
21  double*,
22  double*,int*,double*,int*,double*,int*,int*,int,int);
23 }
24 #endif
25 
26 // Set collision rate from collision strength
28  long nlev,
29  double TeInverse,
30  double **col_str,
31  const double ex[],
32  const double g[],
33  double **CollRate
34 )
35 {
36  resize(nlev);
37  for( long ihi=1; ihi < nlev; ++ihi )
38  {
39  double fac = dsexp((ex[ihi]-ex[ihi-1])*TeInverse);
40  for( long ilo=0; ilo < ihi-1; ++ilo )
41  {
42  excit[ihi][ilo] = fac*excit[ihi-1][ilo];
43  }
44  excit[ihi][ihi-1] = fac;
45  }
46 
47  if( trace.lgTrace && trace.lgTrLevN )
48  {
49  fprintf( ioQQQ, " coll str\n" );
50  fprintf( ioQQQ, " hi lo" );
51  for( long ilo=0; ilo < nlev-1; ilo++ )
52  {
53  fprintf( ioQQQ, "%4ld ", ilo );
54  }
55  fprintf( ioQQQ, " \n" );
56 
57  for( long ihi=1; ihi < nlev; ihi++ )
58  {
59  fprintf( ioQQQ, "%3ld", ihi );
60  for( long ilo=0; ilo < ihi; ilo++ )
61  {
62  fprintf( ioQQQ, "%10.2e", (col_str)[ihi][ilo] );
63  }
64  fprintf( ioQQQ, "\n" );
65  }
66 
67  fprintf( ioQQQ, " excit, te=%10.2e\n", phycon.te );
68  fprintf( ioQQQ, " hi lo" );
69 
70  for( long ilo=0; ilo < (nlev-1); ilo++ )
71  {
72  fprintf( ioQQQ, "%4ld ", ilo );
73  }
74  fprintf( ioQQQ, " \n" );
75 
76  for( long ihi=1; ihi < nlev; ihi++ )
77  {
78  fprintf( ioQQQ, "%3ld", ihi );
79  for( long ilo=0; ilo < ihi; ilo++ )
80  {
81  fprintf( ioQQQ, "%10.2e", excit[ihi][ilo] );
82  }
83  fprintf( ioQQQ, "\n" );
84  }
85  }
86 
87  for( long ihi=1; ihi < nlev; ihi++ )
88  {
89  for( long ilo=0; ilo < ihi; ilo++ )
90  {
91  /* this should be a collision strength */
92  ASSERT( (col_str)[ihi][ilo]>= 0. );
93  /* this is deexcitation rate */
94  (CollRate)[ihi][ilo] = (col_str)[ihi][ilo]/g[ihi]*dense.cdsqte;
95  /* this is excitation rate */
96  (CollRate)[ilo][ihi] = (CollRate)[ihi][ilo]*g[ihi]/g[ilo]*
97  excit[ihi][ilo];
98  }
99  }
100 }
101 
103  long int nLevelCalled,
104  /* ABUND is total abundance of species, used for nth equation
105  * if balance equations are homogeneous */
106  realnum abund,
107  /* G(nlev) is stat weight of levels */
108  const double g[],
109  /* EX(nlev) is excitation potential of levels, deg K or wavenumbers
110  * 0 for lowest level, all are energy rel to ground NOT d(ENER) */
111  const double ex[],
112  /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */
113  char chExUnits,
114  /* populations [cm-3] of each level as deduced here */
115  double pops[],
116  /* departure coefficient, derived below */
117  double depart[],
118  /* net transition rate, A * esc prob, s-1 */
119  double **AulEscp,
120  /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob,
121  * asserts confirm that [ihi][ilo] is zero */
122  double **AulDest,
123  /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero */
124  double **AulPump,
125  /* collision rates (s^-1), evaluated here and returned for cooling by calling function,
126  * unless following flag is true. If true then calling function has already filled
127  * in these rates. CollRate[i][j] is rate from i to j */
128  double **CollRate,
129  /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */
130  const double source[],
131  /* this is an additional destruction rate to continuum, normally zero, units s-1 */
132  const double sink[],
133  /* total cooling and its derivative, set here but nothing done with it*/
134  double *cooltl,
135  double *coolder,
136  /* string used to identify species in case of error */
137  const char *chLabel,
138  /* flag to print matrices input to solvers */
139  const bool lgPrtMatrix,
140  /* nNegPop flag indicating what we have done
141  * positive if negative populations occurred
142  * zero if normal calculation done
143  * negative if too cold, matrix not solved, since highest level had little excitation
144  * (for some atoms other routine will be called in this case) */
145  int *nNegPop,
146  /* true if populations are zero, either due to zero abundance of very low temperature */
147  bool *lgZeroPop,
148  /* option to print debug information */
149  bool lgDeBug,
150  /* option to do atom in LTE, can be omitted, default value false */
151  bool lgLTE,
152  /* array that will be set to the cooling per line, can be omitted, default value NULL */
153  multi_arr<double,2> *Cool,
154  /* array that will be set to the cooling derivative per line, can be omitted, default value NULL */
155  multi_arr<double,2> *dCooldT,
156  /* total excitation rate out of ground state */
157  double *grnd_excit)
158 {
159  long int level,
160  ihi,
161  ilo,
162  j;
163 
164  int32 ner;
165 
166  double cool1,
167  TeInverse,
168  TeConvFac;
169 
170  DEBUG_ENTRY( "atom_levelN()" );
171 
172  *nNegPop = -1;
173  if (grnd_excit != NULL) *grnd_excit = 0.0;
174 
175  /* >>chng 05 dec 14, units of ex[] can be Kelvin (old default) or wavenumbers */
176  if( chExUnits=='K' )
177  {
178  /* ex[] is in temperature units - this will multiply ex[] to
179  * obtain Boltzmann factor */
180  TeInverse = 1./phycon.te;
181  /* this multiplies ex[] to obtain energy difference between levels */
182  TeConvFac = 1.;
183  }
184  else if( chExUnits=='w' )
185  {
186  /* ex[] is in wavenumber units */
187  TeInverse = 1./phycon.te_wn;
188  TeConvFac = T1CM;
189  }
190  else
191  TotalInsanity();
192 
193  avx_ptr<double> arg(1,nLevelCalled), excit_gnd(1,nLevelCalled);
194  for( ihi=1; ihi < nLevelCalled; ++ihi )
195  arg[ihi] = -(ex[ihi]-ex[0])*TeInverse;
196  vexp( arg.ptr0(), excit_gnd.ptr0(), 1, nLevelCalled );
197 
198  long int nlev = nLevelCalled;
199  // decrement number of levels until we have positive excitation rate,
200  while( nlev > 1 && excit_gnd[nlev-1]+AulPump[0][nlev-1] < SMALLFLOAT &&
201  source[nlev-1] == 0. )
202  {
203  pops[nlev-1] = 0.;
204  depart[nlev-1] = 0.;
205  --nlev;
206  }
207 
208  /* exit if zero abundance or all population in ground */
209  ASSERT( abund>= 0. );
210  if( abund == 0. || nlev==1 )
211  {
212  *cooltl = 0.;
213  *coolder = 0.;
214  /* says calc was ok */
215  *nNegPop = 0;
216  *lgZeroPop = true;
217 
218  pops[0] = abund;
219  depart[0] = 1.;
220  for( level=1; level < nlev; level++ )
221  {
222  pops[level] = 0.;
223  depart[level] = 0.;
224  }
225  return;
226  }
227 
228 # ifndef NDEBUG
229  /* excitation temperature of lowest level must be zero */
230  ASSERT( ex[0] == 0. );
231 
232  for( ihi=0; ihi < nlev; ihi++ )
233  {
234  for( ilo=0; ilo < nlev; ilo++ )
235  {
236  /* AulDest[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low
237  * AulEscp[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low */
238  ASSERT( ihi == ilo || (AulDest)[ihi][ilo] >= 0. );
239  ASSERT( ihi == ilo || (AulEscp)[ihi][ilo] >= 0 );
240  }
241  }
242 
243  for( ihi=1; ihi < nlev; ihi++ )
244  {
245  for( ilo=0; ilo < ihi; ilo++ )
246  {
247  ASSERT( (AulPump)[ihi][ilo] >= 0. );
248  }
249  }
250 # endif
251 
252  if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
253  {
254  fprintf( ioQQQ, " atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel, abund);
255  fprintf( ioQQQ, " AulDest\n" );
256 
257  fprintf( ioQQQ, " hi lo" );
258 
259  for( ilo=0; ilo < nlev-1; ilo++ )
260  {
261  fprintf( ioQQQ, "%4ld ", ilo );
262  }
263  fprintf( ioQQQ, " \n" );
264 
265  for( ihi=1; ihi < nlev; ihi++ )
266  {
267  fprintf( ioQQQ, "%3ld", ihi );
268  for( ilo=0; ilo < ihi; ilo++ )
269  {
270  fprintf( ioQQQ, "%10.2e", (AulDest)[ihi][ilo] );
271  }
272  fprintf( ioQQQ, "\n" );
273  }
274 
275  fprintf( ioQQQ, " A*esc\n" );
276  fprintf( ioQQQ, " hi lo" );
277  for( ilo=0; ilo < nlev-1; ilo++ )
278  {
279  fprintf( ioQQQ, "%4ld ", ilo );
280  }
281  fprintf( ioQQQ, " \n" );
282 
283  for( ihi=1; ihi < nlev; ihi++ )
284  {
285  fprintf( ioQQQ, "%3ld", ihi );
286  for( ilo=0; ilo < ihi; ilo++ )
287  {
288  fprintf( ioQQQ, "%10.2e", (AulEscp)[ihi][ilo] );
289  }
290  fprintf( ioQQQ, "\n" );
291  }
292 
293  fprintf( ioQQQ, " AulPump\n" );
294 
295  fprintf( ioQQQ, " hi lo" );
296  for( ilo=0; ilo < nlev-1; ilo++ )
297  {
298  fprintf( ioQQQ, "%4ld ", ilo );
299  }
300  fprintf( ioQQQ, " \n" );
301 
302  for( ihi=1; ihi < nlev; ihi++ )
303  {
304  fprintf( ioQQQ, "%3ld", ihi );
305  for( ilo=0; ilo < ihi; ilo++ )
306  {
307  fprintf( ioQQQ, "%10.2e", (AulPump)[ilo][ihi] );
308  }
309  fprintf( ioQQQ, "\n" );
310  }
311 
312  fprintf( ioQQQ, " coll rate\n" );
313  fprintf( ioQQQ, " hi lo" );
314  for( ilo=0; ilo < nlev-1; ilo++ )
315  {
316  fprintf( ioQQQ, "%4ld ", ilo );
317  }
318  fprintf( ioQQQ, " \n" );
319 
320  for( ihi=1; ihi < nlev; ihi++ )
321  {
322  fprintf( ioQQQ, "%3ld", ihi );
323  for( ilo=0; ilo < ihi; ilo++ )
324  {
325  fprintf( ioQQQ, "%10.2e", (CollRate)[ihi][ilo] );
326  }
327  fprintf( ioQQQ, "\n" );
328  }
329  }
330 
331  /* we will predict populations */
332  *lgZeroPop = false;
333 
334  if( !lgLTE )
335  {
336  bvec.resize(nlev);
337  amat.alloc(nlev,nlev);
338 
339  /* rate equations equal zero */
340  amat.zero();
341 
342  /* eqns for destruction of level
343  * AulEscp[iho][ilo] are A*esc p, CollRate is coll excit in direction
344  * AulPump[low][high] is excitation rate from low to high */
345  for( level=0; level < nlev; level++ )
346  {
347  /* leaving level to below */
348  for( ilo=0; ilo < level; ilo++ )
349  {
350  double rate = (CollRate)[level][ilo] + (AulEscp)[level][ilo] +
351  (AulDest)[level][ilo] + (AulPump)[level][ilo];
352  amat[level][level] += rate;
353  amat[level][ilo] = -rate;
354  }
355 
356  /* leaving level to above */
357  for( ihi=level + 1; ihi < nlev; ihi++ )
358  {
359  double rate = (CollRate)[level][ihi] + (AulPump)[level][ihi];
360  amat[level][level] += rate;
361  amat[level][ihi] = -rate;
362  }
363  }
364  if (grnd_excit != NULL) *grnd_excit = amat[0][0];
365 
368  {
369  double slowrate = -1.0;
370  long slowlevel = -1;
371  // level == 0 is intentionally excluded...
372  for (level=1; level < nlev; ++level)
373  {
374  double rate = dynamics.timestep*amat[level][level];
375  if (rate < slowrate || slowrate < 0.0)
376  {
377  slowrate = rate;
378  slowlevel = level;
379  }
380  }
381  if ( slowrate < 3.0 )
382  {
383  static map<string,double> failures;
384  ostringstream oss;
385  oss << chLabel << "[level=" << slowlevel << ']';
386  string str=oss.str();
387  if (failures.find(str) == failures.end())
388  {
389  failures[str] = slowrate;
390  // >>chng 16 feb 18: demoted the following message from a PROBLEM
391  // to a CAUTION to prevent a flood of messages in minor.txt...
392  fprintf(ioQQQ," CAUTION in atom_levelN: "
393  "non-equilibrium appears important for %s, "
394  "rate * timestep is %g\n",
395  str.c_str(), slowrate);
396  }
397  else
398  {
399  if ( slowrate < failures[str])
400  failures[str] = slowrate;
401  }
402  }
403  }
404 
405  double maxdiag = 0.;
406  for( level=0; level < nlev; level++ )
407  {
408  // source terms from elsewhere
409  bvec[level] = source[level];
410  if( amat[level][level] > maxdiag )
411  maxdiag = amat[level][level];
412  amat[level][level] += sink[level];
413  }
414 
415  // If there are no sources or sinks, then one of the rows of the
416  // linear system is linearly dependent on the others so there is
417  // no matrix inverse. If the sources are sufficiently small,
418  // the answer will be numerically ill-conditioned.
419  //
420  // For strictly zero sources, this may be remedied by applying a
421  // conservation constraint instead of one of the redundant rows.
422  // To handle the broader case, we here add this conservation
423  // constraint scaled to switch in smoothly as the condition
424  // number of the matrix becomes poorer (and hence the accuracy
425  // of the linear system becomes poor).
426  //
427  // The condition number estimate here is quick but very crude.
428  //
429  // Need to specify smallest condition number before we decide
430  // that conservation will get a better answer than the raw
431  // linear solver. Numerical Recipes (3rd edition, S2.6.2)
432  // suggests that ~1e-14 might be appropriate for the solution of
433  // matrices in double precision.
434 
435  const double CONDITION_SCALE = 1e-10;
436  double conserve_scale = maxdiag*CONDITION_SCALE;
437  ASSERT( conserve_scale > 0.0 );
438 
439  for( level=0; level<nlev; ++level )
440  {
441  amat[level][0] += conserve_scale;
442  }
443  bvec[0] += abund*conserve_scale;
444 
445  if( lgDeBug )
446  {
447  fprintf(ioQQQ," amat matrix follows:\n");
448  for( level=0; level < nlev; level++ )
449  {
450  for( j=0; j < nlev; j++ )
451  {
452  fprintf(ioQQQ," %.4e" , amat[level][j]);
453  }
454  fprintf(ioQQQ,"\n");
455  }
456  fprintf(ioQQQ," Vector follows:\n");
457  for( j=0; j < nlev; j++ )
458  {
459  fprintf(ioQQQ," %.4e" , bvec[j] );
460  }
461  fprintf(ioQQQ,"\n");
462  }
463 
464 #ifdef EIGEN
465  if (lgDeBug)
466  {
467  multi_arr<double,2,C_TYPE> mcol(nlev,nlev),mrad(nlev,nlev),
468  vl(nlev,nlev),vr(nlev,nlev);
469  mcol.zero();
470  mrad.zero();
471 
472  for( level=0; level < nlev; level++ )
473  {
474  /* leaving level to below */
475  for( ilo=0; ilo < level; ilo++ )
476  {
477  double rate = (*CollRate)[level][ilo];
478  mcol[level][level] += rate;
479  mcol[level][ilo] = -rate;
480  rate = (*AulEscp)[level][ilo] +
481  (*AulDest)[level][ilo] + (*AulPump)[level][ilo];
482  mrad[level][level] += rate;
483  mrad[level][ilo] = -rate;
484  }
485  /* leaving level to above */
486  for( ihi=level + 1; ihi < nlev; ihi++ )
487  {
488  double rate = (*CollRate)[level][ihi];
489  mcol[level][level] += rate;
490  mcol[level][ihi] = -rate;
491  rate = (*AulPump)[level][ihi];
492  mrad[level][level] += rate;
493  mrad[level][ihi] = -rate;
494  }
495  }
496  multi_arr<double,2,C_TYPE> mtst(nlev,nlev);
497  for( ilo=0; ilo < nlev; ilo++ )
498  {
499  for( ihi=0; ihi < nlev; ++ihi)
500  {
501  mtst[ilo][ihi] = mcol[ilo][ihi];
502  }
503  }
504  const char job[2] = "V";
505  int lwork = 4*nlev,info,inlev=nlev;
506  vector<double> wr(nlev),wi(nlev),work(lwork);
507  dgeev_(job,job,&inlev,get_ptr(mtst.vals()),&inlev,get_ptr(wr),
508  get_ptr(wi),get_ptr(vl.vals()),&inlev,
509  get_ptr(vr.vals()),&inlev,get_ptr(work),&lwork,&info,1,1);
510  fprintf(ioQQQ,"info %d\nW=",info);
511  for( level=0; level < nlev; level++ )
512  {
513  fprintf(ioQQQ,"%ld %15.8g+%15.8gi\n",level,wr[level],wi[level]);
514  }
515 
516  // Reconstruct collision matrix
517  // mcol = u w v^t
518  // mcol[ihi][ilo] = sum(level,vl[level][ilo]*w[level]*vr[level][ihi])
519  // ( u^t mrad v + w ) v^t b = 0
520 
521  // Normalize left eigen-vectors so L & R are inverses
522  for( level=0; level < nlev; level++ )
523  {
524  double tot=0.0;
525  for( ilo=0; ilo < nlev; ilo++ )
526  {
527  tot += vl[level][ilo]*vr[level][ilo];
528  }
529  tot = 1.0/tot;
530  for( ilo=0; ilo < nlev; ilo++ )
531  {
532  vl[level][ilo] *= tot;
533  }
534  }
535 
536  mtst.zero();
537  for( ilo=0; ilo < nlev; ilo++ )
538  {
539  for( ihi=0; ihi < nlev; ihi++ )
540  {
541  for( level=0; level < nlev; level++ )
542  {
543  mtst[ilo][ihi] += mrad[level][ihi] * vr[ilo][level]; // or mcol to test...
544  }
545  }
546  }
547  multi_arr<double,2,C_TYPE> mtst1(nlev,nlev);
548  mtst1.zero();
549  for( ilo=0; ilo < nlev; ilo++ )
550  {
551  for( ihi=0; ihi < nlev; ihi++ )
552  {
553  for( level=0; level < nlev; level++ )
554  {
555  mtst1[ilo][ihi] += vl[ihi][level]*mtst[ilo][level];
556  }
557  }
558  }
559  for( level=0; level < nlev; level++ )
560  {
561  fprintf(ioQQQ,"%ld %15.8g %15.8g\n",level,wr[level],mtst1[level][level]);
562  }
563  }
564 #endif
565 
566  if( lgPrtMatrix )
567  {
568  prt.matrix.prtRates( nlev, amat, bvec );
569  }
570 
571  ner = solve_system(amat.vals(), bvec, nlev, NULL);
572 
573  if( ner != 0 )
574  {
575  fprintf( ioQQQ, " PROBLEM atom_levelN: dgetrs finds singular or ill-conditioned matrix for %s\n",
576  chLabel);
578  }
579 
580  /* set populations */
581  for( level=0; level < nlev; level++ )
582  {
583  /* save bvec into populations */
584  pops[level] = bvec[level];
585  }
586 
587  /* now find total energy exchange rate, normally net cooling and its
588  * temperature derivative */
589  *cooltl = 0.;
590  *coolder = 0.;
591  for( ihi=1; ihi < nlev; ihi++ )
592  {
593  for( ilo=0; ilo < ihi; ilo++ )
594  {
595  /* this is now net cooling rate [erg cm-3 s-1] */
596  cool1 = (pops[ilo]*(CollRate)[ilo][ihi] - pops[ihi]*(CollRate)[ihi][ilo])*
597  (ex[ihi] - ex[ilo])*BOLTZMANN*TeConvFac;
598  *cooltl += cool1;
599  if( Cool != NULL )
600  (*Cool)[ihi][ilo] = cool1;
601  /* derivative wrt temperature - use Boltzmann factor relative to ground */
602  /* >>chng 03 aug 28, use real cool1 */
603  double dcooldT1 = cool1*( (ex[ihi] - ex[0])*thermal.tsq1 - thermal.halfte );
604  *coolder += dcooldT1;
605  if( dCooldT != NULL )
606  (*dCooldT)[ihi][ilo] = dcooldT1;
607  }
608  }
609 
610  /* fill in departure coefficients */
611  if( pops[0] > SMALLFLOAT && excit_gnd[nlev-1] > SMALLFLOAT )
612  {
613  /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */
614  depart[0] = 1.;
615  for( ihi=1; ihi < nlev; ihi++ )
616  {
617  /* these are off by one - lowest index is zero */
618  depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/excit_gnd[ihi];
619  }
620  }
621 
622  else
623  {
624  /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */
625  for( ihi=0; ihi < nlev; ihi++ )
626  {
627  /* Boltzmann factor or abundance too small, set departure coefficient to zero */
628  depart[ihi] = 0.;
629  }
630  depart[0] = 1.;
631  }
632  }
633  else
634  {
635  /* this branch calculates LTE level populations */
636 
637  /* get the value of the partition function and the derivative */
638  valarray<double> help1(nlev), help2(nlev);
639  double partfun = help1[0] = g[0];
640  double dpfdT = help2[0] = 0.; /* help2 is d(help1)/dT */
641  for( level=1; level < nlev; level++ )
642  {
643  help1[level] = g[level]*excit_gnd[level];
644  partfun += help1[level];
645  help2[level] = help1[level]*(ex[level]-ex[0])*TeInverse/phycon.te;
646  dpfdT += help2[level];
647  }
648 
649  /* calculate the level populations and departure coefficients,
650  * as well as the derivative of the level pops */
651  valarray<double> dndT(nlev);
652  for( level=0; level < nlev; level++ )
653  {
654  pops[level] = abund*help1[level]/partfun;
655  dndT[level] = abund*(help2[level]*partfun - dpfdT*help1[level])/pow2(partfun);
656  depart[level] = 1.;
657  }
658 
659  /* now find the net cooling from the atom */
660  *cooltl = 0.;
661  *coolder = 0.;
662  for( ihi=1; ihi < nlev; ihi++ )
663  {
664  for( ilo=0; ilo < ihi; ilo++ )
665  {
666  double deltaE = (ex[ihi] - ex[ilo])*BOLTZMANN*TeConvFac;
667  double cool1 = (pops[ihi]*((AulEscp)[ihi][ilo] + (AulPump)[ihi][ilo]) -
668  pops[ilo]*(AulPump)[ilo][ihi])*deltaE;
669  *cooltl += cool1;
670  if( Cool != NULL )
671  (*Cool)[ihi][ilo] = cool1;
672  double dcooldT1 = (dndT[ihi]*((AulEscp)[ihi][ilo] + (AulPump)[ihi][ilo]) -
673  dndT[ilo]*(AulPump)[ilo][ihi])*deltaE;
674  *coolder += dcooldT1;
675  if( dCooldT != NULL )
676  (*dCooldT)[ihi][ilo] = dcooldT1;
677  }
678  }
679  }
680 
681  /* sanity check for valid solution - non negative populations */
682  *nNegPop = 0;
683  /* the limit we allow the fractional population to go below zero before announcing failure. */
684  const double poplimit = 1.0e-10;
685  for( level=0; level < nlev; level++ )
686  {
687  if( pops[level] < 0. )
688  {
689  if( fabs(pops[level]/abund) > poplimit )
690  {
691  //nNegPop >= 1 leads to a failure
692  *nNegPop = *nNegPop + 1;
693  }
694  else
695  {
696  pops[level] = SMALLFLOAT;
697  //fprintf( ioQQQ, "\n PROBLEM Small negative populations were found in atom = %s . "
698  // "The problem was ignored and the negative populations were set to SMALLFLOAT",chLabel );
699  }
700  }
701  }
702 
703  if( *nNegPop!=0 )
704  {
705  ASSERT( *nNegPop>=1 );
706  // *nNegPop 0 is valid solution, nNegPop> 1 negative populations found
707  fprintf( ioQQQ, "\n PROBLEM atom_levelN found negative population, nNegPop=%i, atom=%s lgSearch=%c\n",
708  *nNegPop , chLabel , TorF( conv.lgSearch ) );
709 
710  for( level=0; level < nlev; level++ )
711  {
712  fprintf( ioQQQ, "%10.2e", pops[level] );
713  }
714 
715  fprintf( ioQQQ, "\n" );
716  for( level=0; level < nlev; level++ )
717  {
718  pops[level] = (double)MAX2(0.,pops[level]);
719  }
720  }
721 
722  if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
723  {
724  fprintf( ioQQQ, "\n atom_leveln absolute population " );
725  for( level=0; level < nlev; level++ )
726  {
727  fprintf( ioQQQ, " %10.2e", pops[level] );
728  }
729  fprintf( ioQQQ, "\n" );
730 
731  fprintf( ioQQQ, " departure coefficients" );
732  for( level=0; level < nlev; level++ )
733  {
734  fprintf( ioQQQ, " %10.2e", depart[level] );
735  }
736  fprintf( ioQQQ, "\n\n" );
737  }
738 
739 # ifndef NDEBUG
740  /* these were reset to non zero values by the solver, but we will
741  * assert that they are zero (for safety) when routine reenters so must
742  * set to zero here, but only if asserts are in place */
743  for( ilo=0; ilo < nlev; ilo++ )
744  {
745  for( ihi=ilo+1; ihi < nlev; ihi++ )
746  {
747  /* zero ots destruction rate */
748  AulDest[ilo][ihi] = 0.;
749  }
750  }
751  for( ihi=1; ihi < nlev; ihi++ )
752  {
753  for( ilo=0; ilo < ihi; ilo++ )
754  {
755  /* both AulDest and AulPump (low, hi) are not used, should be zero */
756  AulPump[ihi][ilo] = 0.;
757  }
758  }
759  for( ilo=0; ilo < nlev; ilo++ )
760  {
761  for( ihi=ilo+1; ihi < nlev; ihi++ )
762  {
763  AulEscp[ilo][ihi] = 0;
764  }
765  }
766 # endif
767 
768  // -1 return had meant too cool and no populations determined, no longer used
769  // since we now decrement until populations can be determined
770  ASSERT( *nNegPop>=0 );
771 
772  return;
773 }
t_thermal thermal
Definition: thermal.cpp:6
T * ptr0()
Definition: vectorize.h:221
T * get_ptr(T *v)
Definition: cddefines.h:1113
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double depart(const genericState &gs)
void resize(long int nlev)
Definition: atoms.h:52
const realnum SMALLFLOAT
Definition: cpu.h:246
void prtRates(const long nlevels_local, const multi_arr< double, 2, C_TYPE > &a, valarray< double > &b)
Definition: prt.cpp:232
double cdsqte
Definition: dense.h:246
char TorF(bool l)
Definition: cddefines.h:753
bool lgAdvection
Definition: dynamics.h:66
bool lgTimeDependentStatic
Definition: dynamics.h:102
t_conv conv
Definition: conv.cpp:5
t_phycon phycon
Definition: phycon.cpp:6
FILE * ioQQQ
Definition: cddefines.cpp:7
t_dynamics dynamics
Definition: dynamics.cpp:42
double dsexp(double x)
Definition: service.cpp:1134
t_dense dense
Definition: global.cpp:15
void vexp(const double x[], double y[], long nlo, long nhi)
long int iteration
Definition: cddefines.cpp:16
bool lgSearch
Definition: conv.h:168
t_trace trace
Definition: trace.cpp:5
t_abund abund
Definition: abund.cpp:5
bool lgTrace
Definition: trace.h:12
double tsq1
Definition: thermal.h:142
double timestep
Definition: dynamics.h:187
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
void operator()(long nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double **AulEscp, double **AulDest, double **AulPump, double **CollRate, const double create[], const double destroy[], double *cooltl, double *coolder, const char *chLabel, bool lgPrtMatrix, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE=false, multi_arr< double, 2 > *Cool=NULL, multi_arr< double, 2 > *dCooldT=NULL, double *grnd_excit=NULL)
#define cdEXIT(FAIL)
Definition: cddefines.h:484
double halfte
Definition: thermal.h:142
t_prt prt
Definition: prt.cpp:14
vector< double * > excit
Definition: atoms.h:48
valarray< T > & vals()
#define ASSERT(exp)
Definition: cddefines.h:617
multi_arr< double, 2, C_TYPE > amat
Definition: atoms.h:82
T pow2(T a)
Definition: cddefines.h:985
t_prt_matrix matrix
Definition: prt.h:238
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
long int n_initial_relax
Definition: dynamics.h:132
double te_wn
Definition: phycon.h:30
double te
Definition: phycon.h:21
bool lgTrLevN
Definition: trace.h:31
void operator()(long nlev, double TeInverse, double **col_str, const double ex[], const double g[], double **CollRate)
Definition: atom_leveln.cpp:27
valarray< double > bvec
Definition: atoms.h:81