cloudy  trunk
 All Data Structures 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-2008 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 "physconst.h"
6 #include "thermal.h"
7 #include "phycon.h"
8 #include "dense.h"
9 #include "trace.h"
10 #include "thirdparty.h"
11 #include "atoms.h"
12 
14  /* nlev is the number of levels to compute*/
15  long int nlev,
16  /* ABUND is total abundance of species, used for nth equation
17  * if balance equations are homogeneous */
18  realnum abund,
19  /* G(nlev) is stat weight of levels */
20  const double g[],
21  /* EX(nlev) is excitation potential of levels, deg K or wavenumbers
22  * 0 for lowest level, all are energy rel to ground NOT d(ENER) */
23  const double ex[],
24  /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */
25  char chExUnits,
26  /* populations [cm-3] of each level as deduced here */
27  double pops[],
28  /* departure coefficient, derived below */
29  double depart[],
30  /* net transition rate, A * esc prob, s-1 */
31  double ***AulEscp,
32  /* col str from up to low */
33  double ***col_str,
34  /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob,
35  * asserts confirm that [ihi][ilo] is zero */
36  double ***AulDest,
37  /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero */
38  double ***AulPump,
39  /* collision rates (s^-1), evaluated here and returned for cooling by calling function,
40  * unless following flag is true. If true then calling function has already filled
41  * in these rates. CollRate[i][j] is rate from i to j */
42  double ***CollRate,
43  /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */
44  const double source[] ,
45  /* this is an additional destruction rate to continuum, normally zero, units s-1 */
46  const double sink[] ,
47  /* flag saying whether CollRate already done, or we need to do it here,
48  * this is stored in data)[ihi][ilo] as either downward rate or collis strength*/
49  bool lgCollRateDone,
50  /* total cooling and its derivative, set here but nothing done with it*/
51  double *cooltl,
52  double *coolder,
53  /* string used to identify calling program in case of error */
54  const char *chLabel,
55  /* nNegPop flag indicating what we have done
56  * positive if negative populations occurred
57  * zero if normal calculation done
58  * negative if too cold (for some atoms other routine will be called in this case) */
59  int *nNegPop,
60  /* true if populations are zero, either due to zero abundance of very low temperature */
61  bool *lgZeroPop ,
62  /* option to print debug information */
63  bool lgDeBug )
64 {
65  bool lgHomogeneous;
66 
67  long int level,
68  ihi,
69  ilo,
70  j;
71 
72  int32 ner;
73 
74  double cool1,
75  TeInverse,
76  TeConvFac,
77  sum;
78 
79  DEBUG_ENTRY( "atom_levelN()" );
80 
81  /* check for zero abundance and exit if so */
82  if( abund <= 0. )
83  {
84  *cooltl = 0.;
85  *coolder = 0.;
86  /* says calc was ok */
87  *nNegPop = false;
88  *lgZeroPop = true;
89 
90  for( level=0; level < nlev; level++ )
91  {
92  pops[level] = 0.;
93  depart[level] = 0.;
94  }
95 
96  depart[0] = 1.;
97 
98  /* there are TWO abort returns in this sub,
99  * this one is for zero abundance */
100  return;
101  }
102 
103  // these are all automatically deallocated when they go out of scope
104  auto_vec<int32> ipiv( new int32[nlev] );
105  auto_vec<double> bvec( new double[nlev] );
106  multi_arr<double,2,C_TYPE> amat(nlev,nlev);
107  multi_arr<double,2> excit(nlev,nlev);
108 
109 # ifndef NDEBUG
110  /* excitation temperature of lowest level must be zero */
111  ASSERT( ex[0] == 0. );
112 
113  for( ihi=1; ihi < nlev; ihi++ )
114  {
115  for( ilo=0; ilo < ihi; ilo++ )
116  {
117  /* following must be zero:
118  * AulDest[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low
119  * AulEscp[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low
120  * AulPump[ihi][ilo] - so that pumping only proceeds from low energy to high */
121  ASSERT( (*AulDest)[ilo][ihi] == 0. );
122  ASSERT( (*AulEscp)[ilo][ihi] == 0 );
123  ASSERT( (*AulPump)[ihi][ilo] == 0. );
124 
125  ASSERT( (*AulEscp)[ihi][ilo] >= 0 );
126  ASSERT( (*AulDest)[ihi][ilo] >= 0 );
127  ASSERT( (*col_str)[ihi][ilo] >= 0 );
128  }
129  }
130 # endif
131 
132  if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
133  {
134  fprintf( ioQQQ, " atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel, abund);
135  fprintf( ioQQQ, " AulDest\n" );
136 
137  fprintf( ioQQQ, " hi lo" );
138 
139  for( ilo=0; ilo < nlev-1; ilo++ )
140  {
141  fprintf( ioQQQ, "%4ld ", ilo );
142  }
143  fprintf( ioQQQ, " \n" );
144 
145  for( ihi=1; ihi < nlev; ihi++ )
146  {
147  fprintf( ioQQQ, "%3ld", ihi );
148  for( ilo=0; ilo < ihi; ilo++ )
149  {
150  fprintf( ioQQQ, "%10.2e", (*AulDest)[ihi][ilo] );
151  }
152  fprintf( ioQQQ, "\n" );
153  }
154 
155  fprintf( ioQQQ, " A*esc\n" );
156  fprintf( ioQQQ, " hi lo" );
157  for( ilo=0; ilo < nlev-1; ilo++ )
158  {
159  fprintf( ioQQQ, "%4ld ", ilo );
160  }
161  fprintf( ioQQQ, " \n" );
162 
163  for( ihi=1; ihi < nlev; ihi++ )
164  {
165  fprintf( ioQQQ, "%3ld", ihi );
166  for( ilo=0; ilo < ihi; ilo++ )
167  {
168  fprintf( ioQQQ, "%10.2e", (*AulEscp)[ihi][ilo] );
169  }
170  fprintf( ioQQQ, "\n" );
171  }
172 
173  fprintf( ioQQQ, " AulPump\n" );
174 
175  fprintf( ioQQQ, " hi lo" );
176  for( ilo=0; ilo < nlev-1; ilo++ )
177  {
178  fprintf( ioQQQ, "%4ld ", ilo );
179  }
180  fprintf( ioQQQ, " \n" );
181 
182  for( ihi=1; ihi < nlev; ihi++ )
183  {
184  fprintf( ioQQQ, "%3ld", ihi );
185  for( ilo=0; ilo < ihi; ilo++ )
186  {
187  fprintf( ioQQQ, "%10.2e", (*AulPump)[ilo][ihi] );
188  }
189  fprintf( ioQQQ, "\n" );
190  }
191 
192  fprintf( ioQQQ, " coll str\n" );
193  fprintf( ioQQQ, " hi lo" );
194  for( ilo=0; ilo < nlev-1; ilo++ )
195  {
196  fprintf( ioQQQ, "%4ld ", ilo );
197  }
198  fprintf( ioQQQ, " \n" );
199 
200  for( ihi=1; ihi < nlev; ihi++ )
201  {
202  fprintf( ioQQQ, "%3ld", ihi );
203  for( ilo=0; ilo < ihi; ilo++ )
204  {
205  fprintf( ioQQQ, "%10.2e", (*col_str)[ihi][ilo] );
206  }
207  fprintf( ioQQQ, "\n" );
208  }
209 
210  fprintf( ioQQQ, " coll rate\n" );
211  fprintf( ioQQQ, " hi lo" );
212  for( ilo=0; ilo < nlev-1; ilo++ )
213  {
214  fprintf( ioQQQ, "%4ld ", ilo );
215  }
216  fprintf( ioQQQ, " \n" );
217 
218  if( lgCollRateDone )
219  {
220  for( ihi=1; ihi < nlev; ihi++ )
221  {
222  fprintf( ioQQQ, "%3ld", ihi );
223  for( ilo=0; ilo < ihi; ilo++ )
224  {
225  fprintf( ioQQQ, "%10.2e", (*CollRate)[ihi][ilo] );
226  }
227  fprintf( ioQQQ, "\n" );
228  }
229  }
230  }
231 
232  /* >>chng 05 dec 14, units of ex[] can be Kelvin (old default) or wavenumbers */
233  if( chExUnits=='K' )
234  {
235  /* ex[] is in temperature units - this will multiply ex[] to
236  * obtain Boltzmann factor */
237  TeInverse = 1./phycon.te;
238  /* this multiplies ex[] to obtain energy difference between levels */
239  TeConvFac = 1.;
240  }
241  else if( chExUnits=='w' )
242  {
243  /* ex[] is in wavenumber units */
244  TeInverse = 1./phycon.te_wn;
245  TeConvFac = T1CM;
246  }
247  else
248  TotalInsanity();
249 
250  /* find set of Boltzmann factors */
251  for( ilo=0; ilo < (nlev - 1); ilo++ )
252  {
253  for( ihi=ilo + 1; ihi < nlev; ihi++ )
254  {
255  /* >>chng 05 dec 14, option to have ex be either Kelvin or wavenumbers */
256  excit[ilo][ihi] = sexp((ex[ihi]-ex[ilo])*TeInverse);
257  }
258  }
259 
260  if( trace.lgTrace && trace.lgTrLevN )
261  {
262  fprintf( ioQQQ, " excit, te=%10.2e\n", phycon.te );
263  fprintf( ioQQQ, " hi lo" );
264 
265  for( ilo=0; ilo < (nlev-1); ilo++ )
266  {
267  fprintf( ioQQQ, "%4ld ", ilo );
268  }
269  fprintf( ioQQQ, " \n" );
270 
271  for( ihi=1; ihi < nlev; ihi++ )
272  {
273  fprintf( ioQQQ, "%3ld", ihi );
274  for( ilo=0; ilo < ihi; ilo++ )
275  {
276  fprintf( ioQQQ, "%10.2e", excit[ilo][ihi] );
277  }
278  fprintf( ioQQQ, "\n" );
279  }
280  }
281 
282  /* punt if total excitation rate is zero */
283  /* >>chng 04 sep 16, add test on non-zero source */
284  if( (excit[0][nlev-1] + (*AulPump)[0][nlev-1] < 1e-13 ) && (source[nlev-1]==0.) )
285  {
286  *cooltl = 0.;
287  *coolder = 0.;
288  /* special flag saying too cool for highest level to be computed.
289  * some routines will call another routine for lower levels in this case */
290  *nNegPop = -1;
291  *lgZeroPop = true;
292 
293  for( level=1; level < nlev; level++ )
294  {
295  pops[level] = 0.;
296  /* these are off by one - lowest index is zero */
297  depart[level] = 0.;
298  }
299 
300  /* everything in ground */
301  pops[0] = abund;
302  depart[0] = 1.;
303 
304  /* there are two error exits from this routine,
305  * previous one for zero abundance, and this one for zero excitation */
306  return;
307  }
308 
309  /* we will predict populations */
310  *lgZeroPop = false;
311 
312  /* already have excitation pumping, now get deexcitation */
313  for( ilo=0; ilo < (nlev - 1); ilo++ )
314  {
315  for( ihi=ilo + 1; ihi < nlev; ihi++ )
316  {
317  /* (*AulPump)[low][ihi] is excitation rate, low -> ihi, due to external continuum,
318  * so derive rate from upper to lower */
319  (*AulPump)[ihi][ilo] = (*AulPump)[ilo][ihi]*g[ilo]/g[ihi];
320  }
321  }
322 
323  /* evaluate collision rates from collision strengths, but only if calling
324  * routine has not already done this */
325  if( !lgCollRateDone )
326  {
327  for( ilo=0; ilo < (nlev - 1); ilo++ )
328  {
329  for( ihi=ilo + 1; ihi < nlev; ihi++ )
330  {
331  /* this should be a collision strength */
332  ASSERT( (*col_str)[ihi][ilo]>= 0. );
333  /* this is deexcitation rate */
334  (*CollRate)[ihi][ilo] = (*col_str)[ihi][ilo]/g[ihi]*dense.cdsqte;
335  /* this is excitation rate */
336  (*CollRate)[ilo][ihi] = (*CollRate)[ihi][ilo]*g[ihi]/g[ilo]*
337  excit[ilo][ihi];
338  }
339  }
340  }
341 
342  /* rate equations equal zero */
343  amat.zero();
344 
345  /* following is column of vector - represents source terms from elsewhere,
346  * if this is zero then matrix is singular and must replace one row with
347  * population sum equation - if sum is non-zero then get total abundance
348  * from source and sink terms */
349  sum = 0.;
350  lgHomogeneous = false;
351  for( level=0; level < nlev; level++ )
352  {
353  bvec[level] = source[level];
354  sum += bvec[level];
355  }
356  if( sum==0. )
357  lgHomogeneous = true;
358 
359  /* eqns for destruction of level
360  * AulEscp[iho][ilo] are A*esc p, CollRate is coll excit in direction
361  * AulPump[low][high] is excitation rate from low to high */
362  for( level=0; level < nlev; level++ )
363  {
364  amat[level][level] = sink[level];
365 
366  /* leaving level to below */
367  for( ilo=0; ilo < level; ilo++ )
368  {
369  amat[level][level] += (*CollRate)[level][ilo] + (*AulEscp)[level][ilo] +
370  (*AulDest)[level][ilo] + (*AulPump)[level][ilo];
371  /* >>chng 97 jul 31, added pumping down
372  * coming to level I from below */
373  amat[ilo][level] = -(*CollRate)[ilo][level] - (*AulPump)[ilo][level];
374  }
375 
376  /* leaving level to above */
377  for( ihi=level + 1; ihi < nlev; ihi++ )
378  {
379  amat[level][level] += (*CollRate)[level][ihi] + (*AulPump)[level][ihi];
380  /* coming to level from above */
381  amat[ihi][level] = -(*CollRate)[ihi][level] - (*AulEscp)[ihi][level] -
382  (*AulDest)[ihi][level] - (*AulPump)[ihi][level];
383  /* >>chng 97 jul 31, added pumping down */
384  }
385  }
386 
387  /* homogeneous case, all source terms add up to zero, so use the population,
388  * system has no total population information,
389  * equation to replace redundant equation */
390  if( lgHomogeneous )
391  {
392  for( level=0; level<nlev; ++level )
393  {
394  amat[level][0] = 1.0;
395  }
396  /* these add up to total abundance */
397  bvec[0] = abund;
398  }
399 
400  if( lgDeBug )
401  {
402  fprintf(ioQQQ," amat matrix follows:\n");
403  for( level=0; level < nlev; level++ )
404  {
405  for( j=0; j < nlev; j++ )
406  {
407  fprintf(ioQQQ," %.4e" , amat[level][j]);
408  }
409  fprintf(ioQQQ,"\n");
410  }
411  if( sum==0. )
412  {
413  fprintf(ioQQQ," Sum creation zero so pop sum used\n");
414  }
415  else
416  {
417  fprintf(ioQQQ," Sum creation non-zero (%e), vector follows:\n",sum);
418  for( j=0; j < nlev; j++ )
419  {
420  fprintf(ioQQQ," %.4e" , bvec[j] );
421  }
422  fprintf(ioQQQ,"\n");
423  }
424  }
425 
426  ner = 0;
427  getrf_wrapper(nlev,nlev,amat.data(),nlev,ipiv.data(),&ner);
428  /* usage DGETRS, 'N' = no transpose
429  * order of matrix,
430  * number of cols in bvec, =1
431  * array
432  * leading dim of array */
433  getrs_wrapper('N',nlev,1,amat.data(),nlev,ipiv.data(),bvec.data(),nlev,&ner);
434 
435  if( ner != 0 )
436  {
437  fprintf( ioQQQ, " atom_levelN: dgetrs finds singular or ill-conditioned matrix\n" );
438  cdEXIT(EXIT_FAILURE);
439  }
440 
441  /* set populations */
442  for( level=0; level < nlev; level++ )
443  {
444  /* save bvec into populations */
445  pops[level] = bvec[level];
446  }
447 
448  /* now find total energy exchange rate, normally net cooling and its
449  * temperature derivative */
450  *cooltl = 0.;
451  *coolder = 0.;
452  for( ihi=1; ihi < nlev; ihi++ )
453  {
454  for( ilo=0; ilo < ihi; ilo++ )
455  {
456  /* this is now net cooling rate [K cm-3 s-1] */
457  cool1 = (pops[ilo]*(*CollRate)[ilo][ihi] - pops[ihi]*(*CollRate)[ihi][ilo])*
458  (ex[ihi] - ex[ilo]);
459  *cooltl += cool1;
460 
461  /* derivative wrt temperature - use Boltzmann factor relative to ground */
462  /* >>chng 03 aug 28, use real cool1 */
463  *coolder += cool1*( (ex[ihi] - ex[0])*thermal.tsq1 - thermal.halfte);
464  }
465  }
466  /* convert from units of ex[] into ergs */
467  /* >>chng 05 dec 14, ex[] may be K or wn, TeConvFac will take care of either case */
468  *cooltl *= BOLTZMANN*TeConvFac;
469  *coolder *= BOLTZMANN*TeConvFac;
470 
471  /* fill in departure coefficients */
472  if( pops[0] > 1e-20 && excit[0][nlev-1] > 1e-20 )
473  {
474  /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */
475  depart[0] = 1.;
476  for( ihi=1; ihi < nlev; ihi++ )
477  {
478  /* these are off by one - lowest index is zero */
479  depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/excit[0][ihi];
480  }
481  }
482 
483  else
484  {
485  /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */
486  for( ihi=0; ihi < nlev; ihi++ )
487  {
488  /* Boltzmann factor or abundance too small, set departure coefficient to zero */
489  depart[ihi] = 0.;
490  }
491  depart[0] = 1.;
492  }
493 
494  /* sanity check for valid solution - non negative populations */
495  *nNegPop = false;
496  for( level=0; level < nlev; level++ )
497  {
498  if( pops[level] < 0. )
499  {
500  *nNegPop = true;
501  }
502  }
503 
504  if( *nNegPop )
505  {
506  fprintf( ioQQQ, "\n PROBLEM atom_levelN found negative population, atom=%s\n",
507  chLabel );
508  fprintf( ioQQQ, " absolute population=" );
509 
510  for( level=0; level < nlev; level++ )
511  {
512  fprintf( ioQQQ, "%10.2e", pops[level] );
513  }
514 
515  fprintf( ioQQQ, "\n" );
516  for( level=0; level < nlev; level++ )
517  {
518  pops[level] = (double)MAX2(0.,pops[level]);
519  }
520  }
521 
522  if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
523  {
524  fprintf( ioQQQ, "\n atom_leveln absolute population " );
525  for( level=0; level < nlev; level++ )
526  {
527  fprintf( ioQQQ, " %10.2e", pops[level] );
528  }
529  fprintf( ioQQQ, "\n" );
530 
531  fprintf( ioQQQ, " departure coefficients" );
532  for( level=0; level < nlev; level++ )
533  {
534  fprintf( ioQQQ, " %10.2e", depart[level] );
535  }
536  fprintf( ioQQQ, "\n\n" );
537  }
538 
539 # ifndef NDEBUG
540  /* these were reset to non zero values by the solver, but we will
541  * assert that they are zero (for safety) when routine reenters so must
542  * set to zero here, but only if asserts are in place */
543  for( ihi=1; ihi < nlev; ihi++ )
544  {
545  for( ilo=0; ilo < ihi; ilo++ )
546  {
547  /* zero ots destruction rate */
548  (*AulDest)[ilo][ihi] = 0.;
549  /* both AulDest and AulPump (low, hi) are not used, should be zero */
550  (*AulPump)[ihi][ilo] = 0.;
551  (*AulEscp)[ilo][ihi] = 0;
552  }
553  }
554 # endif
555  return;
556 }

Generated for cloudy by doxygen 1.8.1.2