cloudy trunk
Loading...
Searching...
No Matches
iso_level.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/*iso_level solve for iso-sequence level populations */
4#include "cddefines.h"
5#include "atmdat.h"
6#include "continuum.h"
7#include "conv.h"
8#include "dense.h"
9#include "dynamics.h"
10#include "elementnames.h"
11#include "grainvar.h"
12#include "he.h"
13#include "heavy.h"
14#include "helike.h"
15#include "hmi.h"
16#include "hydrogenic.h"
17#include "ionbal.h"
18#include "iso.h"
19#include "mole.h"
20#include "opacity.h"
21#include "phycon.h"
22#include "physconst.h"
23#include "rfield.h"
24#include "secondaries.h"
25#include "taulines.h"
26#include "thirdparty.h"
27#include "trace.h"
28#include "yield.h"
29
30/* solve for level populations */
31void iso_level( const long int ipISO, const long int nelem, double &renorm)
32{
33 long int ipHi,
34 ipLo,
35 i,
36 level,
37 level_error;
38
39 t_iso_sp* sp = &iso_sp[ipISO][nelem];
40 const long int numlevels_local = sp->numLevels_local;
41
42 double BigError;
43
44 int32 nerror;
45 double HighestPColOut = 0.;
46 bool lgNegPop=false;
47 valarray<int32> ipiv(numlevels_local);
48 /* this block of variables will be obtained and freed here */
49 double source=0.,
50 sink=0.;
51 valarray<double> PopPerN(sp->n_HighestResolved_local+1);
52
53 DEBUG_ENTRY( "iso_level()" );
54 renorm = 1.;
55
56 /* check that we were called with valid charge */
57 ASSERT( nelem >= ipISO );
58 ASSERT( nelem < LIMELM );
59
60 /* these two collision rates must be the same or we are in big trouble,
61 * since used interchangeably */
62 ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
63 fabs( (sp->fb[0].ColIoniz* dense.EdenHCorr) /
64 SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );
65
66 if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
67 (dense.IonLow[nelem] <= nelem - ipISO) &&
68 ionbal.RateRecomIso[nelem][ipISO] > 0. )
69 {
70 /* get simple estimate of atom to ion ratio */
71 sp->xIonSimple = ionbal.RateIonizTot(nelem,nelem-ipISO)/ionbal.RateRecomIso[nelem][ipISO];
72 }
73 else
74 {
75 sp->xIonSimple = 0.;
76 }
77
78 /* which case atom to solve??? */
79 if( dense.xIonDense[nelem][nelem+1-ipISO] < SMALLFLOAT || sp->xIonSimple < 1e-35 )
80 {
81 /* don't bother if no ionizing radiation */
82 strcpy( sp->chTypeAtomUsed, "zero " );
83 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
84 {
85 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n",
86 ipISO, nelem, sp->xIonSimple, sp->chTypeAtomUsed );
87 }
88
89 /* total ionization will just be the ground state */
90 lgNegPop = false;
91 sp->st[0].Pop() = dense.xIonDense[nelem][nelem-ipISO];
92 for( long n=1; n < numlevels_local; n++ )
93 {
94 sp->st[n].Pop() = 0.;
95 }
96 sp->qTot2S = 0.;
97 }
98 else
99 {
100 /* fill in recombination vector - values were set in iso_ionize_recombine.cpp */
101 valarray<double>
102 creation(numlevels_local),
103 error(numlevels_local),
104 work(numlevels_local),
105 Save_creation(numlevels_local);
106
107 ASSERT( dense.xIonDense[nelem][nelem+1-ipISO] >= 0.f );
108 for( level=0; level < numlevels_local; level++ )
109 {
110 /* total recombination from once more ionized [cm-3 s-1] */
111 creation[level] = sp->fb[level].RateCont2Level * dense.xIonDense[nelem][nelem+1-ipISO];
112 }
113
114 double ctsource=0.0, ctsink=0.0, ctrec=0.0;
115 /* now charge transfer - all into/from ground, two cases, H and not H */
116 if( nelem==ipHYDROGEN )
117 {
118 /* charge transfer, hydrogen onto everything else */
119 /* charge exchange recombination */
120 ctrec += atmdat.CharExcRecTotal[nelem];
121 ctsink += atmdat.CharExcIonTotal[nelem];
122 }
123 else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
124 {
125 /* this is recom of He due to ct with all other gas constituents */
126 ctrec += atmdat.CharExcRecTotal[nelem];
127 ctsink += atmdat.CharExcIonTotal[nelem];
128 }
129 else
130 {
131 for (long nelem1=ipHYDROGEN; nelem1 < t_atmdat::NCX; ++nelem1)
132 {
133 long ipISO=nelem1;
134 ctrec += atmdat.CharExcRecTo[nelem1][nelem][nelem-ipISO]*iso_sp[ipISO][nelem1].st[0].Pop();
135 ctsink += atmdat.CharExcIonOf[nelem1][nelem][nelem-ipISO]*dense.xIonDense[nelem1][1];
136 }
137 }
138 ctsource += ctrec*dense.xIonDense[nelem][nelem+1-ipISO];
139
140 if ( nelem > ipISO )
141 {
142 double ction=0.0;
143 if( nelem==ipHELIUM )
144 {
145 ctsink += atmdat.CharExcRecTotal[nelem];
146 ction += atmdat.CharExcIonTotal[nelem];
147 }
148 else
149 {
150 for (long nelem1=ipHYDROGEN; nelem1 < t_atmdat::NCX; ++nelem1)
151 {
152 long ipISO1=nelem1;
153 ctsink += atmdat.CharExcRecTo[nelem1][nelem][nelem-(ipISO+1)]*iso_sp[ipISO1][nelem1].st[0].Pop();
154 ction += atmdat.CharExcIonOf[nelem1][nelem][nelem-(ipISO+1)]*dense.xIonDense[nelem1][1];
155 }
156 }
157 ctsource += ction * dense.xIonDense[nelem][nelem-(ipISO+1)];
158 }
159
160 /* now do the 2D array */
162 z.alloc(numlevels_local,numlevels_local);
163 z.zero();
164 vector<double> Boltzmann(numlevels_local-1);
165 for (unsigned lev = 0; lev < Boltzmann.size(); ++lev)
166 Boltzmann[lev] = 1.0;
167
168 /* this branch is main solver, full level populations
169 * assert since this code must change if NISO ever increased */
170 ASSERT( NISO == 2 );
171 long SpinUsed[NISO] = {2,1};
172 long indexNmaxP =
173 sp->QuantumNumbers2Index[ sp->n_HighestResolved_local ][1][SpinUsed[ipISO]];
174
175 strcpy( sp->chTypeAtomUsed, "Popul" );
176 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
177 {
178 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n",
179 ipISO, nelem, sp->chTypeAtomUsed );
180 }
181
182 //qList::const_iterator StElm = StatesElemNEW[nelem][nelem-ipISO].begin();
183 qList::const_iterator StElm = sp->st.begin();
184
185 /* master balance equation, use when significant population */
186 for( level=0; level < numlevels_local; level++ )
187 {
188 /* all process depopulating level and placing into the continuum
189 * this does NOT include grain charge transfer ionization, added below */
190 z[level][level] += sp->fb[level].RateLevel2Cont;
191
192 if( level == 1 )
193 /* >>chng 05 dec 21, rm eden to make into rate coefficient */
194 sp->qTot2S = sp->fb[level].ColIoniz;
195
196 //TransitionList::iterator Trans = sp->tr.begin() + sp->ipTrans[level][0];
197
198 if (level != 0)
199 {
200 double arg = (StElm[level].energy().K()-StElm[level-1].energy().K()) /
201 phycon.te;
202 arg = MAX2( -38., arg); fixit(); // Wouldn't need to mask this out if levels were in order
203 double bstep = sexp( arg );
204 for ( ipLo = 0; ipLo < level; ++ipLo )
205 Boltzmann[ipLo] *= bstep;
206 }
207
208 /* all processes populating level from below */
209 for( ipLo=0; ipLo < level; ipLo++ )
210 {
211 double coll_down = sp->trans(level,ipLo).Coll().ColUL( colliders );
212
213 double RadDecay = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().Aul()*
214 (sp->trans(level,ipLo).Emis().Pesc() +
215 sp->trans(level,ipLo).Emis().Pelec_esc() +
216 sp->trans(level,ipLo).Emis().Pdest())*
217 KILL_BELOW_PLASMA(sp->trans(level,ipLo).EnergyRyd()) );
218
219 double pump = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().pump() *
220 KILL_BELOW_PLASMA(sp->trans(level,ipLo).EnergyRyd()) );
221
222 if( iso_ctrl.lgRandErrGen[ipISO] )
223 {
224 coll_down *= sp->ex[level][ipLo].ErrorFactor[IPCOLLIS];
225 RadDecay *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
226 pump *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
227 }
228
229 double coll_up = coll_down *
230 (double)StElm[level].g()/
231 (double)StElm[ipLo].g()*
232 Boltzmann[ipLo];
233
234 z[ipLo][ipLo] += coll_up + pump ;
235 z[ipLo][level] -= coll_up + pump ;
236
237 double pump_down = pump *
238 (double)StElm[ipLo].g()/
239 (double)StElm[level].g();
240
241 z[level][level] += RadDecay + pump_down + coll_down;
242 z[level][ipLo] -= RadDecay + pump_down + coll_down;
243
244 if( level == indexNmaxP )
245 {
246 HighestPColOut += coll_down;
247 }
248 if( ipLo == indexNmaxP )
249 {
250 HighestPColOut += coll_up;
251 }
252
253 /* find total collisions out of 2^3S to singlets. */
254 if( (level == 1) && (ipLo==0) )
255 {
256 sp->qTot2S += coll_down/dense.EdenHCorr;
257 }
258 if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S()==1) )
259 {
260 sp->qTot2S += coll_up/dense.EdenHCorr;
261 }
262 }
263 }
264
265 if( ipISO == nelem )
266 {
267 /* iso_ctrl.lgCritDensLMix[ipISO] is a flag used to print warning if density is
268 * too low for first collapsed level to be l-mixed. Check is if l-mixing collisions
269 * out of highest resolved singlet P are greater than sum of transition probs out. */
270 if( HighestPColOut < 1./sp->st[indexNmaxP].lifetime() )
271 {
272 iso_ctrl.lgCritDensLMix[ipISO] = false;
273 }
274 }
275
277 ASSERT( ipISO <= ipHE_LIKE );
278 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
279 {
280 // induced two photon emission - special because upward and downward are
281 // not related by ratio of statistical weights
282 // iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command
283
284 fixit(); // need Pesc or the equivalent to multiply AulTotal?
285 // downward rate
286 z[tnu->ipHi][tnu->ipLo] -= tnu->AulTotal;
287 z[tnu->ipHi][tnu->ipHi] += tnu->AulTotal;
288
289 }
290
291 if (iso_ctrl.lgInd2nu_On)
292 {
293 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
294 {
295 // induced two photon emission - special because upward and downward are
296 // not related by ratio of statistical weights
297 // iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command
298
299 z[tnu->ipHi][tnu->ipLo] -= tnu->induc_dn;
300 z[tnu->ipHi][tnu->ipHi] += tnu->induc_dn;
301
302 // upward rate
303 z[tnu->ipLo][tnu->ipHi] -= tnu->induc_up;
304 z[tnu->ipLo][tnu->ipLo] += tnu->induc_up;
305 }
306 }
307
308 /* grain charge transfer recombination and ionization to ALL other stages */
309 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
310 {
311 if( ion!=nelem-ipISO )
312 {
313 source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] *
314 dense.xIonDense[nelem][ion];
315 sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
316
317 source += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] *
318 dense.xIonDense[nelem][ion] * atmdat.lgCTOn;
319 sink += mole.xMoleChTrRate[nelem][nelem-ipISO][ion] * atmdat.lgCTOn;
320
321 }
322 }
323
324 /* add in source and sink terms from molecular network. */
325 source += mole.source[nelem][nelem-ipISO];
326 sink += mole.sink[nelem][nelem-ipISO];
327
328#if 1
329 /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
330 /*>>>chng 02 oct 01, only include if lgAdvection is set */
331 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
332 dynamics.Rate != 0.0 &&
333 !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
334 {
335 /* add in advection - these terms normally zero */
336 source += dynamics.Source[nelem][nelem-ipISO];
337 /* >>chng 02 Sep 06 rjrw -- advective term not recombination */
338 sink += dynamics.Rate;
339 }
340#else
341 /*>>>chng 02 oct 01, only include if lgAdvection is set */
342 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
343 dynamics.Rate != 0.0 &&
344 !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
345 {
346 for( level=0; level < numlevels_local; level++ )
347 {
348 creation[level] += dynamics.StatesElem[nelem][nelem-ipISO][level];
349 z[level][level] += dynamics.Rate;
350 }
351 }
352#endif
353
354 /* ionization from/recombination from lower ionization stages */
355 for(long ion_from=dense.IonLow[nelem]; ion_from < MIN2( dense.IonHigh[nelem], nelem-ipISO ) ; ion_from++ )
356 {
357 if( ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] >= 0. )
358 {
359 /* ionization from lower ionization stages, cm-3 s-1 */
360 source += ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] * dense.xIonDense[nelem][ion_from];
361 }
362 /* recombination to next lower ionization stage, s-1 */
363 if( ion_from == nelem-1-ipISO )
364 sink += ionbal.RateRecomTot[nelem][ion_from];
365 }
366
367 ASSERT( source >= 0.f );
368 if (0)
369 {
370 creation[0] += source;
371 for( level=0; level < numlevels_local; level++ )
372 {
373 z[level][level] += sink;
374 }
375 }
376 else
377 {
378 // Try Boltzmann weighting to capture LTE limit correctly
379 t_iso_sp* sp = &iso_sp[ipISO][nelem];
380 double partfun=0.0;
381 for ( level = 0; level < numlevels_local; level++ )
382 {
383 partfun += sp->st[level].Boltzmann()*sp->st[level].g();
384 }
385 source /= partfun;
386 for( level=0; level < numlevels_local; level++ )
387 {
388 creation[level] += source*
389 sp->st[level].Boltzmann()*sp->st[level].g();
390 z[level][level] += sink;
391 }
392 }
393 creation[0] += ctsource;
394 z[0][0] += ctsink;
395
396 /* >>chng 04 nov 30, atom XX-like collisions off turns this off too */
397 if( sp->trans(iso_ctrl.nLyaLevel[ipISO],0).Coll().rate_lu_nontherm() * iso_ctrl.lgColl_excite[ipISO] > 0. )
398 {
399 /* now add on supra thermal excitation */
400 for( level=1; level < numlevels_local; level++ )
401 {
402 double RateUp , RateDown;
403
404 RateUp = sp->trans(level,0).Coll().rate_lu_nontherm();
405 RateDown = RateUp * (double)sp->st[ipH1s].g() /
406 (double)sp->st[level].g();
407
408 /* total rate out of lower level */
409 z[ipH1s][ipH1s] += RateUp;
410
411 /* rate from the upper level to ground */
412 z[level][ipH1s] -= RateDown;
413
414 /* rate from ground to upper level */
415 z[ipH1s][level] -= RateUp;
416
417 z[level][level] += RateDown;
418 }
419 }
420
421 /* ===================================================================
422 *
423 * at this point all matrix elements have been established
424 *
425 * ==================================================================== */
426 /* save matrix, this allocates SaveZ */
427 SaveZ = z;
428
429 for( ipLo=0; ipLo < numlevels_local; ipLo++ )
430 Save_creation[ipLo] = creation[ipLo];
431
432 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
433 {
434 const long numlevels_print = numlevels_local;
435 fprintf( ioQQQ, " pop level others => (iso_level)\n" );
436 for( long n=0; n < numlevels_print; n++ )
437 {
438 fprintf( ioQQQ, " %s %s %2ld", iso_ctrl.chISO[ipISO], elementnames.chElementNameShort[nelem], n );
439 for( long j=0; j < numlevels_print; j++ )
440 {
441 fprintf( ioQQQ,"\t%.9e", z[j][n] );
442 }
443 fprintf( ioQQQ, "\n" );
444 }
445 fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
446 atmdat.CharExcRecTotal[ipHELIUM],
447 findspecieslocal("CO")->den ,
448 atmdat.CharExcRecTo[ipHELIUM][nelem][nelem-ipISO]*iso_sp[ipHE_LIKE][ipHELIUM].st[0].Pop() ,
449 atmdat.CharExcRecTo[ipHYDROGEN][nelem][nelem-ipISO]*iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop() );
450 fprintf(ioQQQ," recomb ");
451 for( long n=0; n < numlevels_print; n++ )
452 {
453 fprintf( ioQQQ,"\t%.9e", creation[n] );
454 }
455 fprintf( ioQQQ, "\n" );
456 }
457
458 nerror = 0;
459
460 getrf_wrapper(numlevels_local,numlevels_local,
461 z.data(),numlevels_local,&ipiv[0],&nerror);
462
463 getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0],
464 &creation[0],numlevels_local,&nerror);
465
466 if( nerror != 0 )
467 {
468 fprintf( ioQQQ, " iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
470 }
471
472 /* check whether solution is valid */
473 /* >>chng 06 aug 28, both of these from numLevels_max to _local. */
474 for( level=ipH1s; level < numlevels_local; level++ )
475 {
476 double qn = 0., qx = 0.;
477 error[level] = 0.;
478 for( long n=ipH1s; n < numlevels_local; n++ )
479 {
480 double q = SaveZ[n][level]*creation[n];
481
482 /* remember the largest size of element in sum to div by below */
483 if ( q > qx )
484 qx = q;
485 else if (q < qn)
486 qn = q;
487
488 error[level] += q;
489 }
490
491 if (-qn > qx)
492 qx = -qn;
493
494 if( qx > 0. )
495 {
496 error[level] = (error[level] - Save_creation[level])/qx;
497 }
498 else
499 {
500 error[level] = 0.;
501 }
502 }
503
504 /* remember largest residual in matrix inversion */
505 BigError = -1.;
506 level_error = -1;
507 /* >>chng 06 aug 28, from numLevels_max to _local. */
508 for( level=ipH1s; level < numlevels_local; level++ )
509 {
510 double abserror;
511 abserror = fabs( error[level]);
512 /* this will be the largest residual in the matrix inversion */
513 if( abserror > BigError )
514 {
515 BigError = abserror;
516 level_error = level;
517 }
518 }
519
520 /* matrix inversion should be nearly as good as the accuracy of a double,
521 * but demand that it is better than epsilon for a float */
522 if( BigError > FLT_EPSILON )
523 {
524 if( !conv.lgSearch )
525 fprintf(ioQQQ," PROBLEM" );
526
527 fprintf(ioQQQ,
528 " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
529 "level was %li Search?%c \n",
530 fnzone,
531 ipISO,
532 elementnames.chElementName[nelem],
533 nelem ,
534 BigError ,
535 level_error,
536 TorF(conv.lgSearch) );
537 }
538
539 // Force level balance to LTE
540 if ( iso_ctrl.lgLTE_levels[ipISO] )
541 {
542 t_iso_sp* sp = &iso_sp[ipISO][nelem];
543 double partfun=0.0;
544 for ( level = 0; level < numlevels_local; level++ )
545 {
546 partfun += sp->st[level].Boltzmann()*sp->st[level].g();
547 }
548 double scale = dense.xIonDense[nelem][nelem-ipISO]/partfun;
549 for ( level = 0; level < numlevels_local; level++ )
550 {
551 creation[level] = sp->st[level].Boltzmann()*sp->st[level].g()*scale;
552 }
553 }
554
555 for( level = numlevels_local-1; level > 0; --level )
556 {
557 /* check for negative populations */
558 if( creation[level] < 0. )
559 lgNegPop = true;
560 }
561
562 if( lgNegPop && dense.lgSetIoniz[nelem] )
563 {
564 // simulation can become unphysical if ionization is fixed.
565 // in this case, just put everything in ground.
566 // It's really the best we can do.
567 for( level = 1; level < numlevels_local; ++level )
568 creation[level] = 0.;
569 creation[0] = dense.xIonDense[nelem][nelem-ipISO];
570 // flip flag back
571 lgNegPop = false;
572 }
573
574 /* put level populations into master array */
575 for( level=0; level < numlevels_local; level++ )
576 {
577 ASSERT( creation[level] >= 0. );
578 sp->st[level].Pop() = creation[level];
579
580 if( sp->st[level].Pop() <= 0 && !conv.lgSearch )
581 {
582 fprintf(ioQQQ,
583 "PROBLEM non-positive level pop for iso = %li, nelem = "
584 "%li = %s, level=%li val=%.3e nTotalIoniz %li\n",
585 ipISO,
586 nelem ,
587 elementnames.chElementSym[nelem],
588 level,
589 sp->st[level].Pop() ,
590 conv.nTotalIoniz);
591 }
592 }
593
594 /* zero populations of unused levels. */
595 for( level=numlevels_local; level < sp->numLevels_max; level++ )
596 {
597 sp->st[level].Pop() = 0.;
598 /* >>chng 06 jul 25, no need to zero this out, fix limit to 3-body heating elsewhere. */
599 /* sp->st[level].PopLTE = 0.; */
600 }
601
602 /* TotalPopExcited is sum of excited level pops */
603 /* renormalize the populations to agree with ion solver */
604 iso_renorm( nelem, ipISO, renorm );
605
606 double TotalPopExcited = 0.;
607 /* create sum of populations */
608 for( level=1; level < numlevels_local; level++ )
609 TotalPopExcited += sp->st[level].Pop();
610 ASSERT( TotalPopExcited >= 0. );
611 double TotalPop = TotalPopExcited + sp->st[0].Pop();
612
613 /* option to force ionization */
614 if( dense.lgSetIoniz[nelem] )
615 {
616 if( !fp_equal( TotalPop, (double)dense.xIonDense[nelem][nelem-ipISO] ) )
617 {
618 if( TotalPopExcited >= dense.xIonDense[nelem][nelem-ipISO] )
619 {
620 for( level=0; level < numlevels_local; level++ )
621 sp->st[level].Pop() *=
622 dense.xIonDense[nelem][nelem-ipISO] / TotalPop;
623 }
624 else
625 {
626 sp->st[0].Pop() =
627 MAX2( 1e-30 * dense.xIonDense[nelem][nelem-ipISO],
628 dense.xIonDense[nelem][nelem-ipISO] - TotalPopExcited );
629 }
630 sp->lgPopsRescaled = true;
631 }
632 ASSERT( sp->st[0].Pop() >= 0. );
633 }
634 }
635 /* all solvers end up here */
636
637 /* check on the sum of the populations */
638 if( lgNegPop || dense.xIonDense[nelem][nelem-ipISO] < 0. )
639 {
640 fprintf( ioQQQ,
641 " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\
642 "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
643 ipISO,
644 nelem,
645 elementnames.chElementSym[nelem],
646 sp->chTypeAtomUsed,
647 dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO]),
648 sp->xIonSimple,
649 phycon.te,
650 nzone );
651
652 fprintf( ioQQQ, " level pop are:\n" );
653 for( i=0; i < numlevels_local; i++ )
654 {
655 fprintf( ioQQQ,PrintEfmt("%8.1e", sp->st[i].Pop() ));
656 if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
657 }
658 fprintf( ioQQQ, "\n" );
659 ContNegative();
660 ShowMe();
662 }
663
664 for( ipHi=1; ipHi<numlevels_local; ++ipHi )
665 {
666 for( ipLo=0; ipLo<ipHi; ++ipLo )
667 {
668 if( sp->trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
669 continue;
670
671 /* population of lower level, corrected for stimulated emission */
672 sp->trans(ipHi,ipLo).Emis().PopOpc() =
673 sp->st[ipLo].Pop() - sp->st[ipHi].Pop()*
674 sp->st[ipLo].g()/sp->st[ipHi].g();
675
676 // don't allow masers from collapsed levels or in Case B
677 if( N_(ipHi) > sp->n_HighestResolved_local || opac.lgCaseB )
678 sp->trans(ipHi,ipLo).Emis().PopOpc() = MAX2( 0., sp->trans(ipHi,ipLo).Emis().PopOpc() );
679 }
680 }
681
682 return;
683}
684
685void iso_set_ion_rates( long ipISO, long nelem)
686{
687 DEBUG_ENTRY("iso_set_ion_rates()");
688 t_iso_sp* sp = &iso_sp[ipISO][nelem];
689 const long int numlevels_local = sp->numLevels_local;
690 /* this is total ionization rate, s-1, of this species referenced to
691 * the total abundance */
692 double TotalPop = 0.;
693 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = 0.;
694 for( long level=0; level < numlevels_local; level++ )
695 {
696 /* sum of all ionization processes from this atom to ion, cm-3 s-1 now,
697 * but is divided by TotalPop below to become s-1 */
698 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] +=
699 sp->st[level].Pop() * sp->fb[level].RateLevel2Cont;
700 TotalPop += sp->st[level].Pop();
701 }
702
703 if( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] > BIGDOUBLE )
704 {
705 fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.",
706 nelem+1, nelem-ipISO);
708 }
709
710 if (TotalPop <= SMALLFLOAT)
711 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = sp->fb[0].RateLevel2Cont;
712 else
713 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] /= TotalPop;
714
715 if( ionbal.RateRecomIso[nelem][ipISO] > 0. )
716 sp->xIonSimple = ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1]/ionbal.RateRecomIso[nelem][ipISO];
717 else
718 sp->xIonSimple = 0.;
719
720 ASSERT( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] >= 0. );
721
722 if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 )
723 {
724 /* find fraction of He0 destructions due to photoionization of 2^3S */
725 double ratio;
726 double rateOutOf2TripS = sp->st[ipHe2s3S].Pop() * sp->fb[ipHe2s3S].RateLevel2Cont;
727 if( rateOutOf2TripS > SMALLFLOAT )
728 {
729 ratio = rateOutOf2TripS /
730 ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
731 }
732 else
733 {
734 ratio = 0.;
735 }
736 if( ratio > he.frac_he0dest_23S )
737 {
738 /* remember zone where this happended and fraction, and frac due to photoionization */
739 he.nzone = nzone;
740 he.frac_he0dest_23S = ratio;
741 rateOutOf2TripS = sp->st[ipHe2s3S].Pop() *sp->fb[ipHe2s3S].gamnc;
742 if( rateOutOf2TripS > SMALLFLOAT )
743 {
744 he.frac_he0dest_23S_photo = rateOutOf2TripS /
745 ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
746 }
747 else
748 {
749 he.frac_he0dest_23S_photo = 0.;
750 }
751 }
752 }
753}
t_atmdat atmdat
Definition atmdat.cpp:6
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
char TorF(bool l)
Definition cddefines.h:710
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
realnum rate_lu_nontherm() const
Definition collision.h:185
realnum ColUL(const ColliderList &colls) const
Definition collision.h:99
double & PopOpc() const
Definition emission.h:603
realnum & Pesc() const
Definition emission.h:523
realnum & Pelec_esc() const
Definition emission.h:533
realnum & Aul() const
Definition emission.h:613
double & pump() const
Definition emission.h:473
realnum & Pdest() const
Definition emission.h:543
CollisionProxy Coll() const
Definition transition.h:424
double EnergyRyd() const
Definition transition.h:83
EmissionList::reference Emis() const
Definition transition.h:408
ProxyIterator< qStateConstProxy, qStateConstProxy > const_iterator
iterator begin()
long int numLevels_max
Definition iso.h:493
multi_arr< extra_tr, 2 > ex
Definition iso.h:449
vector< freeBound > fb
Definition iso.h:452
double qTot2S
Definition iso.h:563
TransitionProxy trans(const long ipHi, const long ipLo)
Definition iso.h:444
long int numLevels_local
Definition iso.h:498
bool lgPopsRescaled
Definition iso.h:484
long int n_HighestResolved_local
Definition iso.h:507
char chTypeAtomUsed[10]
Definition iso.h:556
vector< two_photon > TwoNu
Definition iso.h:586
multi_arr< long, 3 > QuantumNumbers2Index
Definition iso.h:461
qList st
Definition iso.h:453
double xIonSimple
Definition iso.h:465
ColliderList colliders
Definition collision.cpp:57
void ContNegative(void)
t_conv conv
Definition conv.cpp:5
const double BIGDOUBLE
Definition cpu.h:194
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_dynamics dynamics
Definition dynamics.cpp:44
t_elementnames elementnames
GrainVar gv
Definition grainvar.cpp:5
t_he he
Definition he.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
#define IPCOLLIS
Definition iso.h:87
const int ipHe1s1S
Definition iso.h:41
#define IPRAD
Definition iso.h:86
const int ipH1s
Definition iso.h:27
const int ipHe2s3S
Definition iso.h:44
void iso_renorm(long nelem, long ipISO, double &renorm)
#define N_(A_)
Definition iso.h:20
const int ipHE_LIKE
Definition iso.h:63
const int ipH_LIKE
Definition iso.h:62
#define KILL_BELOW_PLASMA(E_)
Definition iso.h:17
void iso_level(const long int ipISO, const long int nelem, double &renorm)
Definition iso_level.cpp:31
void iso_set_ion_rates(long ipISO, long nelem)
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
t_opac opac
Definition opacity.cpp:5
#define S(I_, J_)
t_phycon phycon
Definition phycon.cpp:6
static double * source
Definition species2.cpp:28
static double * g
Definition species2.cpp:28
static double * sink
Definition species2.cpp:28
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
t_trace trace
Definition trace.cpp:5