cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_h2_form.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*mole_H2_form find state specific rates grains and H- form H2 */
4 #include "cddefines.h"
5 #include "grainvar.h"
6 #include "phycon.h"
7 #include "hmi.h"
8 #include "dense.h"
9 #include "h2.h"
10 #include "deuterium.h"
11 #include "mole.h"
12 
13 /*mole_H2_form find state specific rates grains and H- form H2 */
15 {
16  DEBUG_ENTRY( "mole_H2_form()" );
17 
18  /* rate of entry into X from H- and formation on grain surfaces
19  * will one of several distribution functions derived elsewhere
20  * first zero out formation rates and rates others collide into particular level */
21  for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
22  {
23  for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
24  {
25  /* this will be the rate formation (s-1) of H2 due to
26  * both formation on grain surfaces and the H minus route,
27  * also H2+ + H => H2 + H+ into one vJ level */
28  /* units cm-3 s-1 */
29  H2_X_formation[iVib][iRot] = 0.;
30  H2_X_Hmin_back[iVib][iRot] = 0.;
31  }
32  }
33 
34  /* loop over all grain types, finding total formation rate into each ro-vib level,
35  * also keeps trace of total formation into H2 ground and star, as defined by Tielens & Hollenbach,
36  * these are used in the H molecular network */
37  hmi.H2_forms_grains = 0.;
39 
40  double sum_check = 0.;
41 
42  /* >>chng 02 oct 08, resolved grain types */
43  for( size_t nd=0; nd < gv.bin.size(); ++nd )
44  {
45  int ipH2 = (int)gv.which_H2distr[gv.bin[nd]->matType];
46  for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
47  {
48  for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
49  {
50  /* >>chng 02 nov 14, changed indexing into H2_X_grain_formation_distribution and gv.bin, PvH */
51  double one =
52  /* H2_X_grain_formation_distribution is normalized to a summed total of unity */
54  /* units of following are s-1 */
55  gv.bin[nd]->rate_h2_form_grains_used;
56 
57  sum_check += one;
58 
59  /* final units are s-1 */
60  /* units cm-3 s-1 */
61  /* >>chng 04 may 05, added atomic hydrogen density, units cm-3 s-1 */
62  H2_X_formation[iVib][iRot] += (realnum)one*dense.xIonDense[ipHYDROGEN][0];
63 
64  /* save rates for formation into "H2" and "H2*" in the chemical
65  * network - it resolves the H2 into two species, as in
66  * Hollenbach / Tielens work - these rates will be used in the
67  * chemistry solver to get H2 and H2* densities */
68  if( states[ ipEnergySort[0][iVib][iRot] ].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
69  {
70  hmi.H2star_forms_grains += one;
71  }
72  else
73  {
74  /* unit s-1, h2 means h2g*/
75  hmi.H2_forms_grains += one;
76  }
77  }
78  }
79  }
80 
81  ASSERT( fp_equal_tol( sum_check, gv.rate_h2_form_grains_used_total, 1e-5*sum_check + DBL_MIN ) );
82 
83  /* convert to dimensionless factors that add to unity */
84  /* >>chng 02 oct 17, use proper distribution function */
86  hmi.H2_forms_hminus = 0.;
87  double frac_lo , frac_hi;
88  long ipT;
89 
90  /* which temperature point to use? */
91  if( phycon.alogte<=1. )
92  {
93  ipT = 0;
94  frac_lo = 1.;
95  frac_hi = 0.;
96  }
97  else if( phycon.alogte>=4. )
98  {
99  ipT = nTE_HMINUS-2;
100  frac_lo = 0.;
101  frac_hi = 1.;
102  }
103  else
104  {
105  /* find the temp */
106  for( ipT=0; ipT<nTE_HMINUS-1; ++ipT )
107  {
108  if( H2_logte_hminus[ipT+1]>phycon.alogte )
109  break;
110  }
111  frac_hi = (phycon.alogte-H2_logte_hminus[ipT])/(H2_logte_hminus[ipT+1]-H2_logte_hminus[ipT]);
112  frac_lo = 1.-frac_hi;
113  }
114 
115  /* formation of H2 in excited states from H- H minus */
116  /* >>chng 02 oct 17, temp dependent fit to rate, updated reference,
117  * about 40% larger than before */
118  /* rate in has units cm-3 s-1 */
119  double rate = (mole.findrk("H,H-=>H2,e-")+mole.findrk("H,H-=>H2*,e-")) * dense.xIonDense[ipHYDROGEN][0] * findspecieslocal("H-")->den;
120  /* backward rate, s-1 */
121  double rateback = (mole.findrk("H2,e-=>H,H-")+mole.findrk("H2*,e-=>H,H-"))*dense.eden;
122  double oldrate = 0.;
123 
124  /* we now know how to interpolate, now fill in H- formation sites */
125  for( long iVib=0; iVib<=nVib_hi[0]; ++iVib )
126  {
127  for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
128  {
129  /* the temperature-interpolated distribution function, adds up to unity,
130  * dimensionless */
131  double rate_interp =
132  frac_lo*H2_X_hminus_formation_distribution[ipT][iVib][iRot] +
133  frac_hi*H2_X_hminus_formation_distribution[ipT+1][iVib][iRot];
134 
135  /* above rate was set, had dimensions cm-3 s-1
136  * rate is product of parent densities and total formation rate */
137  double one = rate * rate_interp;
138 
139  /* save this rate [s-1] for back reaction in levels solver for v,J */
140  H2_X_Hmin_back[iVib][iRot] = (realnum)(rateback * rate_interp);
141 
142  /* units cm-3 s-1 */
143  H2_X_formation[iVib][iRot] += (realnum)one;
144 
145  oldrate += rate_interp;
146 
147  /* save rates to pass back into molecule network */
148  if( states[ ipEnergySort[0][iVib][iRot] ].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
149  {
150  hmi.H2star_forms_hminus += one;
151  }
152  else
153  {
154  /* unit cm-3s-1, h2 means h2g*/
155  hmi.H2_forms_hminus += one;
156  }
157  }
158  }
159  /* confirm that shape function is normalized correctly */
160  ASSERT( fabs(1.-oldrate)<1e-4 );
161 
162  /* >>chng 03 feb 10, add this population process */
163  /* H2+ + H => H2 + H+,
164  * >>refer H2 population Krstic, P.S., preprint
165  * all goes into v=4 but no J information, assume into J = 0 */
166  /* >>chng 04 may 05, add density at end */
167  rate = mole.findrk("H,H2+=>H+,H2*") * dense.xIonDense[ipHYDROGEN][0] * findspecieslocal("H2+")->den;
168  /* units cm-3 s-1 */
169  H2_X_formation[4][0] += (realnum)rate;
170 
171  fixit("this code is still far too H2-centric. Kick the can a bit");
172  ASSERT( this==&h2 || this==&hd );
173  if( this != &h2 )
174  {
175  for( long iVib=0; iVib<=nVib_hi[0]; ++iVib )
176  {
177  for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
178  {
179  // rescale everything in this routine by D/H
180  // THIS MUST NOT BE KEPT FOR OTHER DIATOMS!!!!
182  }
183  }
184  }
185 
186  return;
187 }
const double ENERGY_H2_STAR
Definition: h2_priv.h:449
bool lgLeiden_Keep_ipMH2s
Definition: hmi.h:219
double H2star_forms_grains
Definition: hmi.h:164
double H2_forms_grains
Definition: hmi.h:164
long int Jlowest[N_ELEC]
Definition: h2_priv.h:480
const int nTE_HMINUS
Definition: h2_priv.h:31
double findrk(const char buf[]) const
t_phycon phycon
Definition: phycon.cpp:6
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:908
double H2_forms_hminus
Definition: hmi.h:164
molezone * findspecieslocal(const char buf[])
t_dense dense
Definition: global.cpp:15
multi_arr< realnum, 2 > H2_X_formation
Definition: h2_priv.h:517
const realnum H2_logte_hminus[nTE_HMINUS]
Definition: h2_priv.h:37
H2_type which_H2distr[MAT_TOP]
Definition: grainvar.h:523
multi_arr< long int, 3 > ipEnergySort
Definition: h2_priv.h:551
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
void mole_H2_form(void)
multi_arr< realnum, 3 > H2_X_hminus_formation_distribution
Definition: h2_priv.h:546
double energy(const genericState &gs)
t_mole_local mole
Definition: mole.cpp:8
float realnum
Definition: cddefines.h:124
multi_arr< realnum, 2 > H2_X_Hmin_back
Definition: h2_priv.h:519
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
realnum gas_phase
Definition: deuterium.h:22
multi_arr< realnum, 3 > H2_X_grain_formation_distribution
Definition: h2_priv.h:540
realnum gas_phase[LIMELM]
Definition: dense.h:76
#define ASSERT(exp)
Definition: cddefines.h:617
qList states
Definition: h2_priv.h:429
double den
Definition: mole.h:421
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
vector< GrainBin * > bin
Definition: grainvar.h:585
diatomics hd("hd", 4100.,&hmi.HD_total, Yan_H2_CS)
t_deuterium deut
Definition: deuterium.cpp:7
double eden
Definition: dense.h:201
valarray< long > nRot_hi[N_ELEC]
Definition: h2_priv.h:477
double alogte
Definition: phycon.h:92
#define fixit(a)
Definition: cddefines.h:416
GrainVar gv
Definition: grainvar.cpp:5
t_hmi hmi
Definition: hmi.cpp:5
const int ipHYDROGEN
Definition: cddefines.h:348
long int nVib_hi[N_ELEC]
Definition: h2_priv.h:475
double rate_h2_form_grains_used_total
Definition: grainvar.h:576
Definition: collision.h:19
double H2star_forms_hminus
Definition: hmi.h:164