cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_solve.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_solve main routine to call iso_level and determine iso level balances */
4 /* HydroRenorm - renormalize H so that it agrees with the chemistry */
5 /* AGN_He1_CS routine to punch table needed for AGN3 - collision strengths of HeI */
6 #include "cddefines.h"
7 #include "atmdat.h"
8 #include "conv.h"
9 #include "dense.h"
10 #include "opacity.h"
11 #include "elementnames.h"
12 #include "h2.h"
13 #include "helike.h"
14 #include "helike_cs.h"
15 #include "hmi.h"
16 #include "hydrogenic.h"
17 #include "ionbal.h"
18 #include "iso.h"
19 #include "phycon.h"
20 #include "secondaries.h"
21 #include "taulines.h"
22 #include "thermal.h"
23 #include "trace.h"
24 
25 /* iso_solve main routine to call iso_level and determine iso level balances */
26 void iso_solve(long ipISO)
27 {
28  long int ipLo,
29  ipHi,
30  nelem,
31  mol,
32  lowsav=-1,
33  ihisav=-1;
34  double coltot,
35  gamtot,
36  sum,
37  error;
38  static bool lgFinitePop[LIMELM];
39  static bool lgMustInit[NISO]={true, true};
40  bool lgH_chem_conv;
41  int loop_H_chem;
42  double solomon_assumed;
43  double *PumpSave=NULL;
44 
45  DEBUG_ENTRY( "iso_solve()" );
46 
47  if( lgMustInit[ipISO] )
48  {
49  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
50  lgFinitePop[nelem] = true;
51  }
52 
53  /*
54  * This is an option to account for intrinsic absorption or emission line the lyman
55  * lines. This is done without changing the incident coarse continuum.
56  * Check if continuum pumping of H Lyman lines to be multiplied by scale factor
57  * hydro.xLymanPumpingScaleFactor is set with atom h-like Lyman pumping scale command
58  * Lya pump rate is always updated - this is simplest way to finese intrinsic absorption
59  * or emission
60  */
61  if( ipISO==ipH_LIKE && hydro.xLymanPumpingScaleFactor!=1.f )
62  {
63  /* >> chng 06 aug 31, from numLevels_max to _local */
64  PumpSave = (double *)MALLOC( (unsigned)iso.numLevels_local[ipISO][ipHYDROGEN]*sizeof(double) );
65  ipLo = 0;
66  /* do not include the highest levels since pumping to them could create
67  * artificial ionizations */
68  /* >> chng 06 aug 31, from numLevels_max to _local */
69  for( ipHi=2; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi )
70  {
71  PumpSave[ipHi] = Transitions[ipISO][ipHYDROGEN][ipHi][ipLo].Emis->pump;
72  Transitions[ipISO][ipHYDROGEN][ipHi][ipLo].Emis->pump *= hydro.xLymanPumpingScaleFactor;
73  }
74  }
75 
76  if( ipISO == ipHE_LIKE )
77  {
78  /* save state of he-like low and high ionization, since must not change here,
79  * since there is a parallel helium atom that must decrement he */
81  fixit(); /* why does this routine not control helium? Shouldn't it? */
82  // can just remove these entirely?, only applies to helium itself anyway.
83  lowsav = dense.IonLow[ipHELIUM];
84  ihisav = dense.IonHigh[ipHELIUM];
85  }
86 
87  for( nelem=ipISO; nelem < LIMELM; nelem++ )
88  {
89  /* do not consider elements that have been turned off */
90  if( dense.lgElmtOn[nelem] )
91  {
92  /* note that nelem scale is totally on c not physical scale, so 0 is h */
93  /* evaluate balance if ionization reaches this high */
94  if( (dense.IonHigh[nelem] >= nelem + 1 - ipISO) &&
95  (dense.IonLow[nelem] <= nelem + 1 - ipISO) )
96  {
97  /* truncate atom if physical conditions limit the maximum principal quantum number of a
98  * bound electron to a number less than the malloc'd size */
99  iso_continuum_lower( ipISO, nelem );
100 
101  /* evaluate collisional rates */
102  iso_collide( ipISO, nelem );
103 
104  /* evaluate photoionization rates */
105  iso_photo( ipISO , nelem );
106 
107  fixit(); /* the order of error calculation is screwed up... rationalize this */
108  /* Generate Gaussian errors if turned on. */
109  if( iso.lgRandErrGen[ipISO] && nzone==0 && !iso.lgErrGenDone[ipISO][nelem] )
110  {
111  iso_error_generation(ipISO, nelem );
112  }
113 
114  /* evaluate recombination rates */
115  iso_radiative_recomb( ipISO , nelem );
116 
117  if( opac.lgRedoStatic )
118  {
119  if( nelem<=ipHELIUM )
120  {
121  iso_collapsed_bnl_set( ipISO, nelem );
122 
123  //iso_collapsed_bnl_print( ipISO, nelem );
124 
125  iso_collapsed_Aul_update( ipISO, nelem );
126 
127  iso_collapsed_lifetimes_update( ipISO, nelem );
128 
129  iso_cascade( ipISO, nelem );
130 
131  iso_radiative_recomb_effective( ipISO, nelem );
132  }
133  else
134  {
135  iso_cascade( ipISO, nelem );
136 
137  iso_radiative_recomb_effective( ipISO, nelem );
138  }
139  }
140 
141  /* evaluate state specific creation and destruction processes,
142  * also define iso.xIonSimple */
143  iso_ionize_recombine( ipISO , nelem );
144 
145  /* solve for the level populations */
146  iso_level( ipISO , nelem );
147 
148  /* this contains a bunch of trace statements and HydroRenorm.
149  * Will eventually come over to iso treatment. */
150  if( ipISO == ipH_LIKE )
151  HydroLevel(nelem);
152 
153  /* say that we have set the populations */
154  lgFinitePop[nelem] = true;
155 
156  }
157  else if( lgFinitePop[nelem] )
158  {
159  /* this branch, pops were set previously, but now are zero,
160  * so must zero them out this time */
161  lgFinitePop[nelem] = false;
162 
163  iso.pop_ion_ov_neut[ipISO][nelem] = 0.;
164  iso.xIonSimple[ipISO][nelem] = 0.;
165 
166  /* zero it out since no population*/
167  StatesElem[ipISO][nelem][0].Pop = 0.;
168  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
169  {
170  StatesElem[ipISO][nelem][ipHi].Pop = 0.;
171  for( ipLo=0; ipLo < ipHi; ipLo++ )
172  {
173  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
174  continue;
175 
176  /* population of lower level rel to ion, corrected for stim em */
177  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = 0.;
178  }
179  }
180  }
181  /* option to force ionization */
182  if( dense.lgSetIoniz[nelem] )
183  {
184  dense.xIonDense[nelem][nelem+1-ipISO] = dense.SetIoniz[nelem][nelem+1-ipISO]*dense.gas_phase[nelem];
185  dense.xIonDense[nelem][nelem-ipISO] = dense.SetIoniz[nelem][nelem-ipISO]*dense.gas_phase[nelem];
186  if( dense.SetIoniz[nelem][nelem+1-ipISO] > SMALLFLOAT )
187  {
188  StatesElem[ipISO][nelem][ipH1s].Pop = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO];
189  }
190  else
191  {
192  StatesElem[ipISO][nelem][ipH1s].Pop = 0.;
193  }
194  ASSERT( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Lo->Pop == StatesElem[ipISO][nelem][ipH1s].Pop );
195  }
196  }
197  }
198 
199  /* setting this flag false means that we have been through this loop at least once already. */
200  lgMustInit[ipISO] = false;
201 
202  if( ipISO == ipHE_LIKE )
203  {
204  dense.IonLow[ipHELIUM] = lowsav;
205  dense.IonHigh[ipHELIUM] = ihisav;
206  }
207 
208  if( ipISO==ipH_LIKE )
209  {
210  /* ============================================================================== */
211  /* rest is for hydrogen only */
212 
213  /* this block appears redundant and could be removed? */
214  /* >> 02 nov 21 rjrw -- xIonDense is used in hmole but only set in ion_solver,
215  * so we need this at least for the first iteration. */
216 
217  /* do molecular balance
218  * hmovh1 will be ratio of molecular to atomic hydrogen
219  * HIonFrac is fraction of H that is ionized, ratio of ion to atom */
220 
221  lgH_chem_conv = false;
222  loop_H_chem = 0;
223  while( loop_H_chem < 5 && !lgH_chem_conv )
224  {
225  /* save solomon rate that was assumed in the above - want to make sure that assumed
226  * solomon rate is stable, and does not change after call to large H2 molecule */
228 
229  /* now do chem, this will reset hmi.H2_Solomon_dissoc_rate_used */
230  hmole();
231  /* now do level populations for H2 */
232  H2_LevelPops();
233  /* >>chng 05 feb 12 - add logic checking that Solomon rate from big molecule is equal to
234  * rate used in H chem */
235  lgH_chem_conv = true;
236  /* check whether solomon rate has changed */
238  {
239  /* >>chng 05 mar 11, go from "total" to H2g for solomon rate */
240  if( fabs( solomon_assumed - hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.H2_rate_destroy) ) >
242  {
243  lgH_chem_conv = false;
244  }
245  }
246  ++loop_H_chem;
247  }
248 
249  {
250  /*@-redef@*/
251  /* often the H- route is the most efficient formation mechanism for H2,
252  * will be through rate called ratach
253  * this debug print statement is to trace h2 oscillations */
254  enum {DEBUG_LOC=false};
255  /*@+redef@*/
256  if(DEBUG_LOC )
257  {
258  fprintf(ioQQQ,"DEBUG \t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
259  fnzone,
260  hmi.H2_total ,
261  hmi.Hmolec[ipMH2g],
267  }
268  }
269  /* >>chng 01 may 09, add option to force abundance, with element name ioniz cmmnd */
271  {
274 
275  /* initialize ground state pop too */
277  }
278  else
279  {
280  /*
281  * >> chng 03 jan 15 rjrw:- terms are now in ion_solver, to allow for
282  * molecular sources and sinks of H and H+. ion_solver renormalizes
283  * to keep the total H abundance correct -- only the molecular
284  * network is allowed to change this.
285  */
286  ion_solver( ipHYDROGEN , false );
287  }
288 
289  /* confirm that species still add up correctly */
290  /* this exposes a "leak" that occurs somewhere, almost certainly in hmole
291  * if DEBUG_LOC is set true in the following there will be comments printing
292  * due to a loss of some protons */
293 
295  for(mol=0;mol<N_H_MOLEC;mol++)
296  {
297  sum += hmi.Hmolec[mol]*hmi.nProton[mol];
298  }
299  /* do not double count H0 and H+ */
300  sum -= hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp];
301 
302  /* >>chng 06 jun 30, remove H in CO network, as much as a few percent can be
303  * in the form of OH, H2O, etc
304  * >>chng 06 jul 03, undo this change so function is as before - need to
305  * think about co protons */
306  error = ( dense.gas_phase[ipHYDROGEN] - sum )/dense.gas_phase[ipHYDROGEN];
307  /*error = ( dense.gas_phase[ipHYDROGEN] - sum - dense.H_sum_in_CO)/dense.gas_phase[ipHYDROGEN];*/
308  {
309  /*@-redef@*/
310  /* often the H- route is the most efficient formation mechanism for H2,
311  * will be through rate called ratach
312  * this debug print statement is to trace h2 oscillations */
313  enum {DEBUG_LOC=false};
314  /*@+redef@*/
315  if(DEBUG_LOC && (fabs(error) > 1e-4) )
316  fprintf(ioQQQ,"PROBLEM hydrogenic zone %li hden %.4e, sum %.4e (h-s)/h %.3e \n", nzone, dense.gas_phase[ipHYDROGEN] , sum ,
317  error );
318  }
319 
320  fixit(); /* this is called in HydroLevel above, is it needed in both places? */
321  /* >>hcng 05 mar 24,
322  * renormalize the populations and emission of H atom to agree with chemistry */
323  HydroRenorm();
324 
325  /* this is test whether collisions are important first define ratio of exits
326  * from ground due to coll, rel to everthing else, then flag is large */
327  if( iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 )
328  {
329  coltot =
330  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL*
336  (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL*dense.eden +
337  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul*
338  (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pesc +
339  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pelec_esc +
340  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest) );
341 
342  /* caution that collisions are important (and so only small changes in
343  * temperature should happen) if more than 20% of the total */
344  if( coltot > 0.2 )
345  {
346  hydro.lgHColionImp = true;
347  }
348  }
349  else
350  {
351  hydro.lgHColionImp = false;
352  }
353 
354  /* remember the ratio of pops of 2p to 1s for possible printout in prtComment
355  * and to obtain Lya excitation temps. the pop of ground is not defined if
356  * NO ionization at all since these pops are relative to ion */
357  /* >>chng 99 jun 03, added MAX2 to protect against totally neutral gas */
360  {
361  hydro.lgHiPop2 = true;
363  hydro.pop2mx);
364  }
365 
367 
368  coltot = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s] +
369  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL*4.*
371 
372  /* if ground state destruction rate is significant, recall different dest procceses */
374  {
377 
378  /* fraction of ionizations of H from ground, due to thermal collisions */
381 
382  /* this flag is used in ConvBase to decide whether we
383  * really need to converge the secondary ionization rates */
386 
387  /* frac of ionizations due to ct */
389  }
390  else
391  {
393  hydro.H_ion_frac_photo = 0.;
394  secondaries.sec2total = 0.;
395  atmdat.HIonFrac = 0.;
396  }
397 
398  if( trace.lgTrace )
399  {
400  fprintf( ioQQQ, " Hydrogenic return %.2f ",fnzone);
401  fprintf(ioQQQ,"H0:%.3e ", dense.xIonDense[ipHYDROGEN][0]);
402  fprintf(ioQQQ,"H+:%.3e ", dense.xIonDense[ipHYDROGEN][1]);
403  fprintf(ioQQQ,"H2:%.3e ", hmi.H2_total);
404  fprintf(ioQQQ,"H-:%.3e ", hmi.Hmolec[ipMHm]);
405  fprintf(ioQQQ,"ne:%.3e ", dense.eden);
406  fprintf( ioQQQ, " REC, COL, GAMT= ");
407  /* recomb rate coef, cm^3 s-1 */
408  fprintf(ioQQQ,"%.2e ", iso.RadRec_effec[ipH_LIKE][ipHYDROGEN] );
409  fprintf(ioQQQ,"%.2e ", coltot);
410  fprintf(ioQQQ,"%.2e ", gamtot);
411  fprintf( ioQQQ, " CSUP=");
413  fprintf( ioQQQ, "\n");
414  }
415 
417  {
418  ipLo = 0;
419  /* do not include the highest levels since pumping to them could create
420  * artificial ionizaitons */
421  /* >> chng 06 aug 31, from numLevels_max to _local */
422  for( ipHi=2; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi )
423  {
424  Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Emis->pump = PumpSave[ipHi];
425  }
426  free(PumpSave);
427  }
428  }
429  return;
430 }
431 
432 /* HydroRenorm - renormalize H so that it agrees with the chemistry */
433 void HydroRenorm( void )
434 {
435  double sum_atom_iso , renorm;
436  long int level,
437  nelem,
438  ipHi ,
439  ipLo;
440 
441  DEBUG_ENTRY( "HydroRenorm()" );
442 
443  /*>>chng 04 mar 23, add this renorm */
444  /* renormalize the state specific populations, so that they
445  * add up to the results that came from ion_solver
446  * units at first is sum div by H+ density, since Pop2Ion defn this way */
447  sum_atom_iso = 0.;
448  /* >> chng 06 aug 31, from numLevels_max to _local */
449  for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; level++ )
450  {
451  sum_atom_iso += StatesElem[ipH_LIKE][ipHYDROGEN][level].Pop;
452  }
453  /* now convert to cm-3 - this is total population in iso solved model */
454  sum_atom_iso *= dense.xIonDense[ipHYDROGEN][1];
455 
456  /* >>chng 04 may 25, humunculus sum_atom_iso is zero */
457  if( sum_atom_iso > SMALLFLOAT )
458  {
459  renorm = dense.xIonDense[ipHYDROGEN][0] / sum_atom_iso;
460  }
461  else
462  {
463  renorm = 0.;
464  }
465  ASSERT( renorm < BIGFLOAT );
466  /*fprintf(ioQQQ,"DEBUG renorm\t%.2f\t%.3e\n",fnzone, renorm);*/
467  /* renormalize populations from iso-model atom so that they agree with ion solver & chemistry */
468  StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop *= renorm;
469  /*fprintf(ioQQQ,"DEBUG h \t%.3e hydrogenic renorm %.3e\n",
470  StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop ,
471  StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/renorm );*/
472 
473  nelem = ipHYDROGEN;
474  /* >> chng 06 aug 31, from numLevels_max to _local */
475  for( ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
476  {
477  StatesElem[ipH_LIKE][ipHYDROGEN][ipHi].Pop *= renorm;
478 
479  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
480  {
481  if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
482  continue;
483 
484  /* population of lower level rel to ion, corrected for stim em */
485  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->PopOpc *= renorm;
486  }
487  }
488 
489  return;
490 }
491 
492 /*iso_prt_pops print out iso sequence populations or departure coefficients */
493 void iso_prt_pops( long ipISO, long nelem, bool lgPrtDeparCoef )
494 {
495  long int in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN;
496  char chPrtType[2][12]={"populations","departure"};
497  /* first dimension is multiplicity */
498  char chSpin[3][9]= {"singlets", "doublets", "triplets"};
499 
500 #define ITEM_TO_PRINT(A_) ( lgPrtDeparCoef ? iso.DepartCoef[ipISO][nelem][A_] : StatesElem[ipISO][nelem][A_].Pop )
501 
502  DEBUG_ENTRY( "iso_prt_pops()" );
503 
504  ASSERT( ipISO < NISO );
505 
506  for( is = 1; is<=3; ++is)
507  {
508  if( ipISO == ipH_LIKE && is != 2 )
509  continue;
510  else if( ipISO == ipHE_LIKE && is != 1 && is != 3 )
511  continue;
512 
513  ipFirstCollapsed= iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem];
514  nResolved = StatesElem[ipISO][nelem][ipFirstCollapsed-1].n;
515  ASSERT( nResolved == iso.n_HighestResolved_local[ipISO][nelem] );
516  ASSERT(nResolved > 0 );
517 
518  /* give element number and spin */
519  fprintf(ioQQQ," %s %s %s %s\n",
520  iso.chISO[ipISO],
521  elementnames.chElementSym[nelem],
522  chSpin[is-1],
523  chPrtType[lgPrtDeparCoef]);
524 
525  /* header with the l states */
526  fprintf(ioQQQ," n\\l=> ");
527  for( i =0; i < nResolved; ++i)
528  {
529  fprintf(ioQQQ,"%2ld ",i);
530  }
531  fprintf(ioQQQ,"\n");
532 
533  /* loop over prin quant numbers, one per line, with l across */
534  for( in = 1; in <= nResolved; ++in)
535  {
536  if( is==3 && in==1 )
537  continue;
538 
539  fprintf(ioQQQ," %2ld ",in);
540 
541  for( il = 0; il < in; ++il)
542  {
543  if( ipISO==ipHE_LIKE && (in==2) && (il==1) && (is==3) )
544  {
545  fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P0) );
546  fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P1) );
547  fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P2) );
548  }
549  else
550  {
551  ipLo = iso.QuantumNumbers2Index[ipISO][nelem][in][il][is];
552  fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipLo) );
553  }
554  }
555  fprintf(ioQQQ,"\n");
556  }
557  }
558  /* above loop was over spin, now do collapsed levels, no spin or ang momen */
559  for( il = ipFirstCollapsed; il < iso.numLevels_local[ipISO][nelem]; ++il)
560  {
561  in = StatesElem[ipISO][nelem][il].n;
562  /* prin quan number of collapsed levels */
563  fprintf(ioQQQ," %2ld ",in);
564  fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(il) );
565  fprintf(ioQQQ,"\n");
566  }
567  return;
568 }
569 
570 /* routine to punch table needed for AGN3 - collision strengths of HeI */
571 void AGN_He1_CS( FILE *ioPun )
572 {
573 
574  long int i;
575 
576  /* list of temperatures where cs will be printed */
577  const int NTE = 5;
578  double TeList[NTE] = {6000.,10000.,15000.,20000.,25000.};
579  double TempSave;
580 
581  DEBUG_ENTRY( "AGN_He1_CS()" );
582 
583  /* put on a header */
584  fprintf(ioPun, "Te\t2 3s 33s\n");
585 
586  /* Restore the original temp when this routine done. */
587  TempSave = phycon.te;
588 
589  for( i=0; i<NTE; ++i )
590  {
591  TempChange(TeList[i] , false);
592 
593  fprintf(ioPun , "%.0f\t",
594  TeList[i] );
595  fprintf(ioPun , "%.2f\t",
597  fprintf(ioPun , "%.2f\t",
599  fprintf(ioPun , "%.2f\t",
601  fprintf(ioPun , "%.3f\t",
603  /*fprintf(ioPun , "%.1f\t%.1f\t%.1f\n", */
604  fprintf(ioPun , "%.1f\n",
608  }
609 
610  /* no need to force update since didn't do above */
611  TempChange(TempSave , false);
612  return;
613 }

Generated for cloudy by doxygen 1.8.1.2