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