cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_level.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 /*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 "helike.h"
14 #include "hydrogenic.h"
15 #include "ionbal.h"
16 #include "iso.h"
17 #include "mole.h"
18 #include "phycon.h"
19 #include "physconst.h"
20 #include "rfield.h"
21 #include "secondaries.h"
22 #include "taulines.h"
23 #include "thirdparty.h"
24 #include "trace.h"
25 
26 /* solve for level populations */
27 void iso_level( const long int ipISO, const long int nelem)
28 {
29  long int ipHi,
30  ipLo,
31  i,
32  level,
33  level_error;
34 
35  const long int numlevels_local = iso.numLevels_local[ipISO][nelem];
36 
37  double BigError;
38 
39  int32 nerror;
40  double HighestPColOut = 0.,
41  collider,
42  sum_popn_ov_ion,
43  TooSmall;
44  bool lgNegPop=false;
45  valarray<int32> ipiv(numlevels_local);
46  /* this block of variables will be obtained and freed here */
47  valarray<double>
48  creation(numlevels_local),
49  error(numlevels_local),
50  work(numlevels_local),
51  Save_creation(numlevels_local);
52  double source=0.,
53  sink=0.;
54  valarray<double> PopPerN(iso.n_HighestResolved_local[ipISO][nelem]+1);
55 
57 
58  DEBUG_ENTRY( "iso_level()" );
59 
60  /* check that we were called with valid charge */
61  ASSERT( nelem >= ipISO );
62  ASSERT( nelem < LIMELM );
63 
64  /* now do the 2D array */
65  z.alloc(numlevels_local,numlevels_local);
66 
67  /* >>chng 05 aug 17, must use real electron density for collisional ionization of H
68  * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
69  * H on H is homonuclear and scaling laws for other elements does not apply
70  * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
71  * EdenHCorr for rest */
72  if( ipISO == ipH_LIKE && nelem==ipHYDROGEN )
73  {
74  /* special version for H0 onto H0 */
76  }
77  else
78  {
80  }
81 
82  /* fill in recombination vector - values were set in iso_ionize_recombine.cpp */
83  for( level=0; level < numlevels_local; level++ )
84  {
85  /* total recombination to level [s-1] */
86  creation[level] = iso.RateCont2Level[ipISO][nelem][level];
87  }
88 
89  /* these two collision rates must be the same or we are in big trouble,
90  * since used interchangably */
91  ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
92  fabs( iso.ColIoniz[ipISO][nelem][0]* collider /
93  SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );
94 
95 
96  if( ipISO == ipH_LIKE )
97  {
98  if( nelem==ipHYDROGEN )
99  TooSmall = 1e-28;
100  else if( nelem==ipHELIUM )
101  TooSmall = 1e-18;
102  else
103  TooSmall = 2e-18;
104  }
105  else if( ipISO == ipHE_LIKE )
106  {
107  /* >>chng 03 apr 30, different limit for He itself,
108  * since so important in molecular regions */
109  if( nelem==ipHELIUM )
110  TooSmall = 1e-20;
111  else
112  TooSmall = 1e-10;
113  }
114  else
115  TotalInsanity();
116 
117  /* which case atom to solve??? */
118  /* >>chng 03 apr 11, from ae-30 bail, to ae-35, will catch it in hydro 2 low */
119  if( ipISO == ipH_LIKE &&
120  ( iso.xIonSimple[ipISO][nelem] < 1e-35 ) )
121  {
122  /* don't bother if no ionizing radiation */
123  strcpy( iso.chTypeAtomUsed[ipISO][nelem], "zero " );
124  if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
125  {
126  fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n",
127  ipISO, nelem, iso.xIonSimple[ipISO][nelem], iso.chTypeAtomUsed[ipISO][nelem] );
128  }
129 
130  /* okay to leave as iso.numLevels_max */
131  for( i=0; i < iso.numLevels_max[ipISO][nelem]; i++ )
132  {
133  iso.DepartCoef[ipISO][nelem][i] = 1.;
134  StatesElem[ipISO][nelem][i].Pop = 0.;
135  }
136 
137  ionbal.RateRecomTot[nelem][nelem-ipISO] = 0.;
138  iso.xIonSimple[ipISO][nelem] = 0.;
139  iso.pop_ion_ov_neut[ipISO][nelem] = 0.;
140  iso.qTot2S[ipISO][nelem] = 0.;
141  }
142  /* >>chng 99 nov 23, very dense model close to lte had very low simple
143  * ionization fraction but moderate ionization since all pops in
144  * excited states - added test for density
145  * second change - for all cases use this routine if ionization is
146  * very low indeed */
147  /* >>chng 00 aug 26, added matrix force option,
148  * logic is to override default options, "none was set in zero */
149 
150  /* NB - this test is mostly bogus since chTypeAtom is "POPU" in zero */
151 
152  /* >>chng 02 mar 13, iso.chTypeAtomSet had been set to POPU in zero, look for
153  * comment with this date. This had the effect of killing this following test.
154  * This also meant that the code had been happily inverting these impossible
155  * matrices for quite some time. This test changed to stringest one, in
156  * light of this */
157  else if( ipISO == ipH_LIKE &&
158  ((strcmp( iso.chTypeAtomSet[ipISO] , "LOWT" )==0) ||
159  (iso.xIonSimple[ipISO][nelem] < TooSmall) ) )
160  {
161  strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Low T" );
162  if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
163  {
164  fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s\n",
165  ipISO, nelem, iso.xIonSimple[ipISO][nelem], iso.chTypeAtomUsed[ipISO][nelem] );
166  }
167 
168  /* the ionization ratio WILL be equal to iso.xIonSimple[ipISO][nelem] */
169  iso.pop_ion_ov_neut[ipISO][nelem] = iso.xIonSimple[ipISO][nelem];
170 
171  /* do simple cascade, should always work,
172  * but not very well at high densities */
173  HydroT2Low( ipISO, nelem );
174  iso.qTot2S[ipISO][nelem] = 0.;
175  }
176  /* branch for very little ionization, use simple estimate but do not do level populations
177  * never set to zero, use simple two-level system, since He+ important for chemistry */
178  else if( ipISO == ipHE_LIKE &&
179  ((strcmp( iso.chTypeAtomSet[ipISO] , "LOWT" )==0) ||
180  (iso.xIonSimple[ipISO][nelem] < TooSmall) ) )
181  {
182  strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Low T" );
183  /* total ionization will just be the ground state */
184  ionbal.RateIonizTot[nelem][nelem-ipISO] = iso.RateLevel2Cont[ipISO][nelem][0];
185  /* >>chng 01 nov 23, increase upper limit by 2 */
186  lgNegPop = false;
187  /* >>chng 02 apr 10, had set to zero,
188  * we will leave the level populations at zero but define an ion ratio,
189  * since He+ abundance cannot go to zero in a molecular cloud,
190  * - He+ charge transfer is the main CO destruction mechanism */
191  iso.pop_ion_ov_neut[ipISO][nelem] = iso.xIonSimple[ipISO][nelem];
192  StatesElem[ipISO][nelem][0].Pop = 1./SDIV(iso.pop_ion_ov_neut[ipISO][nelem]);
193  for( long n=1; n < numlevels_local; n++ )
194  {
195  StatesElem[ipISO][nelem][n].Pop = 0.;
196  }
197  iso.qTot2S[ipISO][nelem] = 0.;
198  }
199  else
200  {
201  ASSERT( NISO == 2 );
202  long SpinUsed[NISO] = {2,1};
203  long indexNmaxP =
204  iso.QuantumNumbers2Index[ipISO][nelem][ iso.n_HighestResolved_local[ipISO][nelem] ][1][SpinUsed[ipISO]];
205 
206  strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Popul" );
207  if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
208  {
209  fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n",
210  ipISO, nelem, iso.chTypeAtomUsed[ipISO][nelem] );
211  }
212 
213  multi_arr<quantumState,3>::const_iterator StElm = StatesElem.begin(ipISO,nelem);
214 
215  /* master balance equation, use when significant population */
216  for( level=0; level < numlevels_local; level++ )
217  {
218  /* all process depopulating level and placing into the continuum
219  * this does NOT include grain charge transfer ionization, added below */
220  z[level][level] = iso.RateLevel2Cont[ipISO][nelem][level];
221 
222  if( level == 1 )
223  /* >>chng 05 dec 21, rm eden to make into rate coefficient */
224  iso.qTot2S[ipISO][nelem] = iso.ColIoniz[ipISO][nelem][level];
225 
226  multi_arr<transition,4>::const_iterator Trans = Transitions.begin(ipISO,nelem,level);
227  md4ci Boltz = iso.Boltzmann.begin(ipISO,nelem,level);
228 
229  /* all processes populating level from below */
230  for( ipLo=0; ipLo < level; ipLo++ )
231  {
232  double coll_down = Trans[ipLo].Coll.ColUL * collider;
233 
234  double RadDecay = MAX2( iso.SmallA, Trans[ipLo].Emis->Aul*
235  (Trans[ipLo].Emis->Pesc +
236  Trans[ipLo].Emis->Pelec_esc +
237  Trans[ipLo].Emis->Pdest)*
238  KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) );
239 
240  double pump = MAX2( iso.SmallA, Trans[ipLo].Emis->pump *
241  KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) );
242 
243  if( iso.lgRandErrGen[ipISO] )
244  {
245  coll_down *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPCOLLIS];
246  RadDecay *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD];
247  pump *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD];
248  }
249 
250  double coll_up = coll_down *
251  (double)StElm[level].g/
252  (double)StElm[ipLo].g*
253  Boltz[ipLo];
254 
255  z[ipLo][ipLo] += coll_up + pump ;
256  z[ipLo][level] = - ( coll_up + pump );
257 
258  double pump_down = pump *
259  (double)StElm[ipLo].g/
260  (double)StElm[level].g;
261 
262  z[level][level] += RadDecay + pump_down + coll_down;
263  z[level][ipLo] = - (RadDecay + pump_down + coll_down);
264 
265  if( level == indexNmaxP )
266  {
267  HighestPColOut += coll_down;
268  }
269  if( ipLo == indexNmaxP )
270  {
271  HighestPColOut += coll_up;
272  }
273 
274  /* find total collisions out of 2^3S to singlets. */
275  if( (level == 1) && (ipLo==0) )
276  {
277  iso.qTot2S[ipISO][nelem] += coll_down/dense.EdenHCorr;
278  }
279  if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S==1) )
280  {
281  iso.qTot2S[ipISO][nelem] += coll_up/dense.EdenHCorr;
282  }
283  }
284  }
285 
286  if( ipISO == nelem )
287  {
288  /* iso.lgCritDensLMix[ipISO] is a flag used to print warning if density is
289  * too low for first collapsed level to be l-mixed. Check is if l-mixing collisions
290  * out of highest resolved singlet P are greater than sum of transition probs out. */
291  if( HighestPColOut < 1./StatesElem[ipISO][nelem][indexNmaxP].lifetime )
292  {
293  iso.lgCritDensLMix[ipISO] = false;
294  }
295  }
296 
298  ASSERT( ipISO <= ipHE_LIKE );
299 
300  /* induced two photon emission - special because upward and downward are
301  * not related by ratio of statistical weights */
302  /* iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command */
303  z[1+ipISO][0] -= iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On;
304  z[0][1+ipISO] -= iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On;
305 
306  /* rates out of 1s, and out of 2s */
307  z[0][0] += iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On;
308  z[1+ipISO][1+ipISO] += iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On;
309 
310  /* grain charge transfer recombination and ionization to ALL other stages */
311  for( long ion=0; ion<=nelem+1; ++ion )
312  {
313  if( ion!=nelem-ipISO )
314  {
315  /* recombination must be multiplied by a ratio of densities to get proper rate. */
316  source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] *
317  dense.xIonDense[nelem][ion] ;
318  /* take ionization out of every level. */
319  sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
320  }
321  }
322 
323  /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
324  /*>>>chng 02 oct 01, only include if lgAdvection is set */
325  if( dynamics.lgAdvection && dynamics.lgISO[ipISO])
326  {
327  /* add in advection - these terms normally zero */
328  /* assume for present that all advection is into ground state */
329  source += dynamics.Source[nelem][nelem-ipISO];
330  /* >>chng 02 Sep 06 rjrw -- advective term not recombination */
331  /* can sink from all components (must do, for conservation) */
332  sink += dynamics.Rate;
333  }
334 
335 #if 0
336  /* add in source and sink terms from molecular network. */
337  CodeReview(); /* Check... */
338  if( nelem == ipISO && conv.nTotalIoniz )
339  {
340  /* these are the external source and sink terms */
341  /* source first */
342  source += mole.source[nelem][nelem-ipISO];
343  sink += mole.sink[nelem][nelem-ipISO];
344  }
345 #endif
346 
347  creation[0] += source/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
348  for( level=0; level < numlevels_local; level++ )
349  {
350  z[level][level] += sink;
351  }
352 
353  /* >>chng 04 nov 30, atom XX-like collisions off turns this off too */
354  if( secondaries.Hx12[ipISO][nelem][iso.nLyaLevel[ipISO]] * iso.lgColl_excite[ipISO] > 0. )
355  {
356  /* now add on supra thermal excitation */
357  for( level=1; level < numlevels_local; level++ )
358  {
359  double RateUp , RateDown;
360 
361  RateUp = secondaries.Hx12[ipISO][nelem][level];
362  RateDown = RateUp * (double)StatesElem[ipISO][nelem][ipH1s].g /
363  (double)StatesElem[ipISO][nelem][level].g;
364 
365  /* total rate out of lower level */
366  z[ipH1s][ipH1s] += RateUp;
367 
368  /* rate from the upper level to ground */
369  z[level][ipH1s] -= RateDown;
370 
371  /* rate from ground to upper level */
372  z[ipH1s][level] -= RateUp;
373 
374  z[level][level] += RateDown;
375  }
376  }
377 
378  /* ===================================================================
379  *
380  * at this point all matrix elements have been established
381  *
382  * ==================================================================== */
383 
384  /* save matrix, this allocates SaveZ */
385  SaveZ = z;
386 
387  for( ipLo=0; ipLo < numlevels_local; ipLo++ )
388  Save_creation[ipLo] = creation[ipLo];
389 
390  if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
391  {
392  fprintf( ioQQQ, " pop level others => (iso_level)\n" );
393  for( long n=0; n < numlevels_local; n++ )
394  {
395  fprintf( ioQQQ, " %s %s %2ld", iso.chISO[ipISO], elementnames.chElementNameShort[nelem], n );
396  for( long j=0; j < numlevels_local; j++ )
397  {
398  fprintf( ioQQQ,"\t%.9e", z[j][n] );
399  }
400  fprintf( ioQQQ, "\n" );
401  }
402  fprintf(ioQQQ," recomb ");
403  for( long n=0; n < numlevels_local; n++ )
404  {
405  fprintf( ioQQQ,"\t%.9e", creation[n] );
406  }
407  fprintf( ioQQQ, "\n" );
408  fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
410  findspecies("CO")->hevmol ,
411  atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] ,
412  atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0] );
413  }
414 
415  nerror = 0;
416 
417  getrf_wrapper(numlevels_local,numlevels_local,
418  z.data(),numlevels_local,&ipiv[0],&nerror);
419 
420  getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0],
421  &creation[0],numlevels_local,&nerror);
422 
423  if( nerror != 0 )
424  {
425  fprintf( ioQQQ, " iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
426  cdEXIT(EXIT_FAILURE);
427  }
428 
429  /* check whether solution is valid */
430  /* >>chng 06 aug 28, both of these from numLevels_max to _local. */
431  for( level=ipH1s; level < numlevels_local; level++ )
432  {
433  double BigSoln= 0.;
434  error[level] = 0.;
435  for( long n=ipH1s; n < numlevels_local; n++ )
436  {
437  /* remember the largest value of the soln matrix to div by below */
438  if( fabs(SaveZ[n][level] ) > BigSoln )
439  BigSoln = fabs(SaveZ[n][level]);
440 
441  error[level] += SaveZ[n][level]*creation[n];
442  }
443 
444  if( BigSoln > 0. )
445  {
446  error[level] = (error[level] - Save_creation[level])/ BigSoln;
447  }
448  else
449  {
450  error[level] = 0.;
451  }
452  }
453 
454  /* remember largest residual in matrix inversion */
455  BigError = -1.;
456  level_error = -1;
457  /* >>chng 06 aug 28, from numLevels_max to _local. */
458  for( level=ipH1s; level < numlevels_local; level++ )
459  {
460  double abserror;
461  abserror = fabs( error[level]);
462  /* this will be the largest residual in the matrix inversion */
463  if( abserror > BigError )
464  {
465  BigError = abserror;
466  level_error = level;
467  }
468  }
469 
470  /* matrix inversion should be nearly as good as the accuracy of a double,
471  * but demand that it is better than epsilon for a float */
472  if( BigError > FLT_EPSILON )
473  {
474  if( !conv.lgSearch )
475  fprintf(ioQQQ," PROBLEM" );
476 
477  fprintf(ioQQQ,
478  " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
479  "level was %li Search?%c \n",
480  fnzone,
481  ipISO,
483  nelem ,
484  BigError ,
485  level_error,
486  TorF(conv.lgSearch) );
487  }
488 
489  /* put departure coefficients and level populations into master array */
490  for( level=0; level < numlevels_local; level++ )
491  {
492  StatesElem[ipISO][nelem][level].Pop = creation[level];
493 
494  /* check for negative level populations */
495  if( StatesElem[ipISO][nelem][level].Pop < 0. )
496  lgNegPop = true;
497 
498  if( StatesElem[ipISO][nelem][level].Pop <= 0 )
499  {
500  fprintf(ioQQQ," non-positive level pop for iso = %li, nelem = %li = %s, level=%li val=%.3e\n",
501  ipISO,
502  nelem ,
503  elementnames.chElementSym[nelem],
504  level,
505  StatesElem[ipISO][nelem][level].Pop );
506  }
507 
508  if( iso.PopLTE[ipISO][nelem][level] > 0. )
509  {
510  /* the LTE departure coefficients */
511  iso.DepartCoef[ipISO][nelem][level] =
512  StatesElem[ipISO][nelem][level].Pop/
513  (iso.PopLTE[ipISO][nelem][level]* dense.eden);
514  }
515  else
516  {
517  iso.DepartCoef[ipISO][nelem][level] = 0.;
518  }
519  }
520 
521  /* zero populations of unused levels. */
522  for( level=numlevels_local; level < iso.numLevels_max[ipISO][nelem]; level++ )
523  {
524  StatesElem[ipISO][nelem][level].Pop = 0.;
525  /* >>chng 06 jul 25, no need to zero this out, fix limit to 3-body heating elsewhere. */
526  /* iso.PopLTE[ipISO][nelem][level] = 0.; */
527  }
528  }
529 
530  /* all solvers end up here */
531  /* sum_popn_ov_ion will become the ratio of iso to parent
532  * species, create sum of level pops per ion first */
533  sum_popn_ov_ion = 0.;
534 
535  /* this is total ionization of this species referenced to the ground state */
536  ionbal.RateIonizTot[nelem][nelem-ipISO] = 0.;
537 
538  for( level=0; level < numlevels_local; level++ )
539  {
540  /* sum of all ionization processes from this atom to ion */
541  ionbal.RateIonizTot[nelem][nelem-ipISO] +=
542  StatesElem[ipISO][nelem][level].Pop * iso.RateLevel2Cont[ipISO][nelem][level];
543 
544  /* create sum of populations relative to ion */
545  sum_popn_ov_ion += StatesElem[ipISO][nelem][level].Pop;
546  }
547 
548  /* convert back to scaled from ground */
549  if( ( ionbal.RateIonizTot[nelem][nelem-ipISO]/MAX2(SMALLDOUBLE , sum_popn_ov_ion) ) > BIGDOUBLE )
550  {
551  fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.",
552  nelem+1, nelem-ipISO);
553  cdEXIT(EXIT_FAILURE);
554  }
555  else
556  ionbal.RateIonizTot[nelem][nelem-ipISO] /= MAX2(SMALLDOUBLE , sum_popn_ov_ion);
557 
558  /* this is sum of all populations in model, div by ion
559  * create ratio of ion pops to total atom */
560  if( sum_popn_ov_ion >= 1e-30 )
561  {
562  iso.pop_ion_ov_neut[ipISO][nelem] = 1./sum_popn_ov_ion;
563  }
564  else
565  {
566  iso.pop_ion_ov_neut[ipISO][nelem] = 0.;
567 
568  /* >>chng 03 aug 16, zero out parent density - had not been done */
569  dense.xIonDense[nelem][nelem+1-ipISO] = 0.;
570 
571  /* reset pointer to one lower stage of ionization so this not
572  * considered again, hydrogenic considered if IonHigh is nelem+2 */
573  dense.IonHigh[nelem] = nelem;
574  ionbal.RateIonizTot[nelem][nelem-ipISO] = 0.;
575 
576  /* now zero this out */
577  /* okay to leave as iso.numLevels_max */
578  for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
579  {
580  /* population of level relative to ion */
581  StatesElem[ipISO][nelem][ipHi].Pop = 0.;
582  }
583  }
584 
585  ASSERT( ionbal.RateIonizTot[nelem][nelem-ipISO] >= 0. );
586 
587  /* check on the sum of the populations */
588  if( lgNegPop || iso.pop_ion_ov_neut[ipISO][nelem] < 0. )
589  {
590  fprintf( ioQQQ,
591  " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\
592  "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
593  ipISO,
594  nelem,
595  elementnames.chElementSym[nelem],
596  iso.chTypeAtomUsed[ipISO][nelem],
597  iso.pop_ion_ov_neut[ipISO][nelem],
598  iso.xIonSimple[ipISO][nelem],
599  phycon.te,
600  nzone );
601 
602  fprintf( ioQQQ, " level pop are:\n" );
603  for( i=0; i < numlevels_local; i++ )
604  {
605  fprintf( ioQQQ,PrintEfmt("%8.1e", StatesElem[ipISO][nelem][i].Pop ));
606  if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
607  }
608  fprintf( ioQQQ, "\n" );
609  ContNegative();
610  ShowMe();
611  cdEXIT(EXIT_FAILURE);
612  }
613 
614  if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 )
615  {
616  /* find fraction of He0 destructions due to photoionization of 2^3S */
617  double ratio;
618  double rateOutOf2TripS = StatesElem[ipISO][nelem][ipHe2s3S].Pop * iso.RateLevel2Cont[ipISO][nelem][ipHe2s3S];
619  if( rateOutOf2TripS > SMALLFLOAT )
620  {
621  ratio = rateOutOf2TripS /
622  ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] );
623  }
624  else
625  {
626  ratio = 0.;
627  }
628  if( ratio > he.frac_he0dest_23S )
629  {
630  /* remember zone where this happended and fraction, and frac due to photoionization */
631  he.nzone = nzone;
632  he.frac_he0dest_23S = ratio;
633  rateOutOf2TripS = StatesElem[ipISO][nelem][ipHe2s3S].Pop *iso.gamnc[ipISO][nelem][ipHe2s3S];
634  if( rateOutOf2TripS > SMALLFLOAT )
635  {
636  he.frac_he0dest_23S_photo = rateOutOf2TripS /
637  ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] );
638  }
639  else
640  {
642  }
643  }
644  }
645 
646  for( ipHi=1; ipHi<numlevels_local; ++ipHi )
647  {
648  for( ipLo=0; ipLo<ipHi; ++ipLo )
649  {
650  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
651  continue;
652 
653  /* population of lower level rel to ion, corrected for stim em */
654  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc =
655  StatesElem[ipISO][nelem][ipLo].Pop - StatesElem[ipISO][nelem][ipHi].Pop*
656  StatesElem[ipISO][nelem][ipLo].g/StatesElem[ipISO][nelem][ipHi].g;
657 
658  /* don't allow masers from collapsed levels */
659  if( N_(ipHi) > iso.n_HighestResolved_local[ipISO][nelem] )
660  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = StatesElem[ipISO][nelem][ipLo].Pop;
661  }
662  }
663 
664  return;
665 }

Generated for cloudy by doxygen 1.8.4