cloudy trunk
Loading...
Searching...
No Matches
atom_leveln.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/*atom_levelN compute an arbitrary N level atom */
4#include "cddefines.h"
5#include "physconst.h"
6#include "thermal.h"
7#include "conv.h"
8#include "phycon.h"
9#include "dense.h"
10#include "trace.h"
11#include "newton_step.h"
12#include "atoms.h"
13#include "dynamics.h"
14
16 /* nlev is the number of levels to compute*/
17 long int nLevelCalled,
18 /* ABUND is total abundance of species, used for nth equation
19 * if balance equations are homogeneous */
21 /* G(nlev) is stat weight of levels */
22 const double g[],
23 /* EX(nlev) is excitation potential of levels, deg K or wavenumbers
24 * 0 for lowest level, all are energy rel to ground NOT d(ENER) */
25 const double ex[],
26 /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */
27 char chExUnits,
28 /* populations [cm-3] of each level as deduced here */
29 double pops[],
30 /* departure coefficient, derived below */
31 double depart[],
32 /* net transition rate, A * esc prob, s-1 */
33 double ***AulEscp,
34 /* col str from up to low */
35 double ***col_str,
36 /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob,
37 * asserts confirm that [ihi][ilo] is zero */
38 double ***AulDest,
39 /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero */
40 double ***AulPump,
41 /* collision rates (s^-1), evaluated here and returned for cooling by calling function,
42 * unless following flag is true. If true then calling function has already filled
43 * in these rates. CollRate[i][j] is rate from i to j */
44 double ***CollRate,
45 /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */
46 const double source[],
47 /* this is an additional destruction rate to continuum, normally zero, units s-1 */
48 const double sink[],
49 /* flag saying whether CollRate already done, or we need to do it here,
50 * this is stored in data)[ihi][ilo] as either downward rate or collis strength*/
51 bool lgCollRateDone,
52 /* total cooling and its derivative, set here but nothing done with it*/
53 double *cooltl,
54 double *coolder,
55 /* string used to identify calling program in case of error */
56 const char *chLabel,
57 /* nNegPop flag indicating what we have done
58 * positive if negative populations occurred
59 * zero if normal calculation done
60 * negative if too cold, matrix not solved, since highest level had little excitation
61 * (for some atoms other routine will be called in this case) */
62 int *nNegPop,
63 /* true if populations are zero, either due to zero abundance of very low temperature */
64 bool *lgZeroPop,
65 /* option to print debug information */
66 bool lgDeBug,
67 /* option to do atom in LTE, can be omitted, default value false */
68 bool lgLTE,
69 /* array that will be set to the cooling per line, can be omitted, default value NULL */
71 /* array that will be set to the cooling derivative per line, can be omitted, default value NULL */
72 multi_arr<double,2> *dCooldT)
73{
74 long int level,
75 ihi,
76 ilo,
77 j;
78
79 int32 ner;
80
81 double cool1,
82 TeInverse,
83 TeConvFac;
84
85 DEBUG_ENTRY( "atom_levelN()" );
86
87 *nNegPop = -1;
88
89 /* >>chng 05 dec 14, units of ex[] can be Kelvin (old default) or wavenumbers */
90 if( chExUnits=='K' )
91 {
92 /* ex[] is in temperature units - this will multiply ex[] to
93 * obtain Boltzmann factor */
94 TeInverse = 1./phycon.te;
95 /* this multiplies ex[] to obtain energy difference between levels */
96 TeConvFac = 1.;
97 }
98 else if( chExUnits=='w' )
99 {
100 /* ex[] is in wavenumber units */
101 TeInverse = 1./phycon.te_wn;
102 TeConvFac = T1CM;
103 }
104 else
106
107 long int nlev = nLevelCalled;
108 // decrement number of levels until we have positive excitation rate,
109 while( ( (dsexp((ex[nlev-1]-ex[0])*TeInverse) + (*AulPump)[0][nlev-1]) < SMALLFLOAT ) &&
110 source[nlev-1]==0. && nlev > 1)
111 {
112 pops[nlev-1] = 0.;
113 depart[nlev-1] = 0.;
114 --nlev;
115 }
116
117 /* exit if zero abundance or all population in ground */
118 ASSERT( abund>= 0. );
119 if( abund == 0. || nlev==1 )
120 {
121 *cooltl = 0.;
122 *coolder = 0.;
123 /* says calc was ok */
124 *nNegPop = 0;
125 *lgZeroPop = true;
126
127 pops[0] = abund;
128 depart[0] = 1.;
129 for( level=1; level < nlev; level++ )
130 {
131 pops[level] = 0.;
132 depart[level] = 0.;
133 }
134 return;
135 }
136
137# ifndef NDEBUG
138 /* excitation temperature of lowest level must be zero */
139 ASSERT( ex[0] == 0. );
140
141 for( ihi=1; ihi < nlev; ihi++ )
142 {
143 for( ilo=0; ilo < ihi; ilo++ )
144 {
145 /* following must be zero:
146 * AulDest[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low
147 * AulEscp[ilo][ihi] - so that spontaneous transitions only proceed from high energy to low
148 * AulPump[ihi][ilo] - so that pumping only proceeds from low energy to high */
149 ASSERT( (*AulDest)[ilo][ihi] == 0. );
150 ASSERT( (*AulEscp)[ilo][ihi] == 0 );
151 ASSERT( (*AulPump)[ihi][ilo] == 0. );
152
153 ASSERT( (*AulEscp)[ihi][ilo] >= 0 );
154 ASSERT( (*AulDest)[ihi][ilo] >= 0 );
155 ASSERT( (*col_str)[ihi][ilo] >= 0 );
156 }
157 }
158# endif
159
160 if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
161 {
162 fprintf( ioQQQ, " atom_levelN trace printout for atom=%s with tot abund %e \n", chLabel, abund);
163 fprintf( ioQQQ, " AulDest\n" );
164
165 fprintf( ioQQQ, " hi lo" );
166
167 for( ilo=0; ilo < nlev-1; ilo++ )
168 {
169 fprintf( ioQQQ, "%4ld ", ilo );
170 }
171 fprintf( ioQQQ, " \n" );
172
173 for( ihi=1; ihi < nlev; ihi++ )
174 {
175 fprintf( ioQQQ, "%3ld", ihi );
176 for( ilo=0; ilo < ihi; ilo++ )
177 {
178 fprintf( ioQQQ, "%10.2e", (*AulDest)[ihi][ilo] );
179 }
180 fprintf( ioQQQ, "\n" );
181 }
182
183 fprintf( ioQQQ, " A*esc\n" );
184 fprintf( ioQQQ, " hi lo" );
185 for( ilo=0; ilo < nlev-1; ilo++ )
186 {
187 fprintf( ioQQQ, "%4ld ", ilo );
188 }
189 fprintf( ioQQQ, " \n" );
190
191 for( ihi=1; ihi < nlev; ihi++ )
192 {
193 fprintf( ioQQQ, "%3ld", ihi );
194 for( ilo=0; ilo < ihi; ilo++ )
195 {
196 fprintf( ioQQQ, "%10.2e", (*AulEscp)[ihi][ilo] );
197 }
198 fprintf( ioQQQ, "\n" );
199 }
200
201 fprintf( ioQQQ, " AulPump\n" );
202
203 fprintf( ioQQQ, " hi lo" );
204 for( ilo=0; ilo < nlev-1; ilo++ )
205 {
206 fprintf( ioQQQ, "%4ld ", ilo );
207 }
208 fprintf( ioQQQ, " \n" );
209
210 for( ihi=1; ihi < nlev; ihi++ )
211 {
212 fprintf( ioQQQ, "%3ld", ihi );
213 for( ilo=0; ilo < ihi; ilo++ )
214 {
215 fprintf( ioQQQ, "%10.2e", (*AulPump)[ilo][ihi] );
216 }
217 fprintf( ioQQQ, "\n" );
218 }
219
220 fprintf( ioQQQ, " coll str\n" );
221 fprintf( ioQQQ, " hi lo" );
222 for( ilo=0; ilo < nlev-1; ilo++ )
223 {
224 fprintf( ioQQQ, "%4ld ", ilo );
225 }
226 fprintf( ioQQQ, " \n" );
227
228 for( ihi=1; ihi < nlev; ihi++ )
229 {
230 fprintf( ioQQQ, "%3ld", ihi );
231 for( ilo=0; ilo < ihi; ilo++ )
232 {
233 fprintf( ioQQQ, "%10.2e", (*col_str)[ihi][ilo] );
234 }
235 fprintf( ioQQQ, "\n" );
236 }
237
238 fprintf( ioQQQ, " coll rate\n" );
239 fprintf( ioQQQ, " hi lo" );
240 for( ilo=0; ilo < nlev-1; ilo++ )
241 {
242 fprintf( ioQQQ, "%4ld ", ilo );
243 }
244 fprintf( ioQQQ, " \n" );
245
246 if( lgCollRateDone )
247 {
248 for( ihi=1; ihi < nlev; ihi++ )
249 {
250 fprintf( ioQQQ, "%3ld", ihi );
251 for( ilo=0; ilo < ihi; ilo++ )
252 {
253 fprintf( ioQQQ, "%10.2e", (*CollRate)[ihi][ilo] );
254 }
255 fprintf( ioQQQ, "\n" );
256 }
257 }
258 }
259
260 // these are all automatically deallocated when they go out of scope
261 multi_arr<double,2> excit(nLevelCalled,nLevelCalled);
262
263 /* find set of Boltzmann factors
264 * evaluate full range of levels since calling routine will use these values
265 * */
266 for( ilo=0; ilo < (nLevelCalled - 1); ilo++ )
267 {
268 for( ihi=ilo + 1; ihi < nLevelCalled; ihi++ )
269 {
270 /* >>chng 05 dec 14, option to have ex be either Kelvin or wavenumbers */
271 excit[ilo][ihi] = dsexp((ex[ihi]-ex[ilo])*TeInverse);
272 }
273 }
274
275 if( trace.lgTrace && trace.lgTrLevN )
276 {
277 fprintf( ioQQQ, " excit, te=%10.2e\n", phycon.te );
278 fprintf( ioQQQ, " hi lo" );
279
280 for( ilo=0; ilo < (nlev-1); ilo++ )
281 {
282 fprintf( ioQQQ, "%4ld ", ilo );
283 }
284 fprintf( ioQQQ, " \n" );
285
286 for( ihi=1; ihi < nlev; ihi++ )
287 {
288 fprintf( ioQQQ, "%3ld", ihi );
289 for( ilo=0; ilo < ihi; ilo++ )
290 {
291 fprintf( ioQQQ, "%10.2e", excit[ilo][ihi] );
292 }
293 fprintf( ioQQQ, "\n" );
294 }
295 }
296
297 /* we will predict populations */
298 *lgZeroPop = false;
299
300 /* already have excitation pumping, now get deexcitation */
301 for( ilo=0; ilo < (nlev - 1); ilo++ )
302 {
303 for( ihi=ilo + 1; ihi < nlev; ihi++ )
304 {
305 /* (*AulPump)[low][ihi] is excitation rate, low -> ihi, due to external continuum,
306 * so derive rate from upper to lower */
307 (*AulPump)[ihi][ilo] = (*AulPump)[ilo][ihi]*g[ilo]/g[ihi];
308 }
309 }
310
311 /* evaluate collision rates from collision strengths, but only if calling
312 * routine has not already done this
313 * evaluate full range of levels since calling routine will use these rates
314 */
315 if( !lgCollRateDone )
316 {
317 for( ilo=0; ilo < (nLevelCalled - 1); ilo++ )
318 {
319 for( ihi=ilo + 1; ihi < nLevelCalled; ihi++ )
320 {
321 /* this should be a collision strength */
322 ASSERT( (*col_str)[ihi][ilo]>= 0. );
323 /* this is deexcitation rate */
324 (*CollRate)[ihi][ilo] = (*col_str)[ihi][ilo]/g[ihi]*dense.cdsqte;
325 /* this is excitation rate */
326 (*CollRate)[ilo][ihi] = (*CollRate)[ihi][ilo]*g[ihi]/g[ilo]*
327 excit[ilo][ihi];
328 }
329 }
330 }
331
332 if( !lgLTE )
333 {
334 // these are all automatically deallocated when they go out of scope
335 valarray<double> bvec(nlev);
337
338 /* rate equations equal zero */
339 amat.zero();
340
341 /* eqns for destruction of level
342 * AulEscp[iho][ilo] are A*esc p, CollRate is coll excit in direction
343 * AulPump[low][high] is excitation rate from low to high */
344 for( level=0; level < nlev; level++ )
345 {
346 /* leaving level to below */
347 for( ilo=0; ilo < level; ilo++ )
348 {
349 double rate = (*CollRate)[level][ilo] + (*AulEscp)[level][ilo] +
350 (*AulDest)[level][ilo] + (*AulPump)[level][ilo];
351 amat[level][level] += rate;
352 amat[level][ilo] = -rate;
353 }
354
355 /* leaving level to above */
356 for( ihi=level + 1; ihi < nlev; ihi++ )
357 {
358 double rate = (*CollRate)[level][ihi] + (*AulPump)[level][ihi];
359 amat[level][level] += rate;
360 amat[level][ihi] = -rate;
361 }
362 }
363
364 if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax)
365 {
366 double slowrate = -1.0;
367 long slowlevel = -1;
368 // level == 0 is intentionally excluded...
369 for (level=1; level < nlev; ++level)
370 {
371 double rate = dynamics.timestep*amat[level][level];
372 if (rate < slowrate || slowrate < 0.0)
373 {
374 slowrate = rate;
375 slowlevel = level;
376 }
377 }
378 if ( slowrate < 3.0 )
379 {
380 fprintf(ioQQQ," PROBLEM in atom_levelN: "
381 "non-equilibrium appears important for %s(level=%ld), "
382 "rate * timestep is %g\n",
383 chLabel, slowlevel, slowrate);
384 }
385 }
386
387 double maxdiag = 0.;
388 for( level=0; level < nlev; level++ )
389 {
390 // source terms from elsewhere
391 bvec[level] = source[level];
392 if( amat[level][level] > maxdiag )
393 maxdiag = amat[level][level];
394 amat[level][level] += sink[level];
395 }
396
397 // If there are no sources or sinks, then one of the rows of the
398 // linear system is linearly dependent on the others so there is
399 // no matrix inverse. If the sources are sufficiently small,
400 // the answer will be numerically ill-conditioned.
401 //
402 // For strictly zero sources, this may be remedied by applying a
403 // conservation constraint instead of one of the redundant rows.
404 // To handle the broader case, we here add this conservation
405 // constraint scaled to switch in smoothly as the condition
406 // number of the matrix becomes poorer (and hence the accuracy
407 // of the linear system becomes poor).
408 //
409 // The condition number estimate here is quick but very crude.
410 //
411 // Need to specify smallest condition number before we decide
412 // that conservation will get a better answer than the raw
413 // linear solver. Numerical Recipes (3rd edition, S2.6.2)
414 // suggests that ~1e-14 might be appropriate for the solution of
415 // matrices in double precision.
416
417 const double CONDITION_SCALE = 1e-10;
418 double conserve_scale = maxdiag*CONDITION_SCALE;
419 ASSERT( conserve_scale > 0.0 );
420
421 for( level=0; level<nlev; ++level )
422 {
423 amat[level][0] += conserve_scale;
424 }
425 bvec[0] += abund*conserve_scale;
426
427 if( lgDeBug )
428 {
429 fprintf(ioQQQ," amat matrix follows:\n");
430 for( level=0; level < nlev; level++ )
431 {
432 for( j=0; j < nlev; j++ )
433 {
434 fprintf(ioQQQ," %.4e" , amat[level][j]);
435 }
436 fprintf(ioQQQ,"\n");
437 }
438 fprintf(ioQQQ," Vector follows:\n");
439 for( j=0; j < nlev; j++ )
440 {
441 fprintf(ioQQQ," %.4e" , bvec[j] );
442 }
443 fprintf(ioQQQ,"\n");
444 }
445
446 ner = solve_system(amat.vals(), bvec, nlev, NULL);
447
448 if( ner != 0 )
449 {
450 fprintf( ioQQQ, " atom_levelN: dgetrs finds singular or ill-conditioned matrix\n" );
452 }
453
454 /* set populations */
455 for( level=0; level < nlev; level++ )
456 {
457 /* save bvec into populations */
458 pops[level] = bvec[level];
459 }
460
461 /* now find total energy exchange rate, normally net cooling and its
462 * temperature derivative */
463 *cooltl = 0.;
464 *coolder = 0.;
465 for( ihi=1; ihi < nlev; ihi++ )
466 {
467 for( ilo=0; ilo < ihi; ilo++ )
468 {
469 /* this is now net cooling rate [erg cm-3 s-1] */
470 cool1 = (pops[ilo]*(*CollRate)[ilo][ihi] - pops[ihi]*(*CollRate)[ihi][ilo])*
471 (ex[ihi] - ex[ilo])*BOLTZMANN*TeConvFac;
472 *cooltl += cool1;
473 if( Cool != NULL )
474 (*Cool)[ihi][ilo] = cool1;
475 /* derivative wrt temperature - use Boltzmann factor relative to ground */
476 /* >>chng 03 aug 28, use real cool1 */
477 double dcooldT1 = cool1*( (ex[ihi] - ex[0])*thermal.tsq1 - thermal.halfte );
478 *coolder += dcooldT1;
479 if( dCooldT != NULL )
480 (*dCooldT)[ihi][ilo] = dcooldT1;
481 }
482 }
483
484 /* fill in departure coefficients */
485 if( pops[0] > SMALLFLOAT && excit[0][nlev-1] > SMALLFLOAT )
486 {
487 /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */
488 depart[0] = 1.;
489 for( ihi=1; ihi < nlev; ihi++ )
490 {
491 /* these are off by one - lowest index is zero */
492 depart[ihi] = (pops[ihi]/pops[0])*(g[0]/g[ihi])/excit[0][ihi];
493 }
494 }
495
496 else
497 {
498 /* >>chng 00 aug 10, loop had been from 1 and 0 was set to total abundance */
499 for( ihi=0; ihi < nlev; ihi++ )
500 {
501 /* Boltzmann factor or abundance too small, set departure coefficient to zero */
502 depart[ihi] = 0.;
503 }
504 depart[0] = 1.;
505 }
506 }
507 else
508 {
509 /* this branch calculates LTE level populations */
510
511 /* get the value of the partition function and the derivative */
512 valarray<double> help1(nlev), help2(nlev);
513 double partfun = help1[0] = g[0];
514 double dpfdT = help2[0] = 0.; /* help2 is d(help1)/dT */
515 for( level=1; level < nlev; level++ )
516 {
517 help1[level] = g[level]*excit[0][level];
518 partfun += help1[level];
519 help2[level] = help1[level]*(ex[level]-ex[0])*TeInverse/phycon.te;
520 dpfdT += help2[level];
521 }
522
523 /* calculate the level populations and departure coefficients,
524 * as well as the derivative of the level pops */
525 valarray<double> dndT(nlev);
526 for( level=0; level < nlev; level++ )
527 {
528 pops[level] = abund*help1[level]/partfun;
529 dndT[level] = abund*(help2[level]*partfun - dpfdT*help1[level])/pow2(partfun);
530 depart[level] = 1.;
531 }
532
533 /* now find the net cooling from the atom */
534 *cooltl = 0.;
535 *coolder = 0.;
536 for( ihi=1; ihi < nlev; ihi++ )
537 {
538 for( ilo=0; ilo < ihi; ilo++ )
539 {
540 double deltaE = (ex[ihi] - ex[ilo])*BOLTZMANN*TeConvFac;
541 double cool1 = (pops[ihi]*((*AulEscp)[ihi][ilo] + (*AulPump)[ihi][ilo]) -
542 pops[ilo]*(*AulPump)[ilo][ihi])*deltaE;
543 *cooltl += cool1;
544 if( Cool != NULL )
545 (*Cool)[ihi][ilo] = cool1;
546 double dcooldT1 = (dndT[ihi]*((*AulEscp)[ihi][ilo] + (*AulPump)[ihi][ilo]) -
547 dndT[ilo]*(*AulPump)[ilo][ihi])*deltaE;
548 *coolder += dcooldT1;
549 if( dCooldT != NULL )
550 (*dCooldT)[ihi][ilo] = dcooldT1;
551 }
552 }
553 }
554
555 /* sanity check for valid solution - non negative populations */
556 *nNegPop = 0;
557 /* the limit we allow the fractional population to go below zero before announcing failure. */
558 const double poplimit = 1.0e-10;
559 for( level=0; level < nlev; level++ )
560 {
561 if( pops[level] < 0. )
562 {
563 if( fabs(pops[level]/abund) > poplimit )
564 {
565 //nNegPop >= 1 leads to a failure
566 *nNegPop = *nNegPop + 1;
567 }
568 else
569 {
570 pops[level] = SMALLFLOAT;
571 //fprintf( ioQQQ, "\n PROBLEM Small negative populations were found in atom = %s . "
572 // "The problem was ignored and the negative populations were set to SMALLFLOAT",chLabel );
573 }
574 }
575 }
576
577 if( *nNegPop!=0 )
578 {
579 ASSERT( *nNegPop>=1 );
580 // *nNegPop 0 is valid solution, nNegPop> 1 negative populations found
581 fprintf( ioQQQ, "\n PROBLEM atom_levelN found negative population, nNegPop=%i, atom=%s lgSearch=%c\n",
582 *nNegPop , chLabel , TorF( conv.lgSearch ) );
583
584 for( level=0; level < nlev; level++ )
585 {
586 fprintf( ioQQQ, "%10.2e", pops[level] );
587 }
588
589 fprintf( ioQQQ, "\n" );
590 for( level=0; level < nlev; level++ )
591 {
592 pops[level] = (double)MAX2(0.,pops[level]);
593 }
594 }
595
596 if( lgDeBug || (trace.lgTrace && trace.lgTrLevN) )
597 {
598 fprintf( ioQQQ, "\n atom_leveln absolute population " );
599 for( level=0; level < nlev; level++ )
600 {
601 fprintf( ioQQQ, " %10.2e", pops[level] );
602 }
603 fprintf( ioQQQ, "\n" );
604
605 fprintf( ioQQQ, " departure coefficients" );
606 for( level=0; level < nlev; level++ )
607 {
608 fprintf( ioQQQ, " %10.2e", depart[level] );
609 }
610 fprintf( ioQQQ, "\n\n" );
611 }
612
613# ifndef NDEBUG
614 /* these were reset to non zero values by the solver, but we will
615 * assert that they are zero (for safety) when routine reenters so must
616 * set to zero here, but only if asserts are in place */
617 for( ihi=1; ihi < nlev; ihi++ )
618 {
619 for( ilo=0; ilo < ihi; ilo++ )
620 {
621 /* zero ots destruction rate */
622 (*AulDest)[ilo][ihi] = 0.;
623 /* both AulDest and AulPump (low, hi) are not used, should be zero */
624 (*AulPump)[ihi][ilo] = 0.;
625 (*AulEscp)[ilo][ihi] = 0;
626 }
627 }
628# endif
629
630 // -1 return had meant too cool and no populations determined, no longer used
631 // since we now decrement until populations can be determined
632 ASSERT( *nNegPop>=0 );
633
634 return;
635}
t_abund abund
Definition abund.cpp:5
static double * amat
void atom_levelN(long int nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double ***AulEscp, double ***col_str, double ***AulDest, double ***AulPump, double ***CollRate, const double source[], const double sink[], bool lgCollRateDone, double *cooltl, double *coolder, const char *chLabel, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE, multi_arr< double, 2 > *Cool, multi_arr< double, 2 > *dCooldT)
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
T pow2(T a)
Definition cddefines.h:931
#define EXIT_FAILURE
Definition cddefines.h:140
char TorF(bool l)
Definition cddefines.h:710
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
NORETURN void TotalInsanity(void)
Definition service.cpp:886
double dsexp(double x)
Definition service.cpp:953
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_conv conv
Definition conv.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_dynamics dynamics
Definition dynamics.cpp:44
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double T1CM
Definition physconst.h:167
static double ** AulEscp
Definition species2.cpp:29
static double * source
Definition species2.cpp:28
static double * ex
Definition species2.cpp:28
static double * pops
Definition species2.cpp:28
static double ** col_str
Definition species2.cpp:29
static double * depart
Definition species2.cpp:28
static double ** AulPump
Definition species2.cpp:29
static double * g
Definition species2.cpp:28
static double ** AulDest
Definition species2.cpp:29
static double * sink
Definition species2.cpp:28
static double ** CollRate
Definition species2.cpp:29
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5