cloudy trunk
Loading...
Searching...
No Matches
mole_reactions.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#include "cdstd.h"
4#include <stdlib.h>
5#include "cddefines.h"
6#include "colden.h"
7#include "physconst.h"
8#include "mole.h"
9#include "version.h"
10#include "mole_priv.h"
11#include "hmi.h"
12#include "conv.h"
13#include "dense.h"
14#include "thermal.h"
15#include "gammas.h"
16#include "grainvar.h"
17#include "ionbal.h"
18#include "opacity.h"
19#include "path.h"
20#include "radius.h"
21#include "thermal.h"
22#include "atmdat.h"
23#include "taulines.h"
24#include "trace.h"
25#include "deuterium.h"
26
27/*
28 * HOWTO:- add a reaction to the new CO network, as at 2006 December 18.
29 *
30 * add a line of the form
31 * O,C2=>CO,C:hmrate:SCALE:B:C # Data source
32 * to the relevant data file for the new reaction. The first component is
33 * the chemical reaction, the second is a string naming the function which
34 * is used to evaluate the rate coefficients.
35 *
36 * SCALE is an overall constant by which the reaction rate is scaled,
37 * the remaining arguments constants used by this function.
38 *
39 * If all the species have previously been defined, and the parser in
40 * mole_co_etc.c can understand the reaction string, then that's it
41 * (unless the rate function is of a different form to those already
42 * defined).
43 *
44 * The sources and sinks to other networks will need to be defined,
45 * as this hasn't been made generic yet.
46 *
47 */
48
49enum {UDFA = 0}; /* UDFA = 1 for: make UDFA comparison and stop */
50
51STATIC void newreact(const char label[], const char fun[], double a, double b, double c);
52STATIC long parse_reaction( count_ptr<mole_reaction>& rate, const char label[] );
53STATIC string canonicalize_reaction_label( const char label[] );
56/* run through all reactions and report which ones do not have a reverse reaction */
58STATIC double mole_get_equilibrium_constant( const char buf[] );
59STATIC double mole_partition_function( const molecule* const sp);
60STATIC void mole_generate_isotopologue_reactions( string atom_old, string atom_new );
61
62/* Save PBM of network coupling matrix fill */
63STATIC void plot_sparsity(void);
64
65/* Check invariants for defined reaction (species and charge totals) */
67
68/* Functions to read UDFA database and compare results with Cloudy's internal reaction set */
69STATIC void read_data(const char file[], void (*parse)(char *s));
70STATIC void parse_base(char *s);
71STATIC void parse_udfa(char *s);
73
74/* Nick Abel, 06 Nov 27 states:
75 "In all the crnu reactions, the factor of two is the albedo factor 1/(1-albedo)."
76 Can this be obtained from the grain physics? */
77static realnum albedo = 0.5;
78
79
80namespace {
81/* Define ratefunc type to make specification of newfunc and findfunc clearer */
82 typedef count_ptr<mole_reaction> ratefunc;
83 template<class T> void newfunc()
84 {
85 ratefunc fun = count_ptr<mole_reaction>(new T);
86 ASSERT( mole_priv::functab.find(fun->name()) == mole_priv::functab.end() );
87 mole_priv::functab[fun->name()] = fun;
88 }
89 ratefunc findfunc(const char name[]);
90
91/* The functions which define the reaction rates -- must all have same template */
92 class mole_reaction_hmrate;
93 double hmrate(const mole_reaction *rate);
94 class mole_reaction_th85rate;
95 double th85rate(const mole_reaction *rate);
96 class mole_reaction_crnurate_noalbedo;
97 double crnurate_noalbedo(const mole_reaction *rate);
98/* >>chng 07 Dec 11, add photodesorption of molecules frozen on grain surfaces */
99/* >>chng 07 Dec 11, add grain surface reactions to chemistry. Now two molecules adsorbed on a grain can interact to form a new molecule */
100 class mole_reaction_grn_react;
101 double grn_react(const mole_reaction *rate);
102 class mole_reaction_h2_collh_excit;
103 double h2_collh_excit(const mole_reaction *rate);
104 class mole_reaction_h2_collh2_excit;
105 double h2_collh2_excit(const mole_reaction *rate);
106 class mole_reaction_h2_collh_deexc;
107 double h2_collh_deexc(const mole_reaction *rate);
108 class mole_reaction_h2_collh2_deexc;
109 double h2_collh2_deexc(const mole_reaction *rate);
110 class mole_reaction_rh2g_dis_h;
111 double rh2g_dis_h(const mole_reaction *rate);
112 class mole_reaction_rh2s_dis_h;
113 double rh2s_dis_h(const mole_reaction *rate);
114 class mole_reaction_rh2g_dis_h2;
115 double rh2g_dis_h2(const mole_reaction *rate);
116 class mole_reaction_rh2s_dis_h2;
117 double rh2s_dis_h2(const mole_reaction *rate);
118 class mole_reaction_rh2s_dis_h2_nodeexcit;
119 double rh2s_dis_h2_nodeexcit(const mole_reaction *rate);
120 class mole_reaction_bh2g_dis_h;
121 double bh2g_dis_h(const mole_reaction *rate);
122 class mole_reaction_bh2s_dis_h;
123 double bh2s_dis_h(const mole_reaction *rate);
124 class mole_reaction_bh2g_dis_h2;
125 double bh2g_dis_h2(const mole_reaction *rate);
126 class mole_reaction_hneut;
127 double hneut(const mole_reaction *rate);
128 class mole_reaction_cionhm;
129 double cionhm(const mole_reaction *rate);
130}
131
132/*
133 * Functions to specify chemical rates -- note that the rate->a overall scale
134 * factor is applied in mole_update_rks
135 *
136 */
137
138#include "phycon.h"
139#include "doppvel.h"
140
141namespace
142{
143 class mole_reaction_null : public mole_reaction
144 {
145 typedef mole_reaction_null T;
146 public:
147 virtual T* Create() const {return new T;}
148 virtual const char* name() {return "null";}
149 double rk() const
150 {
151 ASSERT( false );
152 return 0.0;
153 }
154 };
155}
156
157
158/* Add in non-equilibrium chemistry. This is done by assuming
159 * that turbulence reduces the temperature barrier, thereby
160 * enhancing reaction rates for molecules such as CH+. The
161 * "effective temperature is defined according to
162 * >>refer Federman et al. 1996, MNRAS. L41-46 */
163
164namespace {
165/* The effective temperature is defined as:
166 * T(effective) = T + (1/3)*reduced_mass*turbulent_velocity^2/BOLTZMANN_CONSTANT
167 */
168 STATIC double noneq_offset(const mole_reaction *rate)
169 {
170 /* This logic could be cached by using distinct rate functions in newreact */
171 int nreact, n;
172 bool lgFact;
173
174 DEBUG_ENTRY("noneq_offset()");
175
176 lgFact = false;
177 if(mole_global.lgNonEquilChem)
178 {
179 if(mole_global.lgNeutrals)
180 {
181 lgFact = true;
182 }
183 else
184 {
185 nreact = rate->nreactants;
186 for(n=0;n<nreact;n++)
187 {
188 if (rate->reactants[n]->charge != 0)
189 {
190 lgFact = true;
191 break;
192 }
193 }
194 }
195 }
196
197 if( lgFact )
198 return 0.333f*POW2(DoppVel.TurbVel)/BOLTZMANN*rate->reduced_mass;
199 else
200 return 0.;
201 }
202 double hmrate(const mole_reaction *rate)
203 {
204 double te;
205
206 DEBUG_ENTRY("hmrate()");
207
208 te = phycon.te+noneq_offset(rate);
209
210 if( rate->c < 0. )
211 ASSERT( -rate->c/te < 10. );
212
213 return pow(te/300.,rate->b)*exp(-rate->c/te);
214 }
215
216 class mole_reaction_hmrate_exo : public mole_reaction
217 {
218 typedef mole_reaction_hmrate_exo T;
219 public:
220 virtual T* Create() const {return new T;}
221
222 virtual const char* name() {return "hmrate_exo";}
223
224 double rk() const
225 {
226 double te;
227
228 DEBUG_ENTRY("hmrate_exo()");
229
230 te = phycon.te+noneq_offset(this);
231
232 if( this->c < 0. )
233 te = MIN2( te, -10. * this->c );
234
235 return pow(te/300.,this->b)*exp(-te/this->c);
236 }
237
238 };
239
240 class mole_reaction_hmrate : public mole_reaction
241 {
242 typedef mole_reaction_hmrate T;
243 public:
244 virtual T* Create() const {return new T;}
245 virtual const char* name() {return "hmrate";}
246 double rk() const
247 {
248 return hmrate(this);
249 }
250 };
251
252 class mole_reaction_constrate : public mole_reaction
253 {
254 typedef mole_reaction_constrate T;
255 public:
256 virtual T* Create() const {return new T;}
257 virtual const char* name() {return "constrate";}
258 double rk() const
259 {
260 return 1.;
261 }
262 };
263
264 /* hmi.UV_Cont_rel2_Habing_TH85_depth is field relative to Habing background, dimensionless */
265 /* >>chng 04 apr 01, move from TH85 to DB96, also correct leading coef to use
266 * UMIST database value */
267 /* CO_photo_dissoc_rate = 1.67e-10f*hmi.UV_Cont_rel2_Habing_TH85_depth;*/
268
269 /* TRY MOVING PHOTORATES OUT OF LOOP */
270
271 /* >>chng 02 jul 04 -- The following are dissociation rates for various molecular species
272 For right now we will calculate this rate by the standard form:
273 (alpha)*Go*exp(-Beta*AV)
274 when the command "set Leiden hack UMIST rates" is used. Otherwise we
275 will just let cloudy calculate the value of the UV radiation field */
276}
277#include "rfield.h"
278namespace {
279 double th85rate(const mole_reaction *rate)
280 {
281 double rk;
282
283 DEBUG_ENTRY("th85rate()");
284
285 if (!mole_global.lgLeidenHack || rate->c == 0.0)
286 {
287 rk = hmi.UV_Cont_rel2_Habing_TH85_depth/1.66;
288 }
289 else
290 {
291 rk = hmi.UV_Cont_rel2_Habing_TH85_face/1.66*exp(-(rate->c*rfield.extin_mag_V_point));
292 }
293
294 return rk;
295 }
296 class mole_reaction_th85rate : public mole_reaction
297 {
298 typedef mole_reaction_th85rate T;
299 public:
300 virtual T* Create() const {return new T;}
301 virtual const char* name() {return "th85rate";}
302 double rk() const
303 {
304 return th85rate(this);
305 }
306 };
307}
308
309#include "secondaries.h"
310namespace {
311 /* >> chng aug 24, 05 NPA This is the cosmic ray ionization rate used in the molecular network.
312 * TH85 and the e-mail from Amiel Sternberg has each cosmic ray ionization rate as a
313 * leading coefficient multiplied by a scale factor.
314 * The leading coefficient is defined as the cosmic ray rate for H + CRPHOT = > H+
315 * + e- . For molecules in the heavy element molecular network, this scale
316 * factor is derived by taking the rate for:
317
318 X + CRPHOT => Y + Z
319 and dividing it by the rate:
320 H + CRPHOTP => H+ + e-
321
322 This scale factor is 2.17 for all cosmic ray reactions in this network
323 crnu_rate = secondaries.csupra[ipHYDROGEN][0];*/
324 double crnurate_noalbedo(const mole_reaction *)
325 {
326 return 2.17*secondaries.csupra[ipHYDROGEN][0];
327 }
328 class mole_reaction_crnurate_noalbedo : public mole_reaction
329 {
330 typedef mole_reaction_crnurate_noalbedo T;
331 public:
332 virtual T* Create() const {return new T;}
333
334 virtual const char* name() {return "crnurate_noalbedo";}
335 double rk() const
336 {
337 return crnurate_noalbedo(this);
338 }
339 };
340
341 class mole_reaction_crnurate : public mole_reaction
342 {
343 typedef mole_reaction_crnurate T;
344 public:
345 virtual T* Create() const {return new T;}
346 virtual const char* name() {return "crnurate";}
347 double rk() const
348 {
349 return crnurate_noalbedo(this)/(1.0-albedo);
350 }
351 };
352}
353#include "hextra.h"
354namespace {
355/* Can this be found directly from radiation field?? */
356 class mole_reaction_cryden_ov_bg : public mole_reaction
357 {
358 typedef mole_reaction_cryden_ov_bg T;
359 public:
360 virtual T* Create() const {return new T;}
361 virtual const char* name() {return "cryden_ov_bg";}
362 double rk() const
363 {
365 }
366 };
367
368 class mole_reaction_co_lnu_c_o_lnu : public mole_reaction
369 {
370 typedef mole_reaction_co_lnu_c_o_lnu T;
371 public:
372 virtual T* Create() const {return new T;}
373 virtual const char* name() {return "co_lnu_c_o_lnu";}
374 double rk() const
375 {
376 double val = 0;
377 int ns, ion;
378 /* inner shell photoionization of CO, assume rates are same as K 1s and 2s
379 * shell of C and O */
380 /* >>chng 04 may 26, upper limit should be ns<2, had been ns<2 so picked up
381 * valence shell, which was incorrect */
382
383 DEBUG_ENTRY("co_lnu_c_o_lnu()");
384
385 for( ns=0; ns<2; ++ns )
386 {
387 ion = 0;
388 val += ionbal.PhotoRate_Shell[ipCARBON][ion][ns][0];
389 val += ionbal.PhotoRate_Shell[ipOXYGEN][ion][ns][0];
390 }
391
392 return val;
393 }
394 };
395
396 /******************************** Gas-Grain Chemistry**********************************/
397
398 /* The Gas-Grain Chemistry rates are taken from:
399 >>refer Hasegawa, T. I. & Herbst, E. 1993, MNRAS, 261, 83
400
401 So far only CO depletion onto grains is considered, however, this code can be generalized
402 if desired to other molecules, using the data in the above paper. There are three important reactions
403 to determine the abundance of a molecule on grain surfaces deep in molecular clouds. The
404 rate of accretion of molecule X onto grains is
405
406 R(accretion) = PI*[grain_radius]^2*[thermal velocity]*[density_of_grains]
407
408 Two processes remove molecule X from the grain surface, and that is thermal evaporation, due
409 to the heat of the grain, and cosmic ray deabsorption. The first of these rates come from the
410 above paper, and depends primarily on the dust temperature. The cosmic ray rate is a constant,
411 calculated in Hasegawa and Herbst.
412
413 For each molecule desired, I have created a new species which is the density of that molecule
414 on the grain surface */
415
416
417 /* evaporation rate of molecule on grain is:
418
419 k(evap) = [vibrational absorption frequency]*exp[-binding_energy/dust_temperature]
420
421 The binding energies come from Hasegawa and Herbst, Table 4. The vibrational frequency comes from
422 equation 3 of >>refer Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJSS, 82, 167
423
424 [vibrational absorption frequency] =
425 SQRT[[2*number_of_sites_for_grain_absorption*binding_energy]/[PI^2*mass_of_molecule]]
426
427 **********************************************************************************************/
428
429 class mole_reaction_vib_evap : public mole_reaction
430 {
431 typedef mole_reaction_vib_evap T;
432 public:
433 virtual T* Create() const {return new T;}
434 virtual const char* name() {return "vib_evap";}
435 double rk() const
436 {
437 double binding_energy, exponent, vib_freq;
438
439 DEBUG_ENTRY("vib_evap()");
440
441 exponent = 0.0;
442
443 binding_energy = this->b;
444 double bin_total=0.0;
445 for( size_t nd=0; nd < gv.bin.size() ; nd++ )
446 {
447 double bin_area = gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
448 exponent += exp(-binding_energy/gv.bin[nd]->tedust)*bin_area;
449 bin_total += bin_area;
450 }
451 exponent /= bin_total;
452 const double surface_density_of_sites = 1.5e15;
453
454 /* >> chng 07 dec 08 rjrw -- add 0.3*BOLTZMANN factor from Hasegawa & Herbst */
455 vib_freq = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*binding_energy/(PI*PI*this->reactants[0]->mole_mass));
456 //fprintf(ioQQQ,"Vib freq %g for mass %g\n",vib_freq,rate->reactants[0]->mole_mass);
457
458 /*>>chng 06 jan 11, NPA - In some H+ regions the grain temperatures are so low
459 that molecular freeze out occurs. This should not happen, because the ices
460 should sublimate in such a harsh environment. Therefore, we introduce an
461 artificial sublimation rate to destroy ices. THIS IS NOT A PHYSICAL RATE!!!!
462 only a rate that gets the desired, realistic result */
463 /*>>chng 06 sep 03 rjrw -- include this in standard evaporation rate coeff (the artificial part is the sexp term) */
465 /* Rate comes from Table curve and assumes that rate is high (~1) in H+
466 * region and negligible ( < 1e-30) in molecular cloud - designed to stop
467 * freeze-out above 100 K */
468
469 return vib_freq*exponent +sexp( 555.89/phycon.sqrte - 5.55 );
470 }
471 };
472
473 class mole_reaction_grn_abs : public mole_reaction
474 {
475 typedef mole_reaction_grn_abs T;
476 public:
477 virtual T* Create() const {return new T;}
478 virtual const char* name() {return "grn_abs";}
479 double rk() const
480 {
481 double mass;
482
483 DEBUG_ENTRY("grn_abs()");
484
485 /* Can't rely on order, as it is standardized using mole_cmp */
486 if (this->reactants[0]->n_nuclei() != 0)
487 mass = this->reactants[0]->mole_mass;
488 else
489 mass = this->reactants[1]->mole_mass;
490
491 return sqrt(8.*BOLTZMANN*phycon.te/(PI*mass));
492 }
493 };
494
495 class mole_reaction_grn_react : public mole_reaction
496 {
497 typedef mole_reaction_grn_react T;
498 public:
499 virtual T* Create() const {return new T;}
500 virtual const char* name() {return "grn_react";}
501 double rk() const
502 {
503 return grn_react(this);
504 }
505 };
506 double grn_react(const mole_reaction *rate)
507 {
508 /* This function computes the formation of a new molecule on the surface of
509 * of a grain due to the interaction between two species on the grain surface.
510 * We already do this type of reaction for H + H --> H2, but now we are doing
511 * it for other species as well. The rate depends on:
512
513 1) The ability of each species to move on the surface
514 2) The ability of each species to interact and form a new molecule
515 3) The ability of each species to not evaporate away before interacting
516
517 * The variable "number_of_sites" tells us how many sites there are which
518 * a molecule can attach itself to a grain. quant_barrier tells us the size
519 * of the quantum mechanical barrier between the two atoms. E_i and E_j
520 * relates to how well a species can stay on a grain of temperature T(dust). Finally,
521 * activ_barrier provides information on how well two atoms on the surface of a grain
522 * can actually interact to form a new molecule. The formalism below is taken from
523 * the Hasegawa, Herbst, and Leung paper from 1992 (referenced above)*/
524
525 DEBUG_ENTRY("grn_react()");
526
527 double quant_barrier = 1e-8; /* 1 Angstrom rectangular barrier */
528 double surface_density_of_sites = 1.5e15; /* number of sites available for adsorption */
529 fixit(); // rate->a is _always_ overall scaling factor
530 ASSERT( rate->nreactants == 2 ); // this routine seems coded to only work with 2 reactants
531 double E_i = rate->reactants[0]->form_enthalpy; /* adsorption energy barrier (w/ grain) for first reactant */
532 double E_j = rate->reactants[1]->form_enthalpy; /* adsorption energy barrier (w/ grain)for second reactant */
533 double activ_barrier = rate->c; /* energy barrier between atoms on the grain */
534
535 /* Each species has a vibrational frequency, dependent on its adsorption energy barrier */
536 fixit(); /* Ordering of reactants may have switched */
537 double vib_freq_i = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_i/(PI*PI*rate->reactants[0]->mole_mass));
538 double vib_freq_j = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_j/(PI*PI*rate->reactants[1]->mole_mass));
539
540 double Exp_i = 0.0;
541 double Exp_j = 0.0;
542 double dust_density = 0.0;
543
544 /* Compute thermal evaporation terms along with total dust number density */
545 fixit(); // should there be an area weighting to exponential terms?
546
547 for( size_t nd=0; nd < gv.bin.size(); nd++ )
548 {
549 double bin_density = gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
550 Exp_i += exp(-E_i/gv.bin[nd]->tedust)*bin_density;
551 Exp_j += exp(-E_j/gv.bin[nd]->tedust)*bin_density;
552 dust_density += bin_density/(4*1e-10);
553 }
554
555 ASSERT(fp_equal((realnum)dust_density,(realnum)(mole.grain_area/1e-10)));
556
557 double total_sites = 4.*mole.grain_area*surface_density_of_sites;
558
559 /* rates of migration on surface for each species */
560 double diff_i = vib_freq_i*Exp_i/total_sites;
561 double diff_j = vib_freq_j*Exp_j/total_sites;
562
563 /* Exponential term models the likelyhood that two species, once they meet on the surface, will interact. Larger the barrier, the smaller
564 * the chance of a new molecule */
565 double scale = exp(-2*(quant_barrier/(HBAReV*EN1EV))*sqrt(2.0*rate->reduced_mass*0.3*BOLTZMANN*activ_barrier));
566
567 /* Total Rate */
568 return scale*(diff_i + diff_j)/SDIV(dust_density);
569 }
570
571 class mole_reaction_grn_photo : public mole_reaction
572 {
573 typedef mole_reaction_grn_photo T;
574 public:
575 virtual T* Create() const {return new T;}
576 virtual const char* name() {return "grn_photo";}
577 double rk() const
578 {
579 /* This function models the rate at which UV photons remove
580 * molecules adsorbed on the surface of grains, using the treatment given in:
581 * >>refer Bergin, E. A., Langer, W. D., & Goldsmith, P. F. 1995, ApJ, 441, 222
582 * The following treatment uses their equation 3 */
583
584 DEBUG_ENTRY("grn_photo()");
585
586 /* This is the number of molecules removed from the grain per
587 * incident photon use same value of 1e-4 as used in Bergin, include
588 * conversion to radiation field units to get the dimensions right
589 * (cf the d'Hendercourt reference in Bergin et al) */
590
591 fixit();
592 // 2e-15 is mean adsorbed species cross section, from
593 // >>refer Greenberg, L. T., 1973, IAUS, 52, 413
594 // 1.232e7f * 1.71f is from DB96 referred to below
595 // will be multiplied by yield factor (~1e-4) to get overall photodesorption rate
596 return 2e-15 * hmi.UV_Cont_rel2_Draine_DB96_depth *(1.232e7f * 1.71f);
597 }
598 };
599}
600#include "rt.h"
601namespace {
602 class mole_reaction_th85rate_co : public mole_reaction
603 {
604 typedef mole_reaction_th85rate_co T;
605 public:
606 virtual T* Create() const {return new T;}
607
608 virtual const char* name() {return "th85rate_co";}
609
610 double rk() const
611 {
612 double esc_co, column;
613 /******************************************************************************************
614 * First define the rate, which is of the form:
615 *
616 * R = (Ro)*(Go*exp(-3.2Av))*Beta(tau(CO))
617 *
618 * where:
619 *
620 * Ro = 1.67*e-10
621 * (Go*exp(-3.2Av)) = hmi.UV_Cont_rel2_Habing_TH85_depth
622 * tauCO = 4.4e-15 * findspecies("CO")->column / (DopplerWidth/1e5) /
623 * (1. + phycon.sqrte*0.6019);
624 * tauC = 1.6*e17*N(C)
625 * Beta(tau(CO)) = esca0k2(esc_co)
626 ********************************************************************************************/
627 /* eqn 12 of
628 * >>refer CO dissoc Hollenbach, D.J., Takahashi, T., & Tielens, A. 1991, ApJ, 377, 192
629 * based on
630 * >>refer CO dissoc Black, J.H., & van Dishoeck, E.F. 1988, ApJ, 334, 771 */
631
632 if (this->reactants[0]->n_nuclei() != 0)
633 column = mole.species[ this->reactants[0]->index ].column;
634 else
635 column = mole.species[ this->reactants[1]->index ].column;
636
637 esc_co = 4.4e-15 * column /
638 /* the line width in km/s */
640 /* this term accounts for populations within ground elec state */
641 (1. + phycon.sqrte*0.6019);
642 return esca0k2(esc_co)*th85rate(this);
643 }
644 };
645
646 class mole_reaction_oh_c2h2_co_ch3 : public mole_reaction
647 {
648 typedef mole_reaction_oh_c2h2_co_ch3 T;
649 public:
650 virtual T* Create() const {return new T;}
651 virtual const char* name() {return "oh_c2h2_co_ch3";}
652 double rk() const
653 {
654 /* This rate will blow up if the temperature gets too low, therefore
655 * set this rate for T < 500 equal to the rate at 500 K */
656 if(phycon.te > 500)
657 {
658 return hmrate(this);
659 }
660 else
661 {
662 return 6.3E-18;
663 }
664 }
665 };
666
667
668 class mole_reaction_h_hnc_hcn_h : public mole_reaction
669 {
670 typedef mole_reaction_h_hnc_hcn_h T;
671 public:
672 virtual T* Create() const {return new T;}
673
674 virtual const char* name() {return "h_hnc_hcn_h";}
675 double rk() const
676 {
677 if(phycon.te > 100)
678 {
679 return hmrate(this);
680 }
681 else
682 {
683 return 1e-15;
684 }
685 }
686 };
687}
688#include "iso.h"
689#include "h2.h"
691{
692 /* >>chng 03 sep 09, get ratio of excited to ground state H2 */
693 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
694 {
695 return hmi.H2star_forms_hminus /
696 SDIV(hmi.H2star_forms_hminus+hmi.H2_forms_hminus);
697
698 /* option print statement for above */
699 /* printf( " H2s frac h- %.3e f(H2g) %.3e\n",frac_H2star_hminus ,
700 hmi.H2_forms_hminus/SDIV(hmi.H2star_forms_hminus+hmi.H2_forms_hminus));*/
701 }
702 else
703 {
704 /* the large H2 molecule was not evaluated, so we can't use exact
705 * branching ratios. These are the distribution fractions for around 500K */
706 /*These depend on temperature and distribution function and the definition of ENERGY_H2_STAR.
707 So reset the values properly*/
708 /* >>chng 05 jul 13, TE, with the new definition of H2s these are both 1 */
709 /* >>chng 05 jul 31, activate above print, rest for current 0.5 ev defn */
710 return 1. - 4.938e-6;
711 }
712}
713
714namespace {
715 double frac_H2star_grains(void)
716 {
717 /* >>chng 03 sep 09, get ratio of excited to ground state H2 */
718 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
719 {
720 return hmi.H2star_forms_grains /
721 SDIV(hmi.H2star_forms_grains+hmi.H2_forms_grains);
722
723 /* option print statement for above */
724 /*printf( "DEBUG H2s frac grain %.3e f(H2g) %.3e ",frac_H2star_grains ,
725 hmi.H2_forms_grains/SDIV(hmi.H2star_forms_grains+hmi.H2_forms_grains) ); */
726 }
727 else
728 {
729 /* the large H2 molecule was not evaluated, so we can't use exact
730 * branching ratios. These are the distribution fractions for around 500K */
731 /*These depend on temperature and distribution function and the definition of ENERGY_H2_STAR.
732 So reset the values properly*/
733 /* >>chng 05 jul 13, TE, with the new definition of H2s these are both 1 */
734 /* >>chng 05 jul 31, activate above print, rest for current 0.5 ev defn */
735 return 0.9416;
736 }
737 }
738
739 class mole_reaction_h2ph3p : public mole_reaction
740 {
741 typedef mole_reaction_h2ph3p T;
742 public:
743 virtual T* Create() const {return new T;}
744 virtual const char* name() {return "h2ph3p";}
745 double rk() const
746 {
747 /* H2 + H2+ => H + H3+ */
750
752 return 1.40e-9*(1. - sexp(9940./phycon.te));
753 else
754 return 2.08e-9;
755 }
756 };
757}
758
759#include "opacity.h"
760#include "gammas.h"
761#include "atmdat.h"
762
763namespace {
764 class mole_reaction_hopexch : public mole_reaction
765 {
766 typedef mole_reaction_hopexch T;
767 public:
768 virtual T* Create() const {return new T;}
769 virtual const char* name() {return "hopexch";}
770 double rk() const
771 {
773 }
774 };
775
776 class mole_reaction_hpoexch : public mole_reaction
777 {
778 typedef mole_reaction_hpoexch T;
779 public:
780 virtual T* Create() const {return new T;}
781 virtual const char* name() {return "hpoexch";}
782 double rk() const
783 {
785 }
786 };
787
788 class mole_reaction_hmattach : public mole_reaction
789 {
790 typedef mole_reaction_hmattach T;
791 public:
792 virtual T* Create() const {return new T;}
793 virtual const char* name() {return "hmattach";}
794 double rk() const
795 {
797 }
798 };
799
800 class mole_reaction_hmihph2p : public mole_reaction
801 {
802 typedef mole_reaction_hmihph2p T;
803 public:
804 virtual T* Create() const {return new T;}
805 virtual const char* name() {return "hmihph2p";}
806 double rk() const
807 {
808 /* H- + H+ => H2+ + e
809 * equation (H6) from
810 * >>refer H2 chemistry Galli,D., & Palla, F. 1998, A&A, 335, 403-420
811 * hmihph2p = 6.9e-9f*(Tg)^(-0.35) for Tg<=8000
812 * hmihph2p = 6.9e-9f*(Tg)^(-0.9) for Tg>=8000 */
813 if(phycon.te <= 7891.)
814 {
815 /*hmihph2p = 6.9e-9*pow(phycon.te , -0.35);*/
816 return 6.9e-9 / (phycon.te30 * phycon.te05);
817 }
818 else
819 {
820 /* >>chng 02 nov 18, had typo for leading coefficient here */
821 /*hmihph2p = 9.6e-7*pow(phycon.te , -0.9);*/
822 return 9.6e-7 / phycon.te90;
823 }
824 }
825 };
826
827 class mole_reaction_hmphoto : public mole_reaction
828 {
829 typedef mole_reaction_hmphoto T;
830 public:
831 virtual T* Create() const {return new T;}
832 virtual const char* name() {return "hmphoto";}
833 double rk() const
834 {
835 return hmi.HMinus_photo_rate;
836 }
837 };
838
839 double cionhm(const mole_reaction *)
840 {
841 /* collisional ionization of H-, rate from Janev, Langer et al. */
842 if( phycon.te < 3074. )
843 {
844 return 1.46e-32*(powi(phycon.te,6))*phycon.sqrte*hmi.exphmi;
845 }
846 else if( phycon.te >= 3074. && phycon.te < 30000. )
847 {
848
849 /* >>chng 03 mar 07, from above to below */
850 /*cionhm = 5.9e-19*phycon.te*phycon.te*phycon.sqrte*phycon.te03*
851 phycon.te01*phycon.te01;*/
852 return 5.9e-19*phycon.tesqrd*phycon.sqrte*phycon.te05;
853 }
854 else
855 {
856 return 1.54e-7;
857 }
858
859 }
860 class mole_reaction_cionhm : public mole_reaction
861 {
862 typedef mole_reaction_cionhm T;
863 public:
864 virtual T* Create() const {return new T;}
865 virtual const char* name() {return "cionhm";}
866 double rk() const
867 {
868 return cionhm(this);
869 }
870 };
871
872 double assoc_detach(void)
873 {
874 /* form molecular hydrogen from H minus,
875 * associative detachment: H- + H => H2 + E */
876 /* make H2 from H-
877 * associative detachment; H- + H => H2:
878 * >>referold H2 rates Browne & Dalgarno J PHys B 2, 885 */
879 /* rate coefficient from
880 * >>refer H2 form Launay, J.R., Le Dourneuf, M., & Zeippen, C.J.,
881 * >>refercon 1991, A&A, 252, 842-852*/
882 /* >>chng 02 oct 17, temp dependent fit to rate, updated reference,
883 * about 40% larger than before */
884 double y , x;
885 x = MAX2(10., phycon.te );
886 x = MIN2(1e4, x );
887 y=545969508.1323510+x*71239.23653059864;
888 return 1./y;
889 }
890
891 class mole_reaction_c3bod : public mole_reaction
892 {
893 typedef mole_reaction_c3bod T;
894 public:
895 virtual T* Create() const {return new T;}
896 virtual const char* name() {return "c3bod";}
897 double rk() const
898 {
899 return cionhm(this)*hmi.rel_pop_LTE_Hmin;
900 }
901 };
902
903 class mole_reaction_asdfg : public mole_reaction
904 {
905 typedef mole_reaction_asdfg T;
906 public:
907 virtual T* Create() const {return new T;}
908 virtual const char* name() {return "asdfg";}
909 double rk() const
910 {
911 return assoc_detach()*(1-frac_H2star_hminus());
912 }
913 };
914
915 class mole_reaction_asdfs : public mole_reaction
916 {
917 typedef mole_reaction_asdfs T;
918 public:
919 virtual T* Create() const {return new T;}
920 virtual const char* name() {return "asdfs";}
921 double rk() const
922 {
923 return assoc_detach()*frac_H2star_hminus();
924 }
925 };
926
927 class mole_reaction_asdbg : public mole_reaction
928 {
929 typedef mole_reaction_asdbg T;
930 public:
931 virtual T* Create() const {return new T;}
932 virtual const char* name() {return "asdbg";}
933 double rk() const
934 {
935 double ratio = mole_get_equilibrium_constant("H-,H=>H2,e-");
936 return assoc_detach() * ratio * (1.-frac_H2star_hminus());
937 }
938 };
939
940 class mole_reaction_asdbs : public mole_reaction
941 {
942 typedef mole_reaction_asdbs T;
943 public:
944 virtual T* Create() const {return new T;}
945 virtual const char* name() {return "asdbs";}
946 double rk() const
947 {
948 double ratio = mole_get_equilibrium_constant("H-,H=>H2*,e-");
949 return assoc_detach() * ratio * frac_H2star_hminus();
950 }
951 };
952
953 class mole_reaction_bhneut : public mole_reaction
954 {
955 typedef mole_reaction_bhneut T;
956 public:
957 virtual T* Create() const {return new T;}
958 virtual const char* name() {return "bhneut";}
959 double rk() const
960 {
961 if( phycon.te > 1000. && dense.xIonDense[ipHYDROGEN][0] > 0.0 )
962 {
963 double ratio = mole_get_equilibrium_constant("H-,H+=>H,H");
964 /* HBN(3,1) is defined; when <HydTempLimit then set to 1 */
965 /* mutual neut, mostly into n=3; rates from Janev et al
966 * H + H(n=3) => H- + H+ */
968 /* this is the back reaction, forming H- from Ho */
969 return (hneut(this)*ratio*
970 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH3s].Pop()+
971 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH3p].Pop()+
972 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH3d].Pop()) /
974 }
975 else
976 {
977 return 0.;
978 }
979 }
980 };
981
982 double hneut(const mole_reaction *)
983 {
984 if( phycon.te < 14125. )
985 {
986 /* the fit in Lepp et al. explodes at high temperature,
987 * Te = 14,125 is the temp where the rates reaches its lowest value */
988 return 1.4e-7*pow(phycon.te/300,-0.487)*exp(phycon.te/29300);
989 }
990 else
991 {
992 return 3.4738192887404660e-008;
993 }
994 }
995 class mole_reaction_hneut : public mole_reaction
996 {
997 typedef mole_reaction_hneut T;
998 public:
999 virtual T* Create() const {return new T;}
1000 virtual const char* name() {return "hneut";}
1001
1002 double rk() const
1003 {
1004 return hneut(this);
1005 }
1006 };
1007
1008 class mole_reaction_h2_spon_diss : public mole_reaction
1009 {
1010 typedef mole_reaction_h2_spon_diss T;
1011 public:
1012 virtual T* Create() const {return new T;}
1013 virtual const char* name() {return "h2_spon_diss";}
1014 double rk() const
1015 {
1016 return h2.spon_diss_tot;
1017 }
1018 };
1019
1020 double grnh2tot(const mole_reaction *)
1021 {
1022 /* Should remove factor of dense.gas_phase[ipHYDROGEN] factor from
1023 * derivation of rate_h2_form_grains_used to avoid any need for
1024 * division here */
1025 fixit();
1026 if( mole.grain_area*dense.xIonDense[ipHYDROGEN][0]>0 )
1027 return gv.rate_h2_form_grains_used_total/(mole.grain_area*dense.xIonDense[ipHYDROGEN][0]);
1028 else
1029 return 0.;
1030 }
1031 class mole_reaction_grnh2tot : public mole_reaction
1032 {
1033 public:
1034 virtual const char* name() {return "grnh2tot";}
1035 double rk() const
1036 {
1037 return grnh2tot(this);
1038 }
1039 };
1040
1041 class mole_reaction_grnh2 : public mole_reaction
1042 {
1043 typedef mole_reaction_grnh2 T;
1044 public:
1045 virtual T* Create() const {return new T;}
1046 virtual const char* name() {return "grnh2";}
1047 double rk() const
1048 {
1049 return grnh2tot(this)*(1.-frac_H2star_grains());
1050 }
1051 };
1052
1053 class mole_reaction_grnh2s : public mole_reaction
1054 {
1055 typedef mole_reaction_grnh2s T;
1056 public:
1057 virtual T* Create() const {return new T;}
1058 virtual const char* name() {return "grnh2s";}
1059 double rk() const
1060 {
1061 return grnh2tot(this)*frac_H2star_grains();
1062 }
1063 };
1064
1065 class mole_reaction_radasc : public mole_reaction
1066 {
1067 typedef mole_reaction_radasc T;
1068 public:
1069 virtual T* Create() const {return new T;}
1070 virtual const char* name() {return "radasc";}
1071 double rk() const
1072 {
1073 // excited atom radiative association,
1074 // H(n=2) + H(n=1) => H2 + hnu
1075
1076 if( dense.xIonDense[ipHYDROGEN][0] > 0. )
1077 {
1078 return hmrate(this) *
1080 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop() + iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2s].Pop())/
1082 }
1083 else
1084 return 0.;
1085 }
1086 };
1087
1088 class mole_reaction_assoc_ion : public mole_reaction
1089 {
1090 typedef mole_reaction_assoc_ion T;
1091 public:
1092 virtual T* Create() const {return new T;}
1093 virtual const char* name() {return "assoc_ion";}
1094 double rk() const
1095 {
1096 // excited atom associative ionization,
1097 // H(n=2) + H(n=1) => H2+ + e-
1098 if( dense.xIonDense[ipHYDROGEN][0] > 0. )
1099 {
1100 return hmrate(this) *
1102 (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop() + iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2s].Pop())/
1104 }
1105 else
1106 return 0.;
1107 }
1108 };
1109
1110
1111 double rh2g_dis_h(const mole_reaction *)
1112 {
1113 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1114 {
1115 return h2.Average_collH_dissoc_g;
1116 }
1117 else
1118 {
1119 double corr = MIN2(6.,14.44-phycon.alogte*3.08);
1120
1121 if(corr > 0.)
1122 corr = pow(10.,corr*findspecieslocal("H")->den/(findspecieslocal("H")->den+1.6e4));
1123 else
1124 corr = 1.;
1125 /* must kill H- when UMIST is in place since they do not consider it */
1126 return 1.55e-8/phycon.sqrte*sexp(65107./phycon.te)* corr;
1127 }
1128 }
1129 class mole_reaction_rh2g_dis_h : public mole_reaction
1130 {
1131 typedef mole_reaction_rh2g_dis_h T;
1132 public:
1133 virtual T* Create() const {return new T;}
1134 virtual const char* name() {return "rh2g_dis_h";}
1135 double rk() const
1136 {
1137 return rh2g_dis_h(this);
1138 }
1139 };
1140
1141 double rh2s_dis_h(const mole_reaction *rate)
1142 {
1143 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1144 {
1145 return h2.Average_collH_dissoc_s;
1146 }
1147 else
1148 {
1149 ASSERT( fp_equal( rate->a, 1. ) );
1150 return HMRATE(4.67e-7,-1.,5.5e4);
1151 }
1152 }
1153 class mole_reaction_rh2s_dis_h : public mole_reaction
1154 {
1155 typedef mole_reaction_rh2s_dis_h T;
1156 public:
1157 virtual T* Create() const {return new T;}
1158 virtual const char* name() {return "rh2s_dis_h";}
1159 double rk() const
1160 {
1161 return rh2s_dis_h(this);
1162 }
1163 };
1164
1165 double rh2g_dis_h2(const mole_reaction *rate)
1166 {
1167 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1168 {
1169 return h2.Average_collH2_dissoc_g;
1170 }
1171 else
1172 {
1173 /* >>refer H2 chemistry Palla, F., Salpeter, E.E., & Stahler, S.W., 1983, ApJ,271, 632-641 + detailed balance relation */
1174 ASSERT( fp_equal( rate->a, 1. ) );
1175 return HMRATE(5.5e-29*0.5/(SAHA*3.634e-5)*sqrt(300.),0.5,5.195e4);
1176 }
1177 }
1178 class mole_reaction_rh2g_dis_h2 : public mole_reaction
1179 {
1180 typedef mole_reaction_rh2g_dis_h2 T;
1181 public:
1182 virtual T* Create() const {return new T;}
1183 virtual const char* name() {return "rh2g_dis_h2";}
1184 double rk() const
1185 {
1186 return rh2g_dis_h2(this);
1187 }
1188 };
1189
1190 double rh2s_dis_h2(const mole_reaction *rate)
1191 {
1192 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1193 {
1194 return h2.Average_collH2_dissoc_s;
1195 }
1196 else
1197 {
1198 ASSERT( fp_equal( rate->a, 1. ) );
1199 return HMRATE(1e-11,0.,0.);
1200 }
1201 }
1202 class mole_reaction_rh2s_dis_h2 : public mole_reaction
1203 {
1204 typedef mole_reaction_rh2s_dis_h2 T;
1205 public:
1206 virtual T* Create() const {return new T;}
1207 virtual const char* name() {return "rh2s_dis_h2";}
1208 double rk() const
1209 {
1210 return rh2s_dis_h2(this);
1211 }
1212 };
1213
1214 double rh2s_dis_h2_nodeexcit(const mole_reaction *rate)
1215 {
1216 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1217 {
1218 return h2.Average_collH2_dissoc_s;
1219 }
1220 else
1221 {
1222 ASSERT( fp_equal( rate->a, 1. ) );
1223 return HMRATE(1e-11,0.,2.18e4);
1224 }
1225 }
1226 class mole_reaction_rh2s_dis_h2_nodeexcit : public mole_reaction
1227 {
1228 typedef mole_reaction_rh2s_dis_h2_nodeexcit T;
1229 public:
1230 virtual T* Create() const {return new T;}
1231 virtual const char* name() {return "rh2s_dis_h2_nodeexcit";}
1232 double rk() const
1233 {
1234 return rh2s_dis_h2_nodeexcit(this);
1235 }
1236 };
1237
1238 double bh2g_dis_h(const mole_reaction *rate)
1239 {
1240 return rh2g_dis_h(rate)*hmi.rel_pop_LTE_H2g;
1241 }
1242 class mole_reaction_bh2g_dis_h : public mole_reaction
1243 {
1244 typedef mole_reaction_bh2g_dis_h T;
1245 public:
1246 virtual T* Create() const {return new T;}
1247 virtual const char* name() {return "bh2g_dis_h";}
1248 double rk() const
1249 {
1250 return bh2g_dis_h(this);
1251 }
1252 };
1253
1254 double bh2s_dis_h(const mole_reaction *rate)
1255 {
1256 return rh2s_dis_h(rate)*hmi.rel_pop_LTE_H2s;
1257 }
1258 class mole_reaction_bh2s_dis_h : public mole_reaction
1259 {
1260 typedef mole_reaction_bh2s_dis_h T;
1261 public:
1262 virtual T* Create() const {return new T;}
1263 virtual const char* name() {return "bh2s_dis_h";}
1264 double rk() const
1265 {
1266 return bh2s_dis_h(this);
1267 }
1268 };
1269
1270 double bh2g_dis_h2(const mole_reaction *rate)
1271 {
1272 return rh2g_dis_h2(rate)*hmi.rel_pop_LTE_H2g;
1273 }
1274 class mole_reaction_bh2g_dis_h2 : public mole_reaction
1275 {
1276 typedef mole_reaction_bh2g_dis_h2 T;
1277 public:
1278 virtual T* Create() const {return new T;}
1279 virtual const char* name() {return "bh2g_dis_h2";}
1280 double rk() const
1281 {
1282 return bh2g_dis_h2(this);
1283 }
1284 };
1285
1286 class mole_reaction_bh2s_dis_h2 : public mole_reaction
1287 {
1288 typedef mole_reaction_bh2s_dis_h2 T;
1289 public:
1290 virtual T* Create() const {return new T;}
1291 virtual const char* name() {return "bh2s_dis_h2";}
1292 double rk() const
1293 {
1294 return rh2s_dis_h2(this)*hmi.rel_pop_LTE_H2s;
1295 }
1296 };
1297
1298 class mole_reaction_h2photon : public mole_reaction
1299 {
1300 typedef mole_reaction_h2photon T;
1301 public:
1302 virtual T* Create() const {return new T;}
1303 virtual const char* name() {return "h2photon";}
1304 double rk() const
1305 {
1306 return h2.photoionize_rate;
1307 }
1308 };
1309
1310 class mole_reaction_h2crphot : public mole_reaction
1311 {
1312 typedef mole_reaction_h2crphot T;
1313 public:
1314 virtual T* Create() const {return new T;}
1315 virtual const char* name() {return "h2crphot";}
1316 double rk() const
1317 {
1318 return secondaries.csupra[ipHYDROGEN][0]*2.02;
1319 }
1320 };
1321
1322 class mole_reaction_h2crphh : public mole_reaction
1323 {
1324 typedef mole_reaction_h2crphh T;
1325 public:
1326 virtual T* Create() const {return new T;}
1327 virtual const char* name() {return "h2crphh";}
1328 double rk() const
1329 {
1331 {
1332 return secondaries.x12tot*10.;
1333 }
1334 else
1335 {
1336 return secondaries.x12tot*3.;
1337 }
1338 }
1339 };
1340
1341 class mole_reaction_h2scrphh : public mole_reaction
1342 {
1343 typedef mole_reaction_h2scrphh T;
1344 public:
1345 virtual T* Create() const {return new T;}
1346 virtual const char* name() {return "h2scrphh";}
1347 double rk() const
1348 {
1350 {
1351 return secondaries.x12tot*10.;
1352 }
1353 else
1354 {
1355 return secondaries.x12tot*3.;
1356 }
1357 }
1358 };
1359
1360 class mole_reaction_radath : public mole_reaction
1361 {
1362 typedef mole_reaction_radath T;
1363 public:
1364 virtual T* Create() const {return new T;}
1365 virtual const char* name() {return "radath";}
1366 double rk() const
1367 {
1368 return MAX2(0.,2.325*MIN2(5000.,phycon.te)-1375.)*1e-20;
1369 }
1370 };
1371
1372 class mole_reaction_gamtwo : public mole_reaction
1373 {
1374 typedef mole_reaction_gamtwo T;
1375 public:
1376 virtual T* Create() const {return new T;}
1377 virtual const char* name() {return "gamtwo";}
1378 double rk() const
1379 {
1380 t_phoHeat dummy;
1381 return GammaK(opac.ih2pnt[0],opac.ih2pnt[1],opac.ih2pof,1.,&dummy);
1382 }
1383 };
1384
1385 class mole_reaction_hlike_phot : public mole_reaction
1386 {
1387 typedef mole_reaction_hlike_phot T;
1388 public:
1389 virtual T* Create() const {return new T;}
1390 virtual const char* name() {return "hlike_phot";}
1391 double rk() const
1392 {
1393 /* photoionization by hard photons, crossection =2*HI (wild guess)
1394 * -- rjrw: where do they go???
1395 * -- H3+ + hv => H2+ + H+ + e, best guess (P. Stancil, priv comm) */
1396
1397 // on first pass this has not been set, do it here, but only do it once.
1398 if( !conv.nTotalIoniz )
1400 return iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc;
1401 }
1402 };
1403
1404 class mole_reaction_h2s_sp_decay : public mole_reaction
1405 {
1406 typedef mole_reaction_h2s_sp_decay T;
1407 public:
1408 virtual T* Create() const {return new T;}
1409 virtual const char* name() {return "h2s_sp_decay";}
1410 double rk() const
1411 {
1412 /* >>chng 05 jul 11, TE, rename to use in punch file*/
1413 /* >>chng 05 jul 9, GS, use average A calculated from Big H2 */
1415 {
1416 return h2.Average_A;
1417 }
1418 else
1419 {
1420 return 2e-7;
1421 }
1422 }
1423 };
1424
1425 double h2_collh2_deexc(const mole_reaction *)
1426 {
1427 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1428 {
1429 return h2.Average_collH2_deexcit;
1430 }
1431 else
1432 {
1433 return ((1.4e-12*phycon.sqrte * sexp( 18100./(phycon.te + 1200.) ))/6.);
1434 }
1435 }
1436 class mole_reaction_h2_collh2_deexc : public mole_reaction
1437 {
1438 typedef mole_reaction_h2_collh2_deexc T;
1439 public:
1440 virtual T* Create() const {return new T;}
1441 virtual const char* name() {return "h2_collh2_deexc";}
1442 double rk() const
1443 {
1444 return h2_collh2_deexc(this);
1445 }
1446 };
1447
1448 double h2_collh_deexc(const mole_reaction *)
1449 {
1450 if ( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1451 {
1452 return h2.Average_collH_deexcit;
1453 }
1454 else
1455 {
1456 return ((1e-12*phycon.sqrte * sexp(1000./phycon.te ))/6.);
1457 }
1458 }
1459 class mole_reaction_h2_collh_deexc : public mole_reaction
1460 {
1461 typedef mole_reaction_h2_collh_deexc T;
1462 public:
1463 virtual T* Create() const {return new T;}
1464 virtual const char* name() {return "h2_collh_deexc";}
1465 double rk() const
1466 {
1467 return h2_collh_deexc(this);
1468 }
1469 };
1470
1471 double h2_collh2_excit(const mole_reaction *rate)
1472 {
1473 /* >>chng 05 jul 10, GS, use average collisional rate calculated from Big H2 */
1474 if ( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1475 {
1476 return h2.Average_collH2_excit;
1477 }
1478 else
1479 {
1480 return h2_collh2_deexc(rate)*sexp( 30172./phycon.te); /* deexc_htwo*Boltz_fac_H2_H2star */
1481 }
1482 }
1483 class mole_reaction_h2_collh2_excit : public mole_reaction
1484 {
1485 typedef mole_reaction_h2_collh2_excit T;
1486 public:
1487 virtual T* Create() const {return new T;}
1488 virtual const char* name() {return "h2_collh2_excit";}
1489 double rk() const
1490 {
1491 return h2_collh2_excit(this);
1492 }
1493 };
1494
1495 double h2_collh_excit(const mole_reaction *rate)
1496 {
1497 /* >>chng 05 jul 10, GS, use average collisional rate calculated from Big H2 */
1498 if ( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
1499 {
1500 return h2.Average_collH_excit;
1501 }
1502 else
1503 {
1504 return h2_collh_deexc(rate)*sexp( 30172./phycon.te); /* deexc_hneut*Boltz_fac_H2_H2star */
1505 }
1506 }
1507 class mole_reaction_h2_collh_excit : public mole_reaction
1508 {
1509 typedef mole_reaction_h2_collh_excit T;
1510 public:
1511 virtual T* Create() const {return new T;}
1512 virtual const char* name() {return "h2_collh_excit";}
1513 double rk() const
1514 {
1515 return h2_collh_excit(this);
1516 }
1517 };
1518
1519 class mole_reaction_h2gexcit : public mole_reaction
1520 {
1521 typedef mole_reaction_h2gexcit T;
1522 public:
1523 virtual T* Create() const {return new T;}
1524 virtual const char* name() {return "h2gexcit";}
1525 double rk() const
1526 {
1528 }
1529 };
1530
1531 class mole_reaction_h2sdissoc : public mole_reaction
1532 {
1533 typedef mole_reaction_h2sdissoc T;
1534 public:
1535 virtual T* Create() const {return new T;}
1536 virtual const char* name() {return "h2sdissoc";};
1537 double rk() const
1538 {
1539 /* >>chng 03 sep 11, add this process */
1540 /* photo-destroy H2* by Solomon process at same rate as H2ground dissociation,
1541 see above eqn A12 in TH85 */
1542 /* >>chng 00 nov 25 factor of 0.1, assume pump is total, and 10% destroy H2 */
1543 /* >>chng 03 mar 07, had factor of 0.1 for branching ratio from H2** to H+H,
1544 * but branching is now already included */
1546 }
1547 };
1548
1549 class mole_reaction_h2gdissoc : public mole_reaction
1550 {
1551 typedef mole_reaction_h2gdissoc T;
1552 public:
1553 virtual T* Create() const {return new T;}
1554 virtual const char* name() {return "h2gdissoc";}
1555 double rk() const
1556 {
1557 /* >>chng 03 sep 11, add this process */
1558 /* photo-destroy H2* by Solomon process at same rate as H2ground dissociation,
1559 see above eqn A12 in TH85 */
1560 /* >>chng 00 nov 25 factor of 0.1, assume pump is total, and 10% destroy H2 */
1561 /* >>chng 03 mar 07, had factor of 0.1 for branching ratio from H2** to H+H,
1562 * but branching is now already included */
1564 }
1565 };
1566
1567 class mole_reaction_hd_photodissoc : public mole_reaction
1568 {
1569 typedef mole_reaction_hd_photodissoc T;
1570 public:
1571 virtual T* Create() const {return new T;}
1572 virtual const char* name() {return "hd_photodissoc";}
1573 double rk() const
1574 {
1576 }
1577 };
1578}
1579
1580/*hmirat compute radiative association rate for H- */
1581double hmirat(double te)
1582{
1583 double hmirat_v;
1584
1585 DEBUG_ENTRY( "hmirat()" );
1586
1587 /* fits to radiative association rate coefficients */
1588 if( te < 31.62 )
1589 {
1590 hmirat_v = 8.934e-18*phycon.sqrte*phycon.te003*phycon.te001*
1591 phycon.te001;
1592 }
1593 else if( te < 90. )
1594 {
1595 hmirat_v = 5.159e-18*phycon.sqrte*phycon.te10*phycon.te03*
1596 phycon.te03*phycon.te003*phycon.te001;
1597 }
1598 else if( te < 1200. )
1599 {
1600 hmirat_v = 2.042e-18*te/phycon.te10/phycon.te03;
1601 }
1602 else if( te < 3800. )
1603 {
1604 hmirat_v = 8.861e-18*phycon.te70/phycon.te03/phycon.te01*
1605 phycon.te003;
1606 }
1607 else if( te <= 1.4e4 )
1608 {
1609 /* following really only optimized up to 1e4 */
1610 hmirat_v = 8.204e-17*phycon.sqrte/phycon.te10/phycon.te01*
1611 phycon.te003;
1612 }
1613 else
1614 {
1615 /* >>chng 00 sep 28, add this branch */
1616 hmirat_v = 5.44e-16*phycon.te20/phycon.te01;
1617 }
1618
1619 return( hmirat_v );
1620}
1621STATIC void mole_h2_grain_form(void);
1622
1624STATIC void mole_h_reactions(void);
1625
1626
1627namespace {
1628 class mole_reaction_gamheh : public mole_reaction
1629 {
1630 typedef mole_reaction_gamheh T;
1631 public:
1632 virtual T* Create() const {return new T;}
1633 virtual const char* name() {return "gamheh";}
1634 double rk() const
1635 {
1636 double retval = 0.;
1637 long int limit,
1638 i;
1639 /* photodissociation through 1.6->2.3 continuum */
1640
1641 /* why is this in a look instead of GammaK?
1642 * to fix must set opacities into stack */
1643
1644 limit = MIN2(hmi.iheh2-1 , rfield.nflux );
1645 for( i=hmi.iheh1-1; i < limit; i++ )
1646 {
1647 retval += rfield.flux[0][i] + rfield.ConInterOut[i]+ rfield.outlin[0][i] + rfield.outlin_noplot[i];
1648 }
1649 retval *= 4e-18;
1650
1651 /* hard radiation */
1652 retval += 3.*iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc;
1653
1654 return retval;
1655 }
1656 };
1657
1658 enum {exclude, base, umisthack, federman, lithium, deuterium, ti, misc};
1659 static int source;
1660
1661
1662 static bool lgReactInitialized = false;
1663}
1664
1666{
1667 /* Should adaptively select reactions rather than list them explicitly here */
1668
1669 /* prevent memory leaks */
1670 /* \todo this is a temporary fix for PR14. We should improve the overall design
1671 * of this code to prevent valid pointers being overwritten in a second call to mole_create_react */
1672 if( lgReactInitialized )
1673 return;
1674
1675 lgReactInitialized = true;
1676
1677 DEBUG_ENTRY("mole_create_react()");
1678
1679 if (UDFA)
1680 read_data("rate05_feb17.csv",parse_udfa);
1681
1682 /* Set up registry of rate functions -- could do something intelligent about
1683 * caching rate values by extending this API... */
1684 newfunc<mole_reaction_null>();
1685 newfunc<mole_reaction_hmrate>();
1686 newfunc<mole_reaction_hmrate_exo>();
1687 newfunc<mole_reaction_constrate>();
1688 newfunc<mole_reaction_th85rate>();
1689 newfunc<mole_reaction_crnurate>();
1690 newfunc<mole_reaction_crnurate_noalbedo>();
1691 newfunc<mole_reaction_co_lnu_c_o_lnu>();
1692 newfunc<mole_reaction_vib_evap>();
1693 newfunc<mole_reaction_th85rate_co>();
1694 newfunc<mole_reaction_grn_abs>();
1695 /* >chng 07 Dec 11, add grain surface chemical reactions */
1696 newfunc<mole_reaction_grn_react>();
1697 /* >>chng 07 Dec 10, add photodesorption of grain molecules */
1698 newfunc<mole_reaction_grn_photo>();
1699 newfunc<mole_reaction_oh_c2h2_co_ch3>();
1700 newfunc<mole_reaction_h_hnc_hcn_h>();
1701
1702 newfunc<mole_reaction_gamheh>();
1703 newfunc<mole_reaction_hd_photodissoc>();
1704 newfunc<mole_reaction_h2gdissoc>();
1705 newfunc<mole_reaction_h2sdissoc>();
1706 newfunc<mole_reaction_h2gexcit>();
1707 newfunc<mole_reaction_h2_collh_excit>();
1708 newfunc<mole_reaction_h2_collh2_excit>();
1709 newfunc<mole_reaction_h2_collh_deexc>();
1710 newfunc<mole_reaction_h2_collh2_deexc>();
1711 newfunc<mole_reaction_h2s_sp_decay>();
1712 newfunc<mole_reaction_hlike_phot>();
1713 newfunc<mole_reaction_gamtwo>();
1714 newfunc<mole_reaction_h2ph3p>();
1715 newfunc<mole_reaction_radath>();
1716 newfunc<mole_reaction_cryden_ov_bg>();
1717 newfunc<mole_reaction_h2scrphh>();
1718 newfunc<mole_reaction_h2crphh>();
1719 newfunc<mole_reaction_h2photon>();
1720 newfunc<mole_reaction_h2crphot>();
1721 newfunc<mole_reaction_rh2g_dis_h>();
1722 newfunc<mole_reaction_rh2s_dis_h>();
1723 newfunc<mole_reaction_rh2g_dis_h2>();
1724 newfunc<mole_reaction_rh2s_dis_h2>();
1725 newfunc<mole_reaction_rh2s_dis_h2_nodeexcit>();
1726 newfunc<mole_reaction_bh2g_dis_h>();
1727 newfunc<mole_reaction_bh2s_dis_h>();
1728 newfunc<mole_reaction_bh2g_dis_h2>();
1729 newfunc<mole_reaction_bh2s_dis_h2>();
1730 newfunc<mole_reaction_radasc>();
1731 newfunc<mole_reaction_assoc_ion>();
1732 newfunc<mole_reaction_grnh2>();
1733 newfunc<mole_reaction_grnh2s>();
1734 newfunc<mole_reaction_h2_spon_diss>();
1735 newfunc<mole_reaction_bhneut>();
1736 newfunc<mole_reaction_hneut>();
1737 newfunc<mole_reaction_asdbs>();
1738 newfunc<mole_reaction_asdbg>();
1739 newfunc<mole_reaction_asdfs>();
1740 newfunc<mole_reaction_asdfg>();
1741 newfunc<mole_reaction_c3bod>();
1742 newfunc<mole_reaction_cionhm>();
1743 newfunc<mole_reaction_hmphoto>();
1744 newfunc<mole_reaction_hmihph2p>();
1745 newfunc<mole_reaction_hmattach>();
1746 //newfunc<mole_reaction_hpoexch>();
1747 //newfunc<mole_reaction_hopexch>();
1748
1749 source = base;
1750 read_data("mole_co_base.dat",parse_base);
1751 if (mole_global.lgFederman)
1752 {
1753 source = federman;
1754 read_data("mole_co_federman.dat",parse_base);
1755 }
1756 if (!mole_global.lgLeidenHack)
1757 {
1758 source = umisthack;
1759 read_data("mole_co_umisthack.dat",parse_base);
1760 }
1761
1762 source = lithium;
1763 read_data("mole_lithium.dat",parse_base);
1764
1765 source = deuterium;
1766 read_data("mole_deuterium.dat",parse_base);
1767
1768#if 0
1769 source = ti;
1770 read_data("mole_ti.dat",parse_base);
1771#endif
1772
1773 source = misc;
1774 read_data("mole_misc.dat",parse_base);
1775
1776 /* Load null reaction to delete real rate from database */
1777 if (!mole_global.lgProtElim)
1778 {
1779 source = exclude;
1780 newreact("C+,OH=>CO,H+","hmrate",0.,0.,0.);
1781 }
1782
1783 newreact("H,H,grn=>H2,grn","grnh2",1.,0.,0.);
1784 newreact("H,H,grn=>H2*,grn","grnh2s",1.,0.,0.);
1785 newreact("H-,PHOTON=>H,e-","hmphoto",1.,0.,0.);
1786
1787 /* mutual neutralization with heavies, rate from Dalgarno and McCray
1788 * all charged ions contribute equally,
1789 * H- + A+ => H + A */
1790 fixit(); // this should be atom_list instead of unresolved_atom_list, but we have not defined ionized species of all isotopes yet!!!
1791 //for(vector<count_ptr<chem_atom> >::iterator atom=atom_list.begin(); atom != atom_list.end(); ++atom)
1792 for(ChemAtomList::iterator atom=unresolved_atom_list.begin(); atom != unresolved_atom_list.end(); ++atom)
1793 {
1794 long nelem = (*atom)->el->Z-1;
1795 if( nelem >= ipHELIUM && dense.lgElmtOn[nelem] )
1796 {
1797 char react[32];
1798 sprintf(react,"H-,%s+=>H,%s", (*atom)->label().c_str(), (*atom)->label().c_str() );
1799 newreact(react,"hmrate",4e-6/sqrt(300.),-0.5,0.);
1800 }
1801 }
1802
1803 newreact("H,e-=>H-,PHOTON","hmattach",1.,0.,0.);
1804 newreact("H-,H+=>H2+,e-","hmihph2p",1.,0.,0.);
1805 newreact("H-,e-=>H,e-,e-","cionhm",1.,0.,0.);
1806 newreact("H,e-,e-=>H-,e-","c3bod",1.,0.,0.);
1807 newreact("H,H-=>H2,e-","asdfg",1.,0.,0.);
1808 newreact("H,H-=>H2*,e-","asdfs",1.,0.,0.);
1809 newreact("H2,e-=>H,H-","asdbg",1.,0.,0.);
1810 newreact("H2*,e-=>H,H-","asdbs",1.,0.,0.);
1811 newreact("H-,H+=>H,H","hneut",1.,0.,0.);
1812 newreact("H,H=>H-,H+","bhneut",1.,0.,0.);
1813 newreact("H2*=>H,H","h2_spon_diss",1.,0.,0.);
1814 if (!mole_global.lgLeidenHack)
1815 {
1816 newreact("H,H=>H2,PHOTON","radasc",3e-14,0.,0.);
1817 newreact("H,H=>H2+,e-","assoc_ion",3.27e-10,-0.35,17829.); /* >>refer H2 rates Rawlings J.M.C, Drew J.E, Barlow M. J., 1993, MNRAS 265, 968 */
1818 newreact("H2,H=>H,H,H","rh2g_dis_h",1.,0.,0.);
1819 newreact("H2,H2=>H,H,H2","rh2g_dis_h2",1.,0.,0.);
1820 newreact("H,H,H=>H2,H","bh2g_dis_h",1.,0.,0.); /* back rate, three body recombination, 2H + S => H_2 + S */
1821 newreact("H,H,H2=>H2,H2","bh2g_dis_h2",1.,0.,0.);
1822 //newreact("H,H,H2=>H2,H2","hmrate",5.5e-29/(8*300.),-1.,0.); /* >>refer H2 chemistry Palla, F., Salpeter, E.E., & Stahler, S.W., 1983, ApJ,271, 632-641 */
1823 newreact("H2,PHOTON=>H2+,e-","h2photon",1.,0.,0.);
1824 newreact("H2,CRPHOT=>H2+,e-","h2crphot",1.,0.,0.);
1825 newreact("H2*,PHOTON=>H2+,e-","h2photon",1.,0.,0.);
1826 newreact("H2*,CRPHOT=>H2+,e-","h2crphot",1.,0.,0.);
1827 newreact("H2,CRPHOT=>H,H","h2crphh",1.,0.,0.); /* >>refer H2 k Millar, T.J., et.al, 1997,A&AS, 121, 139 */
1828 newreact("H2,CRPHOT=>H+,H-","cryden_ov_bg",3.9e-21,0.,0.); /* >>refer H2 k Millar, T.J., et.al, 1997,A&AS, 121, 139 */
1829 newreact("H2*,CRPHOT=>H+,H-","cryden_ov_bg",3.9e-21,0.,0.); /* >>refer H2 k Millar, T.J., et.al, 1997,A&AS, 121, 139 */
1830 newreact("H2,CRPHOT=>H+,H,e-","cryden_ov_bg",2.2e-19,0.,0.); /* >>refer H2 k Millar, T.J., et.al, 1997,A&AS, 121, 139 */
1831 newreact("H2*,CRPHOT=>H+,H,e-","cryden_ov_bg",2.2e-19,0.,0.); /* >>refer H2 k Millar, T.J., et.al, 1997,A&AS, 121, 139 */
1832 newreact("H+,H=>H2+,PHOTON","radath",1.,0.,0.);
1833 newreact("H2+,PHOTON=>H,H+","gamtwo",1.,0.,0.);
1834 newreact("H2+,CRPHOT=>H,H+","hlike_phot",1.,0.,0.);
1835 /* >> chng 05 aug 05, NPA comment. This reaction is not in UMIST, for the case
1836 of hard photons. Turn if off for the comparison. */
1837 newreact("H3+,CRPHOT=>H2+,H+,e-","hlike_phot",2.,0.,0.);
1838 newreact("H,H+,e-=>H2+,e-","hmrate",2e-7 * SAHA * (4./(2.*1.)) * 3.634e-5 * pow(300.,-1.50),-1.50,0.);
1839 }
1840 else
1841 {
1842 newreact("H2,CRPHOT=>H2+,e-","hmrate",4.4e-17,0.,0.);
1843 newreact("H2,CRPHOT=>H,H+,e-","hmrate",1e-19,0.,0.);
1844 newreact("H2*,CRPHOT=>H,H+,e-","hmrate",1e-19,0.,0.);
1845 newreact("H2*,CRPHOT=>H2+,e-","hmrate",4.4e-17,0.,0.);
1846 newreact("H2,CRPHOT=>H,H","hmrate",5e-18,0.,0.);
1847 newreact("e-,H3+=>H2,H","hmrate",2.5e-8,-0.3,0.);
1848 newreact("e-,H3+=>H,H,H","hmrate",7.5e-8,-0.3,0.);
1849 newreact("H+,H=>H2+,PHOTON","hmrate",5.3e-19,1.85,0);
1850 /* H2+ + e=> H + H;
1851 * equation (6,table 4) from
1852 * >>refer H2 l Maloney et.al, 1996,ApJ, 466, 561 */
1853 newreact("H2+,e-=>H,H","hmrate",1.6e-8,-0.43,0.);
1854 newreact("H2+,PHOTON=>H,H+","th85rate",5.7e-10,0.,1.9);
1855 }
1856
1857 newreact("H2*,CRPHOT=>H,H","h2scrphh",1.,0.,0.); /* >>refer H2 k Millar, T.J., et.al, 1997,A&AS, 121, 139 */
1858 newreact("H2,H2+=>H,H3+","h2ph3p",1.,0.,0.);
1859 newreact("H2*=>H2,PHOTON","h2s_sp_decay",1.,0.,0.);
1860 newreact("H2*,H2=>H2,H2","h2_collh2_deexc",1.,0.,0.);
1861 newreact("H2*,H=>H2,H","h2_collh_deexc",1.,0.,0.);
1862 newreact("H2,H2=>H2*,H2","h2_collh2_excit",1.,0.,0.);
1863 newreact("H2,H=>H2*,H","h2_collh_excit",1.,0.,0.);
1864
1865 // collisional dissociation of H2*
1866 newreact("H2*,H=>H,H,H","rh2s_dis_h",1.,0.,0.);
1867 newreact("H2*,H2=>H2,H,H","rh2s_dis_h2",1.,0.,0.);
1868 newreact("H2*,H2*=>H2,H,H","rh2s_dis_h2",1.,0.,0.);
1869 newreact("H2*,H2*=>H2*,H,H","rh2s_dis_h2_nodeexcit",1.,0.,0.);
1870
1871#if 0
1872 // more collision dissociation of H2
1873 newreact("H2,He=>He,H,H","rh2g_dis_h",1.,0.,0.);
1874 newreact("H2,H+=>H+,H,H","rh2g_dis_h",1.,0.,0.);
1875 newreact("H2,H3+=>H3+,H,H","rh2g_dis_h",1.,0.,0.);
1876 newreact("H2,e-=>e-,H,H","rh2g_dis_h",1.,0.,0.);
1877 newreact("H2*,He=>He,H,H","rh2s_dis_h",1.,0.,0.);
1878 newreact("H2*,H+=>H+,H,H","rh2s_dis_h",1.,0.,0.);
1879 newreact("H2*,H3+=>H3+,H,H","rh2s_dis_h",1.,0.,0.);
1880 newreact("H2*,e-=>e-,H,H","rh2s_dis_h",1.,0.,0.);
1881
1882 // back reactions
1883 newreact("He,H,H=>H2,He","bh2g_dis_h",1.,0.,0.);
1884 newreact("H+,H,H=>H2,H+","bh2g_dis_h",1.,0.,0.);
1885 newreact("H3+,H,H=>H2,H3+","bh2g_dis_h",1.,0.,0.);
1886 newreact("e-,H,H=>H2,e-","bh2g_dis_h",1.,0.,0.);
1887 newreact("H,H,H=>H2*,H","bh2s_dis_h",1.,0.,0.);
1888 newreact("H2,H,H=>H2*,H2","bh2s_dis_h2",1.,0.,0.);
1889 newreact("H2,H,H=>H2*,H2*","bh2s_dis_h2",1.,0.,0.);
1890 newreact("H2*,H,H=>H2*,H2*","bh2s_dis_h2",1.,0.,0.);
1891 newreact("He,H,H=>H2*,He","bh2s_dis_h",1.,0.,0.);
1892 newreact("H+,H,H=>H2*,H+","bh2s_dis_h",1.,0.,0.);
1893 newreact("H3+,H,H=>H2*,H3+","bh2s_dis_h",1.,0.,0.);
1894 newreact("e-,H,H=>H2*,e-","bh2s_dis_h",1.,0.,0.);
1895#endif
1896
1897 newreact("H2,PHOTON=>H2*","h2gexcit",1.,0.,0.);
1898 newreact("H2*,PHOTON=>H,H","h2sdissoc",1.,0.,0.);
1899 newreact("H2,PHOTON=>H,H","h2gdissoc",1.,0.,0.);
1900 newreact("HeH+,PHOTON=>H,He+","gamheh",1.,0.,0.);
1901
1902 if(0)
1903 {
1904 // grain surface reactions
1905 newreact("OHgrn,Hgrn=>H2Ogrn","grn_react",1.,0.,0.);
1906 }
1907
1908 if (UDFA)
1909 {
1910 fprintf(stderr,"Finished testing against UDFA database\n");
1911 exit(-1);
1912 }
1913
1914 if (0)
1915 plot_sparsity();
1916
1918
1919 if( deut.lgElmtOn )
1921
1922 long index = 0;
1923 for( mole_reaction_i it = mole_priv::reactab.begin(); it != mole_priv::reactab.end(); ++it, ++index )
1924 it->second->index = index;
1925
1926 mole.reaction_rks.resize( index );
1927 mole.old_zone = -1;
1928 memset( &mole.reaction_rks[0], 0, (unsigned long)(index)*sizeof(double) );
1929
1930 // label catalytic, intra-group, and excitition/deexcitation vectors for all reactions
1931 for( mole_reaction_i it = mole_priv::reactab.begin(); it != mole_priv::reactab.end(); ++it )
1932 register_reaction_vectors( it->second );
1933}
1934
1935STATIC void mole_generate_isotopologue_reactions( string atom_old, string atom_new )
1936{
1937 DEBUG_ENTRY("mole_generate_isotopologue_reactions()");
1938
1939 bool lgDebug = false;
1940
1941 // store new reaction strings (and pointer to reaction derived from) and add all to reactab map after following iterator
1942 vector<string> newReactionStrings;
1943 vector<mole_reaction*> oldReactions;
1944
1945 // iterate over all reactions and generate new strings by replacing atom_old with atom_new
1946 for( mole_reaction_ci it = mole_priv::reactab.begin(); it != mole_priv::reactab.end(); ++it )
1947 {
1948 bool lgFound = false;
1949
1950 for( long i=0; i<it->second->nreactants; ++i )
1951 {
1952 // search for atom_old among reactants
1953 for( nAtoms_i atom=it->second->reactants[i]->nAtom.begin(); atom != it->second->reactants[i]->nAtom.end(); ++atom )
1954 {
1955 if( atom->first->label() == atom_old )
1956 {
1957 lgFound = true;
1958 continue;
1959 }
1960 }
1961 if( lgFound )
1962 continue;
1963 }
1964
1965 if( !lgFound )
1966 continue;
1967
1968 if( lgDebug )
1969 fprintf( ioQQQ, "DEBUGGG mole_generate_isotopologue_reactions %s ..........\n", it->first.c_str() );
1970
1971 for( long i=0; i<it->second->nreactants; ++i )
1972 {
1973 // ignore reactants with no nucleons (e-, PHOTON, CRPHOT, ... )
1974 if( it->second->reactants[i]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
1975 continue;
1976
1977 vector<string> react_iso_labels;
1978 ChemAtomList atomsLeftToRight;
1979 vector< int > numAtoms;
1980 string embellishments;
1981 // parse reactant label
1982 bool lgParseOK = parse_species_label( it->second->reactants[i]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
1983 ASSERT( lgParseOK == true );
1984 // generate isotopologue labels
1986 atomsLeftToRight,
1987 numAtoms,
1988 atom_old,
1989 atom_new,
1990 embellishments,
1991 react_iso_labels );
1992 for( unsigned j=0; j<react_iso_labels.size(); ++j )
1993 {
1994 for( long k=0; k<it->second->nproducts; ++k )
1995 {
1996 // ignore products with no nucleons (e-, PHOTON, CRPHOT, ... )
1997 if( it->second->products[k]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
1998 continue;
1999
2000 vector<string> prod_iso_labels;
2001 atomsLeftToRight.clear();
2002 numAtoms.clear();
2003 embellishments.clear();
2004 // parse product label
2005 lgParseOK = parse_species_label( it->second->products[k]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
2006 ASSERT( lgParseOK == true );
2007 // generate isotopologue labels
2009 atomsLeftToRight,
2010 numAtoms,
2011 atom_old,
2012 atom_new,
2013 embellishments,
2014 prod_iso_labels );
2015 for( unsigned l=0; l<prod_iso_labels.size(); ++l )
2016 {
2017 string react_string;
2018 //generate new reaction string
2019 // first write reactants
2020 for( long i1=0; i1<i; ++i1 )
2021 {
2022 react_string += it->second->reactants[i1]->label;
2023 if( i1 != it->second->nreactants-1 )
2024 react_string += ",";
2025 }
2026 react_string += react_iso_labels[j];
2027 if( i != it->second->nreactants-1 )
2028 react_string += ",";
2029 for( long i2=i+1; i2<it->second->nreactants; ++i2 )
2030 {
2031 react_string += it->second->reactants[i2]->label;
2032 if( i2 != it->second->nreactants-1 )
2033 react_string += ",";
2034 }
2035
2036 react_string += "=>";
2037 // now write products
2038 for( long k1=0; k1<k; ++k1 )
2039 {
2040 react_string += it->second->products[k1]->label;
2041 if( k1 != it->second->nproducts-1 )
2042 react_string += ",";
2043 }
2044 react_string += prod_iso_labels[l];
2045 if( k != it->second->nproducts-1 )
2046 react_string += ",";
2047 for( long k2=k+1; k2<it->second->nproducts; ++k2 )
2048 {
2049 react_string += it->second->products[k2]->label;
2050 if( k2 != it->second->nproducts-1 )
2051 react_string += ",";
2052 }
2053
2054 if( lgDebug )
2055 fprintf( ioQQQ, "DEBUGGG mole_generate_isotopologue_reactions .................... %s\n",
2056 react_string.c_str() );
2057
2058 string canon_react_string = canonicalize_reaction_label( react_string.c_str() );
2059 // store new reaction string and pointer to old
2060 newReactionStrings.push_back( canon_react_string );
2061 oldReactions.push_back( it->second.get_ptr() );
2062 }
2063 }
2064 }
2065 }
2066 }
2067
2068 ASSERT( oldReactions.size() == newReactionStrings.size() );
2069
2070 // now declare new reactions
2071 vector<mole_reaction*>::const_iterator it2 = oldReactions.begin();
2072 for( vector<string>::const_iterator it1 = newReactionStrings.begin(); it1 != newReactionStrings.end(); ++it1, ++it2 )
2073 {
2074 fixit(); // make adjustments to a for mass?
2075
2076 // don't overwrite existing reaction with these new auto-generated details
2077 if( mole_priv::reactab.find( it1->c_str() ) == mole_priv::reactab.end() )
2078 newreact( it1->c_str(), (*it2)->name(), (*it2)->a, (*it2)->b, (*it2)->c );
2079 }
2080
2081 return;
2082}
2083
2085{
2086 DEBUG_ENTRY("mole_check_reverse_reactions()");
2087
2088 char chLabel[50], chLabelSave[50];
2089 int exists;
2090
2091 for(mole_reaction_i p=mole_priv::reactab.begin();
2092 p != mole_priv::reactab.end(); ++p)
2093 {
2094 mole_reaction_i q = p;
2095 strcpy( chLabel, p->second->label.c_str() );
2096 strcpy( chLabelSave, p->second->label.c_str() );
2097 char *str = chLabel;
2098 const char *delim = "=>";
2099 char *chNewProducts = strtok( str, delim );
2100 char *chNewReactants = strtok( NULL, delim );
2101 char chNewReaction[50];
2102
2103 strcpy( chNewReaction, chNewReactants );
2104 strcat( chNewReaction, "=>" );
2105 strcat( chNewReaction, chNewProducts );
2106
2107
2108 q = mole_priv::reactab.find(chNewReaction);
2109
2110 exists = (q != mole_priv::reactab.end());
2111 if ( !exists )
2112 {
2113 if( trace.lgTraceMole )
2114 {
2115 fprintf(ioQQQ,"Warning! No reverse reaction for %30s. ", p->second->label.c_str() );
2116 fprintf( ioQQQ,"\n" );
2117 }
2118
2119 fixit();
2120 // NB reverse reactions should be generated here
2121 }
2122 }
2123
2124 return;
2125}
2126
2127STATIC double mole_get_equilibrium_constant( const char buf[] )
2128{
2129 DEBUG_ENTRY("mole_get_equilibrium_constant()");
2130
2131 mole_reaction *rate = mole_findrate_s(buf);
2132
2133 if( !rate )
2134 return 0;
2135
2136 // multiply by reactant partition functions, then divide by products'
2137 double ln_result = 0.;
2138 for( long i=0; i<rate->nreactants; ++i)
2139 {
2140 double fac = mole_partition_function(rate->reactants[i]);
2141 if( fac==0. )
2142 return 0.;
2143 ln_result += log(fac);
2144 }
2145 for( long i=0; i<rate->nproducts; ++i)
2146 {
2147 double fac = mole_partition_function(rate->products[i]);
2148 if (fac <= 0.)
2149 return (double) BIGFLOAT;
2150 //ASSERT( fac > 0. );
2151 ln_result -= log(fac);
2152 }
2153 double result = exp( ln_result );
2154
2155 //fprintf( ioQQQ, "DEBUGGG equilibrium %20s\t%e\n", rate->label.c_str(), result );
2156 return MIN2((double) BIGFLOAT,result);
2157}
2158
2160{
2161 DEBUG_ENTRY("mole_partition_function()");
2162
2163 if( sp->label == "PHOTON" || sp->label == "CRPHOT" )
2164 {
2165 fixit(); // How can we adapt existing structures to have a photon energy or range available here?
2166 fixit(); // include 2hnu^3/c^2. Understand HNU3C2 macro!
2167 return 1.; //sexp( energy_ryd/phycon.te_ryd );
2168 }
2169 else if( sp->label == "CRP" || sp->label == "grn" )
2170 return 1.;
2171
2172 fixit(); // need to figure out stat weight for any given particle
2173 double q = 1.;
2174 // last factors convert kJ/mol to Kelvin/particle
2175 double deltaH = sp->form_enthalpy * (1e10/AVOGADRO/BOLTZMANN);
2176 ASSERT( sp->mole_mass > 0. );
2177 double part_fun = pow(phycon.te*sp->mole_mass/(HION_LTE_POP*ELECTRON_MASS), 1.5) * q * dsexp(deltaH/phycon.te);
2178 ASSERT( part_fun < BIGFLOAT );
2179
2180 return part_fun;
2181}
2182
2184{
2185 DEBUG_ENTRY("mole_cmp_num_in_out_reactions()");
2186
2187 vector<long> numForm, numDest;
2188 numForm.resize( mole_global.num_total );
2189 numDest.resize( mole_global.num_total );
2190
2191 for(mole_reaction_i p=mole_priv::reactab.begin(); p != mole_priv::reactab.end(); p++)
2192 {
2193 count_ptr<mole_reaction> rate = p->second;
2194 for( long i=0; i<rate->nreactants; ++i)
2195 {
2196 ++numDest[ rate->reactants[i]->index ];
2197 }
2198
2199 for( long i=0; i<rate->nproducts; ++i)
2200 {
2201 ++numForm[ rate->products[i]->index ];
2202 }
2203 }
2204
2205 for( unsigned i=0; i<numForm.size(); ++i )
2206 {
2207 if( numForm[i]==0 && numDest[i]==0 )
2208 continue;
2209 if( numForm[i]>1 && numDest[i]>1 )
2210 continue;
2211 if( mole_global.list[i]->isMonatomic() )
2212 continue;
2213 fprintf( ioQQQ, "DEBUGGG mole_cmp_num_in_out_reactions %*s: in %4li out %4li\n", CHARS_SPECIES, mole_global.list[i]->label.c_str(), numForm[i], numDest[i] );
2214 }
2215
2216 return;
2217}
2218
2219STATIC char *getcsvfield(char **s,char c);
2220STATIC void parse_base(char *s)
2221{
2222 char *label, *reactstr, *f;
2223 double a, b, c;
2224 label = getcsvfield(&s,':');
2225 reactstr = getcsvfield(&s,':');
2226 f = getcsvfield(&s,':');
2227 a = atof(f);
2228 f = getcsvfield(&s,':');
2229 b = atof(f);
2230 f = getcsvfield(&s,':');
2231 c = atof(f);
2232
2233 newreact(label,reactstr,a,b,c);
2234
2235}
2236
2237STATIC void newreact(const char label[], const char fun[], double a, double b, double c)
2238{
2239 DEBUG_ENTRY("newreact()");
2240
2241 ratefunc rate = findfunc(fun);
2242 if (rate.get_ptr() == NULL)
2243 {
2244 fprintf(stderr,"Rate function %s not known for reaction %s. Aborting. Sorry.\n",fun,label);
2246 }
2247
2248 rate->label = label;
2249 rate->a = a;
2250 rate->b = b;
2251 rate->c = c;
2252 //rate->rk = 0.0;
2253 rate->source = source;
2254
2255 rate->photon = 0;
2256
2257 if( parse_reaction( rate, label ) == 0 )
2258 return;
2259
2260 canonicalize_reaction( rate );
2261
2262 const char *rateLabelPtr = rate->label.c_str();
2263
2264 ASSERT(lgReactBalance(rate)); /* Verify rate conserves particles and charge */
2265
2266 rate->udfastate = ABSENT;
2267 if (UDFA)
2268 {
2269 compare_udfa(rate);
2270 if (rate->udfastate == ABSENT)
2271 {
2272 fprintf(stderr,"Reaction %s not in UDFA\n",rateLabelPtr);
2273 }
2274 }
2275
2276 /* >> chng 06 Oct 10 rjrw: use 1/(1/m1+1/m2) for reduced mass to prevent underflow */
2277 if (rate->nreactants == 2 && rate->reactants[0]->mole_mass!=0.0 && rate->reactants[1]->mole_mass!=0.0 )
2278 {
2279 rate->reduced_mass = 1./(1./rate->reactants[0]->mole_mass+1./rate->reactants[1]->mole_mass);
2280 }
2281 else
2282 {
2283 rate->reduced_mass = 0.;
2284 }
2285
2286 /* If every thing is OK, can register the reaction */
2287 mole_reaction_i p = mole_priv::reactab.find(rateLabelPtr);
2288 int exists = (p != mole_priv::reactab.end());
2289 // do not comment in release version
2290 if(exists && !t_version::Inst().lgReleaseBranch && !t_version::Inst().lgRelease )
2291 {
2292 /* Replace old rate */
2293 fprintf(ioQQQ,"Warning: duplicate reaction %s -- using new version\n",rateLabelPtr);
2294 }
2295 mole_priv::reactab[rateLabelPtr] = rate;
2296}
2297
2298STATIC long parse_reaction( count_ptr<mole_reaction>& rate, const char label[] )
2299{
2300 DEBUG_ENTRY("parse_reaction()");
2301
2302 for (int i=0;i<MAXREACTANTS;i++)
2303 {
2304 rate->reactants[i] = NULL;
2305 }
2306 rate->nreactants = 0;
2307
2308 for (int i=0;i<MAXPRODUCTS;i++)
2309 {
2310 rate->products[i] = NULL;
2311 }
2312 rate->nproducts = 0;
2313
2314 bool lgProd = false;
2315 string buf = "";
2316 for(int i=0;!i || label[i-1]!='\0';i++)
2317 {
2318 if(label[i] == ',' || label[i] == '=' || label[i] == '\0')
2319 {
2320 molecule *sp = findspecies(buf.c_str());
2321 if( sp == null_mole || !sp->isEnabled )
2322 {
2323 if( trace.lgTraceMole )
2324 fprintf(ioQQQ,"Mole_reactions: ignoring reaction %s (species %s not active)\n",label,buf.c_str());
2325 return 0;
2326 }
2327 buf = "";
2328 if(! lgProd)
2329 {
2330 if (rate->nreactants >= MAXREACTANTS)
2331 {
2332 fprintf(stderr,"Mole_reactions: Too many reactants in %s, only %d allowed\n",label,MAXREACTANTS);
2333 exit(-1);
2334 }
2335 rate->reactants[rate->nreactants] = sp;
2336 rate->nreactants++;
2337 }
2338 else
2339 {
2340 if (rate->nproducts >= MAXPRODUCTS)
2341 {
2342 fprintf(stderr,"Mole_reactions: Too many products in %s, only %d allowed\n",label,MAXPRODUCTS);
2343 exit(-1);
2344 }
2345 rate->products[rate->nproducts] = sp;
2346 rate->nproducts++;
2347 }
2348 if(label[i] == '=')
2349 {
2350 i++; /* skip '>' as well */
2351 if (label[i] != '>')
2352 {
2353 fprintf(ioQQQ,"Format error in reaction %s\n",label);
2355 }
2356 lgProd = true;
2357 }
2358 }
2359 else
2360 {
2361 buf += label[i];
2362 }
2363 }
2364
2365 ASSERT( rate->nreactants );
2366 ASSERT( rate->nproducts );
2367
2368 return 1;
2369}
2370
2371STATIC string canonicalize_reaction_label( const char label[] )
2372{
2373 DEBUG_ENTRY("canonicalize_reaction_label()");
2374
2375 // set up a dummy reaction
2376 ratefunc rate = findfunc("null");
2377 rate->label = label;
2378 parse_reaction( rate, label );
2379 canonicalize_reaction( rate );
2380
2381 //if( !conv.nTotalIoniz && strcmp( label, rate->label.c_str() ) != 0 )
2382 // fprintf( ioQQQ, "DEBUGGG reaction label %20s canonicalized to %20s\n", label, rate->label.c_str() );
2383
2384 return rate->label;
2385}
2386
2388{
2389 DEBUG_ENTRY("canonicalize_reaction()");
2390
2391 /* Put species in canonical order to make sure reactions are unique.
2392 Can cause problems when a consistent ordering of species is
2393 required (look for references to "reactants[0]" for examples) */
2394 t_mole_global::sort(rate->reactants,rate->reactants+rate->nreactants);
2395 t_mole_global::sort(rate->products,rate->products+rate->nproducts);
2396
2397 // now reorder label in same (new) order
2398 string newLabel;
2399 for( long i=0; i<rate->nreactants; ++i )
2400 {
2401 newLabel += rate->reactants[i]->label;
2402 if( i != rate->nreactants-1 )
2403 newLabel += ",";
2404 }
2405 newLabel += "=>";
2406 for( long i=0; i<rate->nproducts; ++i )
2407 {
2408 newLabel += rate->products[i]->label;
2409 if( i != rate->nproducts-1 )
2410 newLabel += ",";
2411 }
2412
2413 // overwrite original label with canonical
2414 rate->label = newLabel;
2415
2416 return;
2417}
2418
2420{
2421 DEBUG_ENTRY("register_reaction_vectors()");
2422
2423 for (long k=0;k<rate->nreactants;k++)
2424 {
2425 rate->rvector[k] = NULL;
2426 rate->rvector_excit[k] = NULL;
2427 }
2428
2429 for (long k=0;k<rate->nproducts;k++)
2430 {
2431 rate->pvector[k] = NULL;
2432 rate->pvector_excit[k] = NULL;
2433 }
2434
2435 /* Label catalytic species */
2436 for (long i=0;i<rate->nproducts;i++)
2437 {
2438 if (rate->pvector[i] == NULL)
2439 {
2440 for (long k=0;k<rate->nreactants;k++)
2441 {
2442 if (rate->rvector[k] == NULL)
2443 {
2444 if (rate->products[i] == rate->reactants[k])
2445 {
2446 rate->rvector[k] = rate->products[i];
2447 rate->pvector[i] = rate->reactants[k];
2448 break;
2449 }
2450 }
2451 }
2452 }
2453 }
2454
2455 /* Label other intra-group transfers */
2456 for (long i=0;i<rate->nproducts;i++)
2457 {
2458 if (rate->pvector[i] == NULL)
2459 {
2460 for (long k=0;k<rate->nreactants;k++)
2461 {
2462 if (rate->rvector[k] == NULL)
2463 {
2464 if (rate->products[i]->groupnum != -1 &&
2465 rate->products[i]->groupnum ==
2466 rate->reactants[k]->groupnum)
2467 {
2468 rate->rvector[k] = rate->products[i];
2469 rate->pvector[i] = rate->reactants[k];
2470 break;
2471 }
2472 }
2473 }
2474 }
2475 }
2476
2477 /* Label excited/deexcited pairs */
2478 for (long i=0;i<rate->nproducts;i++)
2479 {
2480 if (rate->pvector[i] == NULL && rate->pvector_excit[i] == NULL)
2481 {
2482 for (long k=0;k<rate->nreactants;k++)
2483 {
2484 if (rate->rvector[k] == NULL && rate->rvector_excit[k] == NULL )
2485 {
2486 if ( lgDifferByExcitation( *rate->products[i], *rate->reactants[k] ) )
2487 {
2488 rate->rvector_excit[k] = rate->products[i];
2489 rate->pvector_excit[i] = rate->reactants[k];
2490 break;
2491 }
2492 }
2493 }
2494 }
2495 }
2496
2497 return;
2498}
2499
2500
2502{
2503 FILE *sparsefp;
2504 int i, j, nb, ch;
2505 long int ratej;
2506 double **c;
2507
2508 c = (double **)MALLOC((size_t)mole_global.num_total*sizeof(double *));
2509 c[0] = (double *)MALLOC((size_t)mole_global.num_total*
2510 mole_global.num_total*sizeof(double));
2511
2512 for(i=1;i<mole_global.num_total;i++)
2513 {
2514 c[i] = c[i-1]+mole_global.num_total;
2515 }
2516
2517 for ( j=0; j < mole_global.num_total; j++ )
2518 {
2519 for ( i=0; i < mole_global.num_total; i++ )
2520 {
2521 c[j][i] = 0.;
2522 }
2523 }
2524
2525 for(mole_reaction_i p=mole_priv::reactab.begin();
2526 p != mole_priv::reactab.end(); ++p)
2527 {
2528 mole_reaction &rate = *p->second;
2529
2530 for (j=0;j<rate.nreactants;j++)
2531 {
2532 ratej = rate.reactants[j]->index;
2533 for (i=0;i<rate.nreactants;i++)
2534 {
2535 if (rate.rvector[i] == NULL)
2536 c[ratej][rate.reactants[i]->index] = 1.0;
2537 }
2538 for (i=0;i<rate.nproducts;i++)
2539 {
2540 if (rate.pvector[i] == NULL)
2541 c[ratej][rate.products[i]->index] = 1.0;
2542 }
2543 }
2544 }
2545
2546 sparsefp = fopen("sparse.pbm","w");
2547 fprintf(sparsefp,"P4\n%d %d\n",
2548 mole_global.num_total,mole_global.num_total);
2549
2550 for ( j=0; j < mole_global.num_total; j++ )
2551 {
2552 nb = ch = 0;
2553 for ( i=0; i < mole_global.num_total; i++ )
2554 {
2555 ch = (ch << 1) | (c[i][j] != 0.0);
2556 nb++;
2557 if (nb == 8)
2558 {
2559 fputc(ch,sparsefp);
2560 nb = ch = 0;
2561 }
2562 }
2563 if (nb != 0)
2564 {
2565 ch <<= 8-nb;
2566 fputc(ch,sparsefp);
2567 }
2568 }
2569 fclose(sparsefp);
2570 free(c[0]);
2571 free(c);
2572}
2573
2575{
2577 int dcharge, n, sign;
2578 bool lgOK = true;
2579
2580 dcharge = 0;
2581 for (n=0;n<rate->nreactants;n++)
2582 {
2583 for( nAtoms_i it = rate->reactants[n]->nAtom.begin(); it != rate->reactants[n]->nAtom.end(); ++it )
2584 nel[it->first] += it->second;
2585 dcharge += rate->reactants[n]->charge;
2586 }
2587 for (n=0;n<rate->nproducts;n++)
2588 {
2589 for( nAtoms_i it = rate->products[n]->nAtom.begin(); it != rate->products[n]->nAtom.end(); ++it )
2590 nel[it->first] -= it->second;
2591 dcharge -= rate->products[n]->charge;
2592 }
2593 if (dcharge != 0)
2594 {
2595 fprintf(stderr,"Reaction %s charge out of balance by %d\n",
2596 rate->label.c_str(),dcharge);
2597 lgOK = false;
2598 }
2599
2600 for( nAtoms_i it = nel.begin(); it != nel.end(); ++it )
2601 {
2602 if(it->second != 0)
2603 {
2604 if(it->second > 0)
2605 sign = 1;
2606 else
2607 sign = -1;
2608 fprintf(stderr,"Error: reaction %s %s %d of element %s\n",
2609 rate->label.c_str(),sign==1?"destroys":"creates",
2610 sign*it->second,
2611 it->first->label().c_str() );
2612 lgOK = false;
2613 }
2614 }
2615 return lgOK;
2616}
2617
2618#ifndef ZLIB_H
2619#define gzopen(file,mode) fopen(file,mode)
2620#define gzgets(fp,buf,siz) fgets(buf,siz,fp)
2621#define gzclose(fp) fclose(fp)
2622#define gzFile FILE
2623#endif
2624
2625enum {BUFSIZE=256};
2626
2627namespace
2628{
2629 class formula_species {
2630 public:
2631 molecule *reactants[MAXREACTANTS], *products[MAXPRODUCTS];
2632 };
2633
2634 bool operator< (const formula_species &a, const formula_species &b)
2635 {
2636 int i;
2637 for (i=0;i<MAXREACTANTS;i++)
2638 {
2639 if (a.reactants[i]<b.reactants[i])
2640 return true;
2641 else if (a.reactants[i]>b.reactants[i])
2642 return false;
2643 }
2644 for (i=0;i<MAXPRODUCTS;i++)
2645 {
2646 if (a.products[i]<b.products[i])
2647 return true;
2648 else if (a.products[i]>b.products[i])
2649 return false;
2650 }
2651 return false;
2652 }
2653
2654 struct udfa_reaction_s {
2655 int index;
2656 formula_species l;
2657 char source; /* Calculated, Estimated, Literature compilation, Measured */
2658 double a, b, c, trange[2]; /* Overall scale and valid temperature range */
2659 };
2660}
2661
2662static map <formula_species,struct udfa_reaction_s *> udfatab;
2663
2664STATIC void read_data(const char file[], void (*parse)(char *s))
2665{
2666 FILE *fp;
2667 char buf[BUFSIZE];
2668
2669 fp = open_data(file,"r");
2670 if (!fp)
2671 {
2672 fprintf(stderr,"Error, could not read %s\n",file);
2673 exit(-1);
2674 }
2675
2676 fixit(); // this seg-faults if file ends in blank line!
2677 while(fgets(buf,BUFSIZE,fp))
2678 {
2679 if( buf[0] == '#' )
2680 continue;
2681 parse(buf);
2682 }
2683 fclose(fp);
2684}
2685#define FLTEQ(a,b) (fabs((a)-(b)) <= 1e-6*fabs((a)+(b)))
2686STATIC void parse_udfa(char *s)
2687{
2688 char *f;
2689 struct udfa_reaction_s *r;
2690 map <formula_species, struct udfa_reaction_s *>::iterator p;
2691 unsigned int havespecies = 1, i, n;
2692 /* lgCRPHOT is true for CRPHOT reactions as we change the data
2693 * format for them */
2694 int lgCRPHOT=0;
2695
2696 r = (struct udfa_reaction_s *)MALLOC(sizeof(struct udfa_reaction_s));
2697 f = getcsvfield(&s,',');
2698 r->index = atoi(f);
2699
2700 /* Load reactants */
2701 for (n=0;n<MAXREACTANTS;n++)
2702 {
2703 r->l.reactants[n] = NULL;
2704 }
2705 i = 0;
2706 for (n=0;n<MAXREACTANTS;n++)
2707 {
2708 f = getcsvfield(&s,',');
2709 if (f[0] != '\0')
2710 {
2711 i++;
2712 r->l.reactants[n] = findspecies(f);
2713 if (r->l.reactants[n] == null_mole)
2714 havespecies = 0;
2715 if (!strncmp(f,"CRPHOT",6))
2716 lgCRPHOT = 1;
2717 }
2718 }
2719 t_mole_global::sort(r->l.reactants,r->l.reactants+i);
2720
2721 /* Load products */
2722 for (n=0;n<MAXPRODUCTS;n++)
2723 {
2724 r->l.products[n] = NULL; /* Sentinel */
2725 }
2726 i = 0;
2727 for (n=0;n<MAXPRODUCTS;n++)
2728 {
2729 f = getcsvfield(&s,',');
2730 if (f[0] != '\0')
2731 {
2732 i++;
2733 r->l.products[n] = findspecies(f);
2734 if (r->l.products[n] == null_mole)
2735 havespecies = 0;
2736 }
2737 }
2738
2739 t_mole_global::sort(r->l.products,r->l.products+i);
2740
2741 /* Load rate parameters */
2742 f = getcsvfield(&s,',');
2743 r->a = atof(f);
2744 f = getcsvfield(&s,',');
2745 r->b = atof(f);
2746 f = getcsvfield(&s,',');
2747 r->c = atof(f);
2748
2749 if (lgCRPHOT)
2750 {
2751 /* UDFA has a uniform value for the cosmic ray field independent of
2752 circumstances which we correct -- verify they haven't changed
2753 anything and move the multiplicative constant into our usual place
2754 for it. */
2755 ASSERT (FLTEQ(r->a,1.3e-17));
2756 r->a = r->c;
2757 r->c = 0.;
2758 }
2759
2760 /* Load data confidence and range of temperature validity */
2761 f = getcsvfield(&s,',');
2762 r->source = f[0]?f[0]:'?';
2763 for (n=0;n<2;n++) {
2764 f = getcsvfield(&s,',');
2765 r->trange[n] = atof(f);
2766 }
2767
2768 if (!havespecies)
2769 {
2770 free(r); /* We don't handle some of the species present (yet) */
2771 }
2772 else
2773 {
2774 p = udfatab.find(r->l);
2775 if (p != udfatab.end())
2776 {
2777 fprintf(stderr,"Duplicate reaction\n");
2778 }
2779 udfatab[r->l] = r;
2780 }
2781}
2782STATIC char *getcsvfield(char **s, char c)
2783{
2784 char *sv, *f;
2785
2786 sv = strchr(*s,c);
2787 if (sv) {
2788 *sv++ = '\0';
2789 }
2790 f = *s;
2791 *s = sv;
2792 return f;
2793}
2795{
2796 formula_species s;
2797 unsigned int n;
2798 map<formula_species, struct udfa_reaction_s *>::iterator p;
2799 struct udfa_reaction_s *u;
2800
2801 for (n=0;n<MAXREACTANTS;n++)
2802 {
2803 s.reactants[n] = rate->reactants[n];
2804 }
2805 for (n=0;n<MAXPRODUCTS;n++)
2806 {
2807 s.products[n] = rate->products[n];
2808 }
2809 p = udfatab.find(s);
2810 if (p == udfatab.end() )
2811 {
2812 /* fprintf(stdout,"Did not find reaction %s\n",rate->label); */
2813 return;
2814 }
2815 else
2816 {
2817 u = p->second;
2818 if (FLTEQ(rate->a,u->a) && FLTEQ(rate->b,u->b) && FLTEQ(rate->c,u->c))
2819 {
2820 rate->udfastate = CORRECT;
2821 /* fprintf(stdout,"Reaction %s agrees\n",rate->label); */
2822 }
2823 else
2824 {
2825 rate->udfastate = CONFLICT;
2826 /* fprintf(stdout,"Reaction %18.18s clashes: a %9.2e %9.2e|b %9.2e %9.2e|c %9.2e %9.2e\n",
2827 rate->label,rate->a,u->a,rate->b,u->b,rate->c,u->c); */
2828 }
2829 }
2830}
2831
2832namespace
2833{
2834 ratefunc findfunc (const char name[])
2835 {
2836 return count_ptr<mole_reaction>(mole_priv::functab[name]->Create());
2837 }
2838}
2839
2841{
2842 enum { DEBUG_MOLE = false };
2843
2844 DEBUG_ENTRY("mole_update_rks()");
2845
2847
2849
2850 for (mole_reaction_i p
2851 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
2852 {
2853 mole_reaction &rate = *p->second;
2854 long index = rate.index;
2855 realnum oldrk = (realnum)mole.reaction_rks[index];
2856 realnum newrk = rate.a*rate.rk();
2857 mole.reaction_rks[index] = newrk;
2858 if (DEBUG_MOLE)
2859 {
2860 if (fabs(newrk-oldrk) > 0.1*newrk)
2861 fprintf(ioQQQ,"%s: %15.8g => %15.8g\n",
2862 rate.label.c_str(),oldrk,newrk);
2863 }
2864 }
2865}
2867{
2868 enum { DEBUG_MOLE = false };
2869
2870 DEBUG_ENTRY("mole_rk_bigchange()");
2871
2872 if ( mole.old_reaction_rks.size() == 0 )
2873 {
2874 mole.old_zone = -1;
2875 mole.old_reaction_rks.resize(mole.reaction_rks.size());
2876 }
2877
2878 if (nzone > 1)
2879 {
2880 ASSERT(mole.old_zone == nzone - 1);
2881 double bigchange = 0.;
2882 unsigned long bigindex = ULONG_MAX;
2883 for (unsigned long index=0; index<mole.reaction_rks.size(); ++index)
2884 {
2885 double oldrk = mole.old_reaction_rks[index],
2886 newrk = mole.reaction_rks[index],
2887 sum = oldrk+newrk, diff = newrk-oldrk;
2888 if (sum > 0.)
2889 {
2890 double change = fabs(diff)/sum;
2891 if (change > bigchange)
2892 {
2893 bigchange = change;
2894 bigindex = index;
2895 }
2896 }
2897 }
2898
2899 for (mole_reaction_i p
2900 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
2901 {
2902 mole_reaction &rate = *p->second;
2903 if (bigindex == (unsigned) rate.index)
2904 {
2905 double oldrk = mole.old_reaction_rks[bigindex],
2906 newrk = mole.reaction_rks[bigindex];
2907 double pct = 0.;
2908 if (oldrk > 0.)
2909 pct = 100.*(newrk-oldrk)/oldrk;
2910 fprintf(ioQQQ, "Zone %ld, big chemistry rate change %s:"
2911 " %15.8g => %15.8g (%6.2g%%)\n",
2912 nzone,rate.label.c_str(),oldrk,newrk,pct);
2913 break;
2914 }
2915 }
2916 }
2917
2918 mole.old_zone = nzone;
2919 for (unsigned long index=0; index<mole.reaction_rks.size(); ++index)
2920 {
2921 mole.old_reaction_rks[index] = mole.reaction_rks[index];
2922 }
2923}
2924
2926{
2927 double T_ortho_para_crit;
2928 realnum AveVelH = GetAveVelocity( dense.AtomicWeight[ipHYDROGEN] );
2929 realnum AveVelH2 = GetAveVelocity( 2.f * dense.AtomicWeight[ipHYDROGEN] );
2930
2931 /* H2 formation on grains;
2932 * rate from
2933 * >>refer H2 grain formation Hollenbach, D., & McKee, C.F., 1979, ApJS, 41, 555 eq 3.4 3.8 */
2934 if( gv.lgDustOn() )
2935 {
2936
2937# ifndef IGNORE_QUANTUM_HEATING
2938 /* hmole is called before grains, so assure that all the grain stuff is properly initialized */
2939 GrainDrive();
2940# endif
2941
2942 /* these are rates (s-1) H2 will be deactivated by collisions with grains
2943 * will be incremented below
2944 * H2 ortho - para conversion on grain surface */
2945 h2.rate_grain_op_conserve = 0.;
2946 /* rate (s-1) v=0, J=1 level goes to 0 */
2947 h2.rate_grain_J1_to_J0 = 0.;
2948
2949 /* loop over all grain species */
2950 for( size_t nd=0; nd < gv.bin.size(); nd++ )
2951 {
2952# ifndef IGNORE_QUANTUM_HEATING
2953 long k, qnbin;
2954 double *qtemp, *qprob;
2955 bool lgUseQHeat = gv.lgGrainPhysicsOn && gv.bin[nd]->lgQHeat;
2956# endif
2957 /* >>chng 02 feb 15, removed check tedust > 1.01, change in GrainsInit
2958 * guarantees that all relevant parameters are initialized, PvH */
2959
2960 /* sticking probability, 2H + grain equation 3.7 of
2961 * >>refer grain phys Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555,
2962 * fraction of H impacts on grain surface that stick */
2963 /* this sticking probability is used for both HM79 and CT02 */
2964 double sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) +
2965 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
2966
2967# ifndef IGNORE_QUANTUM_HEATING
2968 /* >>chng 04 feb 21, included quantum heating in calculation of formation rate, PvH */
2969 if( lgUseQHeat )
2970 {
2971 qtemp = (double*)MALLOC((size_t)(NQGRID*sizeof(double)));
2972 qprob = (double*)MALLOC((size_t)(NQGRID*sizeof(double)));
2973
2974 qheat(qtemp,qprob,&qnbin,nd);
2975
2976 if( gv.bin[nd]->lgUseQHeat )
2977 {
2978 ASSERT( qnbin > 0 );
2979 }
2980 else
2981 {
2982 qnbin = 1;
2983 qprob[0] = 1.;
2984 qtemp[0] = gv.bin[nd]->tedust;
2985 }
2986
2987 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;
2988
2989 for( k=0; k < qnbin; k++ )
2990 {
2991 /* fraction of impacts that produce H2 before evaporation from grain surface.
2992 * this is equation 3.4 of
2993 * >>refer grain phys Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555
2994 * 1e4 is ratio of total absorption sites to appropriate sites
2995 * 920 is D_H and chosen to get f_a = 0.5 at 100 K.
2996 * factor of 0.6252 needed to obtain std ism rate to be 3e-17 at 100 K,
2997 * the value deduced by
2998 * >>refer H2 grain physics Jura, M., 1974, ApJ, 197, 581 */
2999 double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./qtemp[k]));
3000 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) +
3001 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
3002
3003 gv.bin[nd]->rate_h2_form_grains_HM79 += qprob[k] * sticking_probability_H *
3004 conversion_efficiency_HM79;
3005 }
3006
3007 /* NB IntArea is total, not projected, area, must div by 4 */
3008 /* gv.bin[nd]->rate_h2_form_grains_HM79 has units s^-1 since gv.bin[nd]->cnv_H_pCM3 has units cm-3 */
3009 /* cnv_H_pCM3 converts <unit>/H (default depletion) -> <unit>/cm^3 (actual depletion), units are cm-3 */
3010 gv.bin[nd]->rate_h2_form_grains_HM79 *= 0.5 * AveVelH *
3011 gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;
3012
3013 ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
3014 }
3015 else
3016# endif
3017 {
3018 /* fraction of impacts that produce H2 before evaporation from grain surface.
3019 * this is equation 3.4 of
3020 * >>refer grain phys Hollenbach, D.J., & McKee, C.F., 1979, ApJS, 41, 555
3021 * 1e4 is ratio of total absorption sites to appropriate sites
3022 * 920 is D_H and chosen to get f_a = 0.5 at 100 K.
3023 * factor of 0.6252 needed to obtain std ism rate to be 3e-17 at 100 K,
3024 * the value deduced by
3025 * >>refer H2 grain physics Jura, M., 1974, ApJ, 197, 581 */
3026 double conversion_efficiency_HM79 = 1/(1. + 1e4*sexp(920./gv.bin[nd]->tedust));
3027
3028 /* NB IntArea is total area per H for default abundances, not projected area, must div by 4
3029 * units s^-1 since gv.bin[nd]->cnv_H_pCM3 has units H cm-3
3030 * final units are cm s-1*/
3031 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.5 * AveVelH * gv.bin[nd]->IntArea/4. *
3032 /* cnv_H_pCM3 converts <unit>/H (default depletion) -> <unit>/cm^3 (actual depletion), units are cm-3 */
3033 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * conversion_efficiency_HM79;
3034 ASSERT( gv.bin[nd]->rate_h2_form_grains_HM79 > 0. );
3035 }
3036
3037# ifndef IGNORE_QUANTUM_HEATING
3038 if( lgUseQHeat )
3039 {
3040 /* H2 formation on grains from
3041 * >>refer H2 form Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */
3042 /* number of monolayers per second - only affects efficiency at very low or high temperatures */
3043 double f = 1e-10;
3044 /* equation 17
3045 double sqrt_term = POW2( 1. + sqrt( (10000.-200.)/(600.-200.) ) );*/
3046 double sqrt_term = 35.399494936611667;
3047
3048 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;
3049
3050 for( k=0; k < qnbin; k++ )
3051 {
3052 double beta_alpha = 0.25 * sqrt_term *sexp(200./qtemp[k] );
3053 /* equation 16 */
3054 double xi = 1./ (1. + 1.3e13*sexp(1.5*1e4/qtemp[k])*sqrt_term/(2.*f) );
3055 /* expression for beta comes from just after equation 5 */
3056 double beta = 3e12 * sexp( 320. / qtemp[k] );
3057 /* recombination efficiency given by their equation 15, they call
3058 * this epsilon_H2 */
3059 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );
3060 sticking_probability_H = 1./(1. + 0.04*sqrt(qtemp[k]+phycon.te) +
3061 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
3062
3063 /* printf( " k %ld Td %.6e re*sp %.6e\n", k, qtemp[k], recombination_efficiency_CT02* */
3064 /* sticking_probability_H ); */
3065
3066 gv.bin[nd]->rate_h2_form_grains_CT02 += qprob[k] * sticking_probability_H *
3067 recombination_efficiency_CT02;
3068 }
3069
3070 /* gv.bin[nd]->IntArea integrated grain surface area Int(4pi*a^2), normalized per H, in cm^2/H,
3071 * so x/4 is projected area of circle */
3072 /* gv.bin[nd]->cnv_H_pCM3 is H density [cm-3] times grain depletion factor */
3073 /* gv.bin[nd]->rate_h2_form_grains_CT02 units s-1 */
3074 gv.bin[nd]->rate_h2_form_grains_CT02 *= 0.5 * AveVelH *
3075 gv.bin[nd]->IntArea/4. * gv.bin[nd]->cnv_H_pCM3;
3076
3077 ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );
3078
3079 free(qtemp);
3080 free(qprob);
3081 }
3082 else
3083# endif
3084 {
3085 /* H2 formation on grains from
3086 * >>refer H2 form Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */
3087 /* number of monolayers per second - only affects efficiency at very low or high temperatures */
3088 double f = 1e-10;
3089 /* equation 17
3090 double sqrt_term = POW2( 1. + sqrt( (10000.-200.)/(600.-200.) ) );*/
3091 double sqrt_term = 35.399494936611667;
3092 double beta_alpha = 0.25 * sqrt_term *sexp(200./gv.bin[nd]->tedust );
3093 /* equation 16 */
3094 double xi = 1./ (1. + 1.3e13*sexp(1.5*1e4/ gv.bin[nd]->tedust )*sqrt_term/(2.*f) );
3095 /* expression for beta comes from just after equation 5 */
3096 double beta = 3e12 * sexp( 320. / gv.bin[nd]->tedust );
3097 /* recombination efficiency given by their equation 15, they call
3098 * this epsilon_H2 */
3099 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./SDIV(beta) + beta_alpha );
3100
3101 // hack for Leiden 2012 meeting, per description from Markus Roellig
3102 fixit(); // add a command-line option and refactor much of this and similar coding to remove spaghetti
3103#if 0
3104 /* carbonaceous grains */
3105 if( gv.bin[nd]->matType == MAT_CAR || gv.bin[nd]->matType == MAT_CAR2 )
3106 {
3107 double Td = gv.bin[nd]->tedust;
3108 recombination_efficiency_CT02 = 1./(
3109 ( 1. + 4.609569629185726e24*sexp(45000./Td) ) *
3110 ( 1. + sexp(800./Td) / (0.5389970511202651 * sexp(540./Td) + 5.6333909478365e-14*sqrt(Td)) )
3111 );
3112 //fprintf( ioQQQ, "DEBUGGG mole_reactions -- Setting recombination efficiency -- %3li\t%.3e\tcarbonaceous\n", nd, recombination_efficiency_CT02 );
3113 }
3114 /* silicate grains */
3115 else if( gv.bin[nd]->matType == MAT_SIL || gv.bin[nd]->matType == MAT_SIL2 )
3116 {
3117 double Td = gv.bin[nd]->tedust;
3118 recombination_efficiency_CT02 = 1./(
3119 ( 1. + 6.998161265697586e24*sexp(45000./Td) ) *
3120 ( 1. + sexp(450./Td) / (0.4266153643741715 * sexp(340./Td) + 2.5335919594255e-14*sqrt(Td)) )
3121 );
3122 //fprintf( ioQQQ, "DEBUGGG mole_reactions -- Setting recombination efficiency -- %3li\t%.3e\tsilicate\n", nd, recombination_efficiency_CT02 );
3123 }
3124#endif
3125
3126 /* gv.bin[nd]->IntArea integrated grain surface area Int(4pi*a^2), normalized per H, in cm^2/H,
3127 * so x/4 is projected area of circle */
3128 /* gv.bin[nd]->cnv_H_pCM3 is H density [cm-3] times grain depletion factor */
3129 /* units s-1 */
3130 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.5 * AveVelH * gv.bin[nd]->IntArea/4. *
3131 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_CT02;
3132 ASSERT( gv.bin[nd]->rate_h2_form_grains_CT02 > 0. );
3133 }
3134
3135# ifndef IGNORE_QUANTUM_HEATING
3136 /* reset sticking probability for code below */
3137 sticking_probability_H = 1./(1. + 0.04*sqrt(gv.bin[nd]->tedust+phycon.te) +
3138 0.002*phycon.te + 8e-6*phycon.te*phycon.te);
3139# endif
3140
3141 /* rate (s-1) all H2 v,J levels go to 0 or 1, preserving nuclear spin */
3142 /* ortho to para on grain surfaces, taken from
3143 *>refer H2 sticking Le Bourlot, J., 2000, A&A, 360, 656-662
3144 * >chng 05 apr 30, GS, hmi.H2_total/dense.gas_phase[ipHYDROGEN] is removed
3145 * This is used in h2.c.
3146 * NB IntArea is total are per H, not projected area, must div by 4
3147 * gv.bin[nd]->cnv_H_pCM3 has units H cm-3 to product with above
3148 * is cm2/H H/cm3 or cm-1 or an opacity
3149 * multiply by velocity of H2, cm s-1, so product
3150 * h2.rate_grain_op_conserve has units s^-1 */
3151 h2.rate_grain_op_conserve += AveVelH2 * gv.bin[nd]->IntArea/4. *
3152 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H;
3153
3154 /* ortho to para on grain surfaces, taken from
3155 *>refer H2 sticking Le Bourlot, J., 2000, A&A, 360, 656-662
3156 * For all grain temperatures, this process corresponds to high J going to
3157 * either 0 or 1 preserving nuclear spin. All ortho go to 1 and para go to 0.
3158 * When the dust temperature is below Tcrit all 1 go to 0 and so all J go to 0.
3159
3160 * this temperature depends on grain composition, discussion left column of page 657,
3161 * this is for a bare grain */
3169
3170 /* AveVelH2 is average speed of H2 molecules
3171 * for now assume that sticking probability for H2 on the grain is equal to
3172 * that for H
3173 * efficiency factor efficiency_opr is vary fast function of t dust -
3174 * large at low Td and small at Td > T_ortho_para_crit
3175 * start evaluating just above the critical temperature
3176 * T_ortho_para_crit this is roughly 24.345 K,GS */
3177 T_ortho_para_crit = 2. * hmi.Tad / log( POW2(60. *1.1e11)*hmi.Tad);
3178 if( gv.bin[nd]->tedust < T_ortho_para_crit )
3179 {
3180 double efficiency_opr = sexp(60.*1.1e11*sqrt(hmi.Tad)*sexp(hmi.Tad/gv.bin[nd]->tedust));
3181 /* rate (s-1) all v,J levels go to 0, regardless of nuclear spin
3182 * see above discussion for how units work out */
3183 h2.rate_grain_J1_to_J0 += AveVelH2 * gv.bin[nd]->IntArea/4. *
3184 gv.bin[nd]->cnv_H_pCM3 * sticking_probability_H * efficiency_opr;
3185 }
3186 }
3187 /*fprintf(ioQQQ," H2 grain form rate HM79 %.2e %.2e CT02 %.2e %.2e O-P grn %.2e %.2e\n",
3188 gv.bin[nd]->rate_h2_form_grains_HM79 ,
3189 gv.bin[nd]->rate_h2_form_grains_HM79 ,
3190 gv.bin[nd]->rate_h2_form_grains_CT02 ,
3191 gv.bin[nd]->rate_h2_form_grains_CT02 ,
3192 h2.rate_grain_J1_to_J0,
3193 hmi.rate_h2_allX_2_J1_grains
3194 );*/
3195 /* options to turn off grain collision with atom h2 collisions grains off command */
3196 h2.rate_grain_op_conserve *= h2.lgH2_grain_deexcitation;
3197 h2.rate_grain_J1_to_J0 *= h2.lgH2_grain_deexcitation;
3198
3199 }
3200 else
3201 {
3202 /* grains are not enabled, set these to zero */
3203 for( size_t nd=0; nd < gv.bin.size(); nd++ )
3204 {
3205 gv.bin[nd]->rate_h2_form_grains_CT02 = 0.;
3206 gv.bin[nd]->rate_h2_form_grains_HM79 = 0.;
3207 }
3208 /* rate all H2 goes to either 0 or 1 depending on ortho/para */
3209 h2.rate_grain_op_conserve = 0.;
3210 /* at low temp, rate all H2 goes to J=0 */
3211 h2.rate_grain_J1_to_J0 = 0.;
3212 }
3213
3214 /* the H2 catalysis rate on grains that is actually used in calculations
3215 * hmi.ScaleJura is scale factor set with set Jura scale command
3216 * units are s-1
3217 * default is 'C' Cazaux & Tielens */
3218 gv.rate_h2_form_grains_used_total = 0.;
3219 for( size_t nd=0; nd < gv.bin.size(); nd++ )
3220 {
3221 if( hmi.chJura == 'C' )
3222 {
3223 /* use the new rate by
3224 * >>refer H2 form Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29
3225 * units are s-1*/
3226 gv.bin[nd]->rate_h2_form_grains_used =
3227 gv.bin[nd]->rate_h2_form_grains_CT02*hmi.ScaleJura;
3228 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
3229 }
3230 else if( hmi.chJura == 'T' )
3231 {
3232 /* rate from Hollenbach & McKee 1979 */
3233 gv.bin[nd]->rate_h2_form_grains_used =
3234 gv.bin[nd]->rate_h2_form_grains_HM79*hmi.ScaleJura;
3235 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
3236 }
3237 else if( hmi.chJura == 'S' )
3238 {
3239 /* H2 formation rate from
3240 * >>refer H2 form Sternberg, A. & Neufeld, D.A. 1999, ApJ, 516, 371 */
3241 gv.bin[nd]->rate_h2_form_grains_used =
3242 3e-18 * phycon.sqrte / gv.bin.size() * dense.gas_phase[ipHYDROGEN]*hmi.ScaleJura;
3243 /* this is simple rate from Sternberg & Neufeld 99 */
3244 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
3245 }
3246 /*>>chng 07 jan 10, this had been C for constant, and so could never have been triggered.
3247 * caught by robin Williams, fixed by nick Abel, error was in sense that any set jura rate
3248 * would use Cazaux & Tielens */
3249 else if( hmi.chJura == 'F' )
3250 {
3251 /* command "set H2 rate" - enters log of Jura rate - C for constant,
3252 * no dependence on grain properties */
3253 gv.bin[nd]->rate_h2_form_grains_used = hmi.rate_h2_form_grains_set*dense.gas_phase[ipHYDROGEN] / gv.bin.size();
3254 gv.rate_h2_form_grains_used_total += gv.bin[nd]->rate_h2_form_grains_used;
3255 }
3256 }
3257 ASSERT( gv.rate_h2_form_grains_used_total >= 0. );
3258
3259# ifndef IGNORE_QUANTUM_HEATING
3260 printf( " fnzone %.2f H2 rate %.4e\n", fnzone, gv.rate_h2_form_grains_used_total );
3261# endif
3262
3263 /* print rate coefficient */
3264 /*fprintf(ioQQQ," total grain h2 form rate %.3e\n",gv.rate_h2_form_grains_used_total);*/
3265
3266}
3267/*mole_h_reactions update mole reactions for H */
3269{
3270 static double teused=-1;
3271 double exph2,
3272 exph2s,
3273 exphp,
3274 ex3hp;
3275 long i;
3276 double h2esc,
3277 th2,
3278 cr_H2s ,
3279 cr_H2dis,
3280 cr_H2dis_ELWERT_H2g,
3281 cr_H2dis_ELWERT_H2s;
3282
3283 DEBUG_ENTRY( "hmole_reactions()" );
3284
3285 /* everything here depends only on temperature - don't do anything if we don't
3286 * need to */
3287
3288 bool need_update = ! fp_equal( phycon.te, teused );
3289
3290 if( need_update )
3291 {
3292 teused = phycon.te;
3293
3294 /* get LTE populations */
3295 /* related to population of H- in LTE
3296 * IP is 0.754 eV */
3297 hmi.exphmi = sexp(8.745e3/phycon.te);
3298 if( hmi.exphmi > 0. )
3299 {
3300 /* these are ratio n*(H-)/[ n*(ne) n*(Ho) ] */
3301 hmi.rel_pop_LTE_Hmin = SAHA/(phycon.te32*hmi.exphmi)*(1./(2.*2.));
3302 }
3303 else
3304 {
3305 hmi.rel_pop_LTE_Hmin = 0.;
3306 }
3307
3308 /* population of H2+ in LTE, hmi.rel_pop_LTE_H2p is H_2+/H / H+
3309 * dissociation energy is 2.647 */
3310 exphp = sexp(3.072e4/phycon.te);
3311 if( exphp > 0. )
3312 {
3313 /* stat weight of H2+ is 4
3314 * last factor was put in ver 85.23, missing before */
3315 hmi.rel_pop_LTE_H2p = SAHA/(phycon.te32*exphp)*(4./(2.*1.))*3.634e-5;
3316 }
3317 else
3318 {
3319 hmi.rel_pop_LTE_H2p = 0.;
3320 }
3321
3322 /* related to population of H3+ in LTE, hmi.rel_pop_LTE_H3p is H_3+/( H2+ H+ )
3323 * dissociation energy is 2.647 */
3324 ex3hp = sexp(1.882e4/phycon.te);
3325 if( ex3hp > 0. )
3326 {
3327 /* stat weight of H2+ is 4
3328 * last factor was put in ver 85.23, missing before */
3329 hmi.rel_pop_LTE_H3p = SAHA/(phycon.te32*ex3hp)*(4./(2.*1.))*3.634e-5;
3330 }
3331 else
3332 {
3333 hmi.rel_pop_LTE_H3p = 0.;
3334 }
3335 }
3336 /* end constant temperature - */
3337
3338 // Big H2 rates are dependent on population as well as temperature
3339 /* population of H2 in LTE
3340 * dissociation energy of H2g is 4.477eV, for TH85 model */
3341 /* >>chng 05 oct 17, GS, averaged in big H2 in order to consider correct statistical weight*/
3342 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
3343 {
3344 /* the terms on the right were computed in the large molecule */
3345 hmi.rel_pop_LTE_H2g = h2.rel_pop_LTE_g;
3346 hmi.rel_pop_LTE_H2s = h2.rel_pop_LTE_s;
3347 }
3348 else
3349 {
3350 if (need_update)
3351 {
3352 /* H2 ground */
3353 exph2 = sexp((5.195e4)/phycon.te);
3354 /* >>chng 05 oct 17, GS, note that statical weight of H2g is assumed to be 1 if big H2 is not turned on*/
3355
3356 if( exph2 > 0. )
3357 {
3358 /* >>chng 05 oct 17, GS, note that statical weight of H2g is assumed to be 1 if big H2 is not turned on*/
3359 hmi.rel_pop_LTE_H2g = SAHA/(phycon.te32*exph2)*(1./(2.*2.))*3.634e-5;
3360 }
3361 else
3362 {
3363 hmi.rel_pop_LTE_H2g = 0.;
3364 }
3365
3366 /* H2 star */
3367 /* population of H2s in LTE
3368 * dissociation energy is 1.877eV, if h2s = 2.6eV, assumed for TH85 model */
3369 /* >>chng 05 oct 17, GS, averaged in big H2 in order to consider correct statistical weight*/
3370 exph2s = sexp(2.178e4/phycon.te);
3371
3372 if( exph2s > 0. )
3373 {
3374 /* >>chng 05 oct 17, GS, note that statical weight of H2s is assumed to be 1 if big H2 is not turned on*/
3375 hmi.rel_pop_LTE_H2s = SAHA/(phycon.te32*exph2s)*(1./(2.*2.))*3.634e-5;
3376 }
3377 else
3378 {
3379 hmi.rel_pop_LTE_H2s = 0.;
3380 }
3381 }
3382 }
3383 {
3384 /*@-redef@*/
3385 /* often the H- route is the most efficient formation mechanism for H2,
3386 * will be through rate called Hmolec_old[ipMH]*hmi.assoc_detach
3387 * this debug print statement is to trace h2 oscillations */
3388 enum {DEBUG_LOC=false};
3389 /*@+redef@*/
3390 if( DEBUG_LOC && nzone>187&& iteration > 1/**/)
3391 {
3392 /* rapid increase in H2 density caused by rapid increase in hmi.rel_pop_LTE_H2g */
3393 fprintf(ioQQQ,"ph2lteee\t%.2e\t%.1e\t%.1e\n",
3394 hmi.rel_pop_LTE_H2g,
3395 sexp(2.178e4/phycon.te),
3396 phycon.te);
3397 }
3398 }
3399
3400 /* cooling due to radiative attachment */
3401 hmi.hmicol = hmirat(phycon.te)*EN1RYD*phycon.te*1.15e-5;
3402
3403 fixit(); // Wasted cycles if we don't use Stancil's rates below? Why not put this down there/if used?
3404 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
3405 {
3406 if( (*diatom)->lgEnabled && mole_global.lgStancil )
3407 (*diatom)->Mol_Photo_Diss_Rates();
3408 }
3409
3410 /*fprintf(ioQQQ,"%.2e %.2e %.2e %.2e\n", phycon.te, hmi.hminus_rad_attach , hmi.hmicol,
3411 hmi.hmicol/(hmi.hminus_rad_attach*EN1RYD*phycon.te*1.15e-5) );*/
3412
3413 /* get per unit vol */
3414 hmi.hmicol *= dense.eden*findspecieslocal("H")->den; /* was dense.xIonDense[ipHYDROGEN][0]; */
3415
3416 /* ================================================================= */
3417 /* evaluate H- photodissociation rate, induced rec and rec cooling rates */
3418 /* >>chng 00 dec 24, add test so that photo rate only reevaluated two times per zone.
3419 * in grain-free models this was sometimes dominated by Lya and so oscillated.
3420 * especially bad in primal.in - change 2 to 4 and primal.in will stop due to Lya oscil */
3423 /* >>chng 02 feb 16, add damper on H- photo rate, wild oscillations in Lya photo rate in
3424 * grain free models */
3425
3426 t_phoHeat photoHeat;
3427
3428 /*hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , nhe1Com.nhe1[0] , opac.iphmop ,
3429 0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling );*/
3430 hmi.HMinus_photo_rate = GammaBn( hmi.iphmin-1 , iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon , opac.iphmop ,
3431 0.055502 , &hmi.HMinus_induc_rec_rate , &hmi.HMinus_induc_rec_cooling, &photoHeat );
3432
3433 /* save H- photodissociation heating */
3434 hmi.HMinus_photo_heat = photoHeat.HeatNet;
3435
3436 {
3437 /* following should be set true to print populations */
3438 /*@-redef@*/
3439 enum {DEBUG_LOC=false};
3440 /*@+redef@*/
3441 if( DEBUG_LOC)
3442 {
3443 fprintf(ioQQQ,"hminphoto\t%li\t%li\t%.2e\n", nzone, conv.nPres2Ioniz , hmi.HMinus_photo_rate );
3444 }
3445 }
3446
3447 /* induced recombination */
3448 hmi.HMinus_induc_rec_rate *= hmi.rel_pop_LTE_Hmin*dense.eden;
3449
3450 /* induced recombination cooling per unit volume */
3452 hmi.HMinus_induc_rec_cooling *= hmi.rel_pop_LTE_Hmin*dense.eden*findspecieslocal("H")->den; /* dense.xIonDense[ipHYDROGEN][0]; */
3453
3454 {
3455 /* following should be set true to debug H- photoionization rates */
3456 /*@-redef@*/
3457 enum {DEBUG_LOC=false};
3458 /*@+redef@*/
3459 if( DEBUG_LOC && nzone>400/*&& iteration > 1*/)
3460 {
3461 fprintf(ioQQQ,"hmoledebugg %.2f ",fnzone);
3462 GammaPrt(
3463 hmi.iphmin-1 , iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon , opac.iphmop ,
3464 /* io unit we will write to */
3465 ioQQQ,
3466 /* total photo rate from previous call to GammaK */
3467 hmi.HMinus_photo_rate,
3468 /* we will print contributors that are more than this rate */
3469 hmi.HMinus_photo_rate*0.05);
3470 }
3471 }
3472 /* add on high energy ionization, assume hydrogen cross section
3473 * n.b.; HGAMNC with secondaries */
3474 /* >>chng 00 dec 24, above goes to HeI edge, no need for this, and was not important*/
3475 /*hmi.HMinus_photo_rate += iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].gamnc;*/
3476
3477 /* ================================================================= */
3478 /* photodissociation by Lyman band absorption: esc prob treatment,
3479 * treatment based on
3480 * >>refer HI abs Tielens & Hollenbach 1985 ApJ 291, 722. */
3481 /* do up to carbon photo edge if carbon is turned on */
3482 /* >>>chng 00 apr 07, add test for whether element is turned on */
3483 hmi.UV_Cont_rel2_Habing_TH85_depth = 0.;
3484 /* >>chng 00 apr 07 from explicit ipHeavy to ipLo */
3485 /* find total intensity over carbon-ionizing continuum */
3486 /* >>chng 03 jun 09, use exact bounds rather than CI photo threshold for lower bound */
3487 /*for( i=ipLo; i < iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].ipIsoLevNIonCon; i++ )*/
3488 /* the integral is from 6eV to 13.6, or 2060 - 912 Ang */
3489 for( i=rfield.ipG0_TH85_lo; i < rfield.ipG0_TH85_hi; ++i )
3490 {
3491 hmi.UV_Cont_rel2_Habing_TH85_depth += ((rfield.flux[0][i-1]) + (rfield.ConInterOut[i-1])+
3492 (rfield.outlin[0][i-1])+ (rfield.outlin_noplot[i-1]))*rfield.anu[i-1];
3493 }
3494
3495 /* now convert to Habing ISM units
3496 * UV_Cont_rel2_Habing_TH85_face is FUV continuum relative to Habing value
3497 * 1.6e-3 ergs em-2 s-1 is the Habing 1968 field, quoted on page 723, end of first
3498 * full paragraph on left */
3499 hmi.UV_Cont_rel2_Habing_TH85_depth = (realnum)(hmi.UV_Cont_rel2_Habing_TH85_depth*EN1RYD/1.6e-3);
3500 /* if start of calculation remember G0 at illuminated face */
3501 if( nzone<=1 )
3502 {
3503 hmi.UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_depth;
3504 }
3505
3506
3507 /* >>chng 05 jan 09, add special ver of G0 */
3508 hmi.UV_Cont_rel2_Habing_spec_depth = 0.;
3509 for( i=rfield.ipG0_spec_lo; i < rfield.ipG0_spec_hi; ++i )
3510 {
3511 hmi.UV_Cont_rel2_Habing_spec_depth += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
3512 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1])*rfield.anu[i-1];
3513 }
3514 hmi.UV_Cont_rel2_Habing_spec_depth = (realnum)(hmi.UV_Cont_rel2_Habing_spec_depth*EN1RYD/1.6e-3);
3515
3516 /* the Draine & Bertoldi version of the same thing, defined over their energy range */
3517 /* >>chng 04 feb 07, only evaluate at the illuminated face */
3518 if( !conv.nTotalIoniz )
3519 {
3520 hmi.UV_Cont_rel2_Draine_DB96_face = 0.f;
3521 /* this is sum of photon number between 912A and 1110, as per BD96 */
3522 for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; ++i )
3523 {
3524 hmi.UV_Cont_rel2_Draine_DB96_face += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
3525 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1]);
3526 }
3527 /* Convert into scaled ISM background field, total number of photons over value for 1 ISM field,
3528 * the coefficient 1.232e7 is the number of photons over this wavelength range for 1H and is
3529 * given in BD96, page 225, 4th line from start of section 4, also page 272, table 1, 2nd line
3530 * from bottom */
3531 /* >>chng 04 mar 16, introduce the 1.71 */
3532 /* equation 20 of DB96 gives chi as flux over 1.21e7, to produce one Habing field.
3533 * to get the Draine field we need to further divide by 1.71 as stated on the first
3534 * line after equation 23 */
3535 hmi.UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face/(1.232e7f * 1.71f);
3536 }
3537
3538 /* escape prob takes into account line shielding,
3539 * next is opacity then optical depth in H2 UV lines, using eqn A7 of TH85 */
3540 hmi.H2Opacity = (realnum)(1.2e-14*(1e5/GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])));
3541 /* the typical Lyman -Werner H2 line optical depth eq A7 of TH85a */
3542 th2 = (colden.colden[ipCOL_H2g]+ colden.colden[ipCOL_H2s])*hmi.H2Opacity;
3543 /* the escape probability - chance that continuum photon will penetrate to
3544 * this depth to pump the Lyman Werner bands */
3545 h2esc = esc_PRD_1side(th2,1e-4);
3546
3547 /* cross section is eqn A8 of
3548 * >>refer H2 dissoc Tielens, A.G.G.M., & Hollenbach, D., 1985, ApJ, 291, 722
3549 * branching ratio of 10% is already included, so 10x smaller than their number
3550 * 10% lead to dissociation through H_2 + h nu => 2H */
3551 /* >>chng 05 mar 10, by TE, break into 2g and 2s
3552 * note use of same shielding column in below - can do far better */
3553 hmi.H2_Solomon_dissoc_rate_TH85_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
3554 hmi.H2_Solomon_dissoc_rate_TH85_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
3555 hmi.H2_H2g_to_H2s_rate_TH85 = hmi.H2_Solomon_dissoc_rate_TH85_H2g*9.;
3556
3557 /* these are Burton et al. 1990 rates */
3558 hmi.H2_Solomon_dissoc_rate_BHT90_H2g = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
3559 hmi.H2_Solomon_dissoc_rate_BHT90_H2s = 3.4e-11 * hmi.UV_Cont_rel2_Habing_TH85_depth * h2esc;
3560 hmi.H2_H2g_to_H2s_rate_BHT90 = hmi.H2_Solomon_dissoc_rate_BHT90_H2g*9.;
3561
3562 {
3563 /* the following implements Drain & Bertoldi 1996, equation 37 from
3564 * >>refer H2 destruction Draine, B.T., & Bertoldi, F., 1996, ApJ, 468, 269-289
3565 * but the constant 4.6e-11 comes from Bertoldi & Draine equation 23,
3566 * this is the normalized H2 column density */
3567 double x = (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) / 5e14;
3568 double sqrtx = sqrt(1.+x);
3569 /* Doppler with of H2 in km/s */
3570 double b5 = GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])/1e5;
3571 /* the molecular hydrogen line self-shielding factor */
3572 double fshield = 0.965/POW2(1.+x/b5) + 0.035/sqrtx *
3573 sexp(8.5e-4*sqrtx);
3574
3575 /*double fshield = pow( MAX2(1.,colden.colden[ipCOLH2]/1e14) , -0.75 );*/
3576 /* this is the Bertoldi & Draine version of the Habing field,
3577 * with dilution of radiation and extinction due to grains */
3578 /* >>chng 04 apr 18, moved fshield, the line shielding factor, from this defn to
3579 * the following defn of dissociation rate, since following should measure continuum */
3580 hmi.UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_face *
3581 (realnum)(sexp( opac.TauAbsFace[rfield.ip1000A-1] )/radius.r1r0sq);
3582
3583 /* the following comes from Bertoldi & Draine 1996, equation 23,
3584 * hmi.UV_Cont_rel2_Draine_DB96_depth already includes a factor of 1.71 to correct back from TH85 */
3585 /* >>chng 05 mar 10, TE, break into 2s and 2s */
3586 if( !mole_global.lgLeidenHack )
3587 {
3588 /* this is default, when set Leiden hack UMIST rates not entered */
3589 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
3590 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 4.6e-11 * hmi.UV_Cont_rel2_Draine_DB96_depth * fshield;
3591 }
3592 else
3593 {
3594 /* done when set Leiden hack UMIST rates command entered */
3595 hmi.H2_Solomon_dissoc_rate_BD96_H2g = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
3596 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
3597 hmi.H2_Solomon_dissoc_rate_BD96_H2s = 5.18e-11* (hmi.UV_Cont_rel2_Habing_TH85_face/1.66f)
3598 *sexp(3.02*rfield.extin_mag_V_point)* fshield;
3599 }
3600
3601 /* BD do not give an excitation rate, so used 9 times the dissociation
3602 * rate by analogy with 90% going back into X, as per TH85 */
3603 /*>>chng 04 feb 07, had used 90% relax into X from TH85,
3604 * BD say that 15% dissociate, so 85/15 = 5.67 is factor */
3605 hmi.H2_H2g_to_H2s_rate_BD96 = 5.67* hmi.H2_Solomon_dissoc_rate_BD96_H2g;
3606 }
3607
3608 /* do Elwert approximation to the dissociation rate */
3609 if( hmi.UV_Cont_rel2_Draine_DB96_face > SMALLFLOAT )
3610 {
3611 /* this evaluates the new H2 dissociation rate by Torsten Elwert */
3612 /* chng 05 jun 23, TE
3613 * >>chng 05 sep 13, update master source with now approximation */
3614
3615 /* Doppler with of H2 in km/s */
3616 double b5 = GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN])/1e5;
3617
3618 /* split the Solomon rates in H2g and H2s */
3619 /* >>chng 05 oct 19,
3620 * >>chng 05 dec 05, TE, define new approximation for the heating due to the destruction of H2
3621 * use this approximation for the specified cloud parameters, otherwise
3622 * use limiting cases for 1 <= G0, G0 >= 1e7, n >= 1e7, n <= 1 */
3623
3624 double x_H2g, x_H2s,
3625 fshield_H2g, fshield_H2s,
3626 f_H2s;
3627 static double a_H2g, a_H2s,
3628 e1_H2g, e1_H2s,
3629 e2_H2g,
3630 b_H2g,
3631 sl_H2g, sl_H2s,
3632 k_f_H2s,
3633 k_H2g_to_H2s,
3634 log_G0_face = -1;
3635
3636 /* define parameter range for the new approximation
3637 * test for G0
3638 *BUGFIX - this tested on lg_G0_face < 0 for initialization needed so did not work
3639 * in grids - change to evaluate in zone 0 */
3640 /* >>chng 07 feb 24, BUGFIX crash when G0=0 at start and radiation
3641 * field builds up due to diffuse fields - soln is to always reevaluate */
3642 /*if( !nzone )*/
3643 {
3644 if(hmi.UV_Cont_rel2_Draine_DB96_face <= 1.)
3645 {
3646 log_G0_face = 0.;
3647 }
3648 else if(hmi.UV_Cont_rel2_Draine_DB96_face >= 1e7)
3649 {
3650 log_G0_face = 7.;
3651 }
3652 else
3653 {
3654 log_G0_face = log10(hmi.UV_Cont_rel2_Draine_DB96_face);
3655 }
3656
3657 /* terms only dependent on G0_face */
3658
3659 /* coefficients and exponents */
3660 a_H2g = 0.06 * log_G0_face + 1.32;
3661 a_H2s = 0.375 * log_G0_face + 2.125;
3662
3663 e1_H2g = -0.05 * log_G0_face + 2.25;
3664 e1_H2s = -0.125 * log_G0_face + 2.625;
3665
3666 e2_H2g = -0.005 * log_G0_face + 0.625;
3667
3668 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
3669
3670 /* scalelength for H2g and H2s */
3671 sl_H2g = 4.0e14;
3672 sl_H2s = 9.0e15;
3673
3674 /* coefficient for 2nd term of Solomon H2s */
3675 k_f_H2s = MAX2(0.1,2.375 * log_G0_face - 1.875 );
3676
3677 /* coefficient for branching ratio */
3678 k_H2g_to_H2s = MAX2(1.,-1.75 * log_G0_face + 11.25);
3679
3680 /*fprintf( ioQQQ, "e1_H2g%.2e, e1_H2s%.2e, e2_H2g%.2e, b_H2g%.2e, a_H2g%.2e, a_H2s%.2e,sl_H2g: %.2e,sl_H2s: %.2e\n",
3681 e1_H2g, e1_H2s, e2_H2g, b_H2g, a_H2g, a_H2s, sl_H2g, sl_H2s);
3682 */
3683 }
3684
3685 /* Solomon H2s ~G0^0.2 at large depth*/
3686 f_H2s = k_f_H2s * pow((double)hmi.UV_Cont_rel2_Draine_DB96_depth, 0.2 );
3687
3688 /* scale length for absorption of UV lines */
3689 x_H2g = (colden.colden[ipCOL_H2g]) / sl_H2g;
3690 x_H2s = (colden.colden[ipCOL_H2s]) / sl_H2s;
3691
3692 /* the molecular hydrogen line self-shielding factor */
3693 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
3694 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
3695
3696 /* the Elwert Solomon rates for H2g and H2s hmi.chH2_small_model_type == 'E' */
3697 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = a_H2g * 4.6e-11 * fshield_H2g * hmi.UV_Cont_rel2_Draine_DB96_depth;
3698 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = a_H2s * 4.6e-11 * fshield_H2s * (hmi.UV_Cont_rel2_Draine_DB96_depth + f_H2s);
3699
3700 /* assume branching ratio dependent on G0*/
3701 hmi.H2_H2g_to_H2s_rate_ELWERT = k_H2g_to_H2s * hmi.H2_Solomon_dissoc_rate_ELWERT_H2g;
3702
3703 /* use G0_BD96 as this definition declines faster with depth which is physical as
3704 * the longer wavelengths in the definition of G0_TH85 cannot dissociate
3705 * H2s directly */
3706 hmi.H2_photodissoc_ELWERT_H2s = hmi.UV_Cont_rel2_Draine_DB96_depth*1e-11;
3707 hmi.H2_photodissoc_ELWERT_H2g = hmi.H2_photodissoc_ELWERT_H2s * 1.0e-10;
3708 }
3709 else
3710 {
3711 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g = 0.;
3712 hmi.H2_Solomon_dissoc_rate_ELWERT_H2s = 0.;
3713 hmi.H2_photodissoc_ELWERT_H2s = 0.;
3714 hmi.H2_photodissoc_ELWERT_H2g = 0.;
3715 }
3716
3717 /* this is rate of photodissociation of H2*, A12 of TH85 */
3718 hmi.H2_photodissoc_TH85 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
3719
3720 /* dissociation rate from Burton et al. 1990 */
3721 hmi.H2_photodissoc_BHT90 = hmi.UV_Cont_rel2_Habing_TH85_depth*1e-11;
3722
3723 /* rates for cosmic ray excitation of singlet bound electronic bound excited states
3724 * only add this to small molecule since automatically included in large
3725 *>>refer H2 cr excit Dalgarno, A., Yan, Min, & Liu, Weihong 1999, ApJS, 125, 237
3726 * this is excitation of H2* */
3727 /* >>chng 05 sep 13, do not include this process when Leiden hacks are in place */
3728 cr_H2s = secondaries.x12tot*0.9 / 2. * hmi.lgLeidenCRHack;
3729 /* this is the fraction that dissociate */
3730 /* >>chng 05 sep 13, do not include this process when Leiden hacks are in place */
3731 cr_H2dis = secondaries.x12tot*0.1 / 2. * hmi.lgLeidenCRHack;
3732
3733 /* >>chng 05 sep 13, TE, synchronize treatment of CR */
3734 /* cosmic ray rates for dissociation of ground and H2s
3735 * two factors done to agree with large H2 deep in the cloud where
3736 * cosmic rays are important */
3737 cr_H2dis_ELWERT_H2g = secondaries.x12tot*5e-8 * hmi.lgLeidenCRHack;
3738 cr_H2dis_ELWERT_H2s = secondaries.x12tot*4e-2 * hmi.lgLeidenCRHack;
3739
3740 /* at this point there are two or three independent estimates of the H2 dissociation rate.
3741 * if the large H2 molecule is on, then H2 Solomon rates has been defined in the last
3742 * call to the large molecule. Just above we have defined hmi.H2_Solomon_dissoc_rate_TH85,
3743 * the dissociation rate from Tielens & Hollenbach 1985, and hmi.H2_Solomon_dissoc_rate_BD96,
3744 * the rate from Bertoldi & Draine 1996. We can use any defined rate. If the big H2
3745 * molecule is on, use its rate. If not, for now use the TH85 rate, since that is the
3746 * rate we always have used in the past.
3747 * The actual rate we will use is given by hmi.H2_Solomon_dissoc_rate_used
3748 */
3749 /* this is the Solomon process dissociation rate actually used */
3750 if( h2.lgEnabled && h2.lgEvaluated && hmi.lgH2_Chemistry_BigH2 )
3751 {
3752 /* only update after big H2 molecule has been evaluated,
3753 * when very little H2 and big molecule not used, leave at previous (probably TH85) value,
3754 * since that value is always known */
3755
3756 /* Solomon process rate from X into the X continuum with units s-1
3757 * rates are total rate, and rates from H2g and H2s */
3758 hmi.H2_Solomon_dissoc_rate_used_H2g = h2.Solomon_dissoc_rate_g;
3759 ASSERT( hmi.H2_Solomon_dissoc_rate_used_H2g >= 0. );
3760
3761 hmi.H2_Solomon_dissoc_rate_used_H2s = h2.Solomon_dissoc_rate_s;
3762 ASSERT( hmi.H2_Solomon_dissoc_rate_used_H2s >= 0. );
3763
3764 /* photoexcitation from H2g to H2s */
3765 hmi.H2_H2g_to_H2s_rate_used = h2.gs_rate();
3766 ASSERT( hmi.H2_H2g_to_H2s_rate_used >= 0. );
3767
3768 /* add up H2s + hnu (continuum) => 2H + KE, continuum photodissociation,
3769 * this is not the Solomon process, true continuum, units s-1 */
3770 /* only include rates from H2s since this is only open channel, this process is well
3771 * shielded against Lyman continuum destruction by atomic hydrogen */
3772 hmi.H2_photodissoc_used_H2s = h2.photodissoc_BigH2_H2s;
3773 /* NPA - 07/24/09 - logic to use Stancil photodissociation rate for H2s */
3774 if( mole_global.lgStancil && h2.lgEnabled )
3775 hmi.H2_photodissoc_used_H2s = h2.Cont_Dissoc_Rate_H2s;
3776 ASSERT( hmi.H2_photodissoc_used_H2s >= 0. );
3777
3778 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
3779 * for unfavorable wavelength range of G0*/
3780 hmi.H2_photodissoc_used_H2g = h2.photodissoc_BigH2_H2g;
3781 /* NPA - 07/24/09 - logic to use Stancil photodissociation rate for H2g */
3782 if( mole_global.lgStancil && h2.lgEnabled )
3783 hmi.H2_photodissoc_used_H2g = h2.Cont_Dissoc_Rate_H2g;
3784 ASSERT( hmi.H2_photodissoc_used_H2g >= 0. );
3785 }
3786 else if( hmi.chH2_small_model_type == 'T' )
3787 {
3788 /* the TH85 rate */
3789 /*>>chng 05 jun 23, add cosmic rays */
3790 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_TH85_H2g + cr_H2dis;
3791 /* >>chng 05 sep 13, cr_H2dis was not included */
3792 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_TH85_H2s + cr_H2dis;
3793 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_TH85 + cr_H2s;
3794
3795 /* continuum photodissociation H2s + hnu => 2H, ,
3796 * this is not the Solomon process, true continuum, units s-1 */
3797 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
3798 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
3799 * for unfavorable wavelength range of G0*/
3800 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
3801 }
3802
3803 else if( hmi.chH2_small_model_type == 'H' )
3804 {
3805 /* the improved H2 formalism given by
3806 *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M
3807 >>refcon 1990, ApJ, 365, 620 */
3808 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BHT90_H2g + cr_H2dis;
3809 /* >>chng 05 sep 13, cr_H2dis was not included */
3810 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BHT90_H2s + cr_H2dis;
3811 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BHT90 + cr_H2s;
3812
3813 /* continuum photodissociation H2s + hnu => 2H, ,
3814 * this is not the Solomon process, true continuum, units s-1 */
3815 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_BHT90;
3816 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
3817 * for unfavorable wavelength range of G0*/
3818 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_BHT90*1.0e-10f;
3819 }
3820
3821 else if( hmi.chH2_small_model_type == 'B' )
3822 {
3823 /* the Bertoldi & Draine rate - this is the default */
3824 /*>>chng 05 jun 23, add cosmic rays */
3825 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_BD96_H2g + cr_H2dis;
3826 /* >>chng 05 sep 13, cr_H2dis was not included */
3827 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_BD96_H2s + cr_H2dis;
3828 /* they did not do the excitation or dissoc rate, so use TH85 */
3829 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_BD96 + cr_H2s;
3830
3831
3832 /* continuum photodissociation H2s + hnu => 2H, ,
3833 * this is not the Solomon process, true continuum, units s-1 */
3834 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_TH85;
3835 /* >>chng 05 mar 24, TE, continuum photodissociation rate of H2g, small correction factor accounts
3836 * for unfavorable wavelength range of G0*/
3837 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_TH85*1.0e-10f;
3838 }
3839 else if( hmi.chH2_small_model_type == 'E' )
3840 {
3841 /* the Elwert et al. rate
3842 *>>chng 05 jun 23, add cosmic rays */
3843 hmi.H2_Solomon_dissoc_rate_used_H2g = hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + cr_H2dis_ELWERT_H2g;
3844 hmi.H2_Solomon_dissoc_rate_used_H2s = hmi.H2_Solomon_dissoc_rate_ELWERT_H2s + cr_H2dis_ELWERT_H2s;
3845 hmi.H2_H2g_to_H2s_rate_used = hmi.H2_H2g_to_H2s_rate_ELWERT + cr_H2s;
3846
3847
3848 /* continuum photodissociation H2s + hnu => 2H, ,
3849 * this is not the Solomon process, true continuum, units s-1 */
3850 hmi.H2_photodissoc_used_H2s = hmi.H2_photodissoc_ELWERT_H2s;
3851 hmi.H2_photodissoc_used_H2g = hmi.H2_photodissoc_ELWERT_H2g;
3852 }
3853 else
3854 TotalInsanity();
3855
3856 {
3857 /*@-redef@*/
3858 enum {DEBUG_LOC=false};
3859 /*@+redef@*/
3860 if( DEBUG_LOC && h2.lgEnabled )
3861 {
3862 fprintf(ioQQQ," Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
3863 hmi.H2_Solomon_dissoc_rate_TH85_H2g,
3864 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
3865 h2.Solomon_dissoc_rate_g,
3866 hmi.H2_H2g_to_H2s_rate_TH85 ,
3867 h2.gs_rate() );
3868 }
3869 }
3870 return;
3871}
3873{
3874 DEBUG_ENTRY("mole_findrate_s()");
3875
3876 string newbuf = canonicalize_reaction_label(buf);
3877
3878 mole_reaction_i p = mole_priv::reactab.find(newbuf);
3879
3880 if(p != mole_priv::reactab.end())
3881 return &(*p->second);
3882 else
3883 return NULL;
3884}
3885
3886double t_mole_local::findrk(const char buf[]) const
3887{
3888 DEBUG_ENTRY("t_mole_local::findrk()");
3889
3890 mole_reaction *rate = mole_findrate_s(buf);
3891
3892 if(!rate)
3893 return 0.;
3894
3895 /* check for NaN */
3896 ASSERT( !isnan( reaction_rks[ rate->index ] ) );
3897
3898 return reaction_rks[ rate->index ];
3899}
3900double t_mole_local::findrate(const char buf[]) const
3901{
3902 double ret;
3903 int i;
3904
3905 DEBUG_ENTRY("t_mole_local::findrate()");
3906
3907 mole_reaction *rate = mole_findrate_s(buf);
3908 if(!rate)
3909 {
3910 return 0.;
3911 }
3912
3913 ret = reaction_rks[ rate->index ];
3914 for(i=0;i<rate->nreactants;i++)
3915 ret *= species[ rate->reactants[i]->index ].den;
3916
3917 return ret;
3918}
3919/* Calculate rate at which molecular network abstracts species */
3920
3921/* Need to check reactants vs reactant behaviour */
3922double t_mole_local::sink_rate_tot(const char chSpecies[]) const
3923{
3924 DEBUG_ENTRY("t_mole_local::sink_rate_tot()");
3925
3926 const molecule* const sp = findspecies(chSpecies);
3927 double ratev = sink_rate_tot(sp);
3928
3929 return ratev;
3930}
3931double t_mole_local::sink_rate_tot(const molecule* const sp) const
3932{
3933 DEBUG_ENTRY("t_mole_local::sink_rate_tot()");
3934 double ratev = 0;
3935
3936 for(mole_reaction_i p=mole_priv::reactab.begin();
3937 p != mole_priv::reactab.end(); ++p)
3938 {
3939 mole_reaction &rate = *p->second;
3940 ratev += sink_rate( sp, rate );
3941 }
3942
3943 return ratev;
3944}
3945
3946double t_mole_local::sink_rate(const molecule* const sp, const char buf[]) const
3947{
3948 const mole_reaction* const rate = mole_findrate_s(buf);
3949 return sink_rate( sp, *rate );
3950}
3951
3952double t_mole_local::sink_rate(const molecule* const sp, const mole_reaction& rate) const
3953{
3954 DEBUG_ENTRY("t_mole_local::sink_rate()");
3955
3956 int ipthis = -1;
3957 for(int i=0;i<rate.nreactants && ipthis == -1;i++)
3958 {
3959 if(rate.reactants[i] == sp && rate.rvector[i]==NULL && rate.rvector_excit[i]==NULL )
3960 {
3961 ipthis = i;
3962 }
3963 }
3964 if(ipthis != -1)
3965 {
3966 double ratevi = rate.a * rate.rk();
3967 for(int i=0;i<rate.nreactants;i++)
3968 {
3969 if(i!=ipthis)
3970 {
3971 ratevi *= species[ rate.reactants[i]->index ].den;
3972 }
3973 }
3974 return ratevi;
3975 }
3976 else
3977 return 0.;
3978}
3979
3980#ifdef PRINTREACT
3981STATIC void printreact(mole_reaction *rate)
3982{
3983 int i;
3984
3985 DEBUG_ENTRY("printreact()");
3986
3987 for(i=0;i<rate->nreactants;i++) {
3988 fprintf(stderr,"%s,",rate->reactants[i]->label);
3989 }
3990 fprintf(stderr,"=>");
3991 for(i=0;i<rate->nproducts;i++) {
3992 fprintf(stderr,"%s,",rate->products[i]->label);
3993 }
3994 fprintf(stderr,"\n");
3995
3996}
3997#endif
3999double t_mole_local::dissoc_rate(const char chSpecies[]) const
4000{
4001 DEBUG_ENTRY("t_mole_local::dissoc_rate()");
4002
4004 if (sp == null_mole)
4005 return 0.0;
4006 ASSERT(sp->isMonatomic());
4007 const chem_atom *tgt = sp->nAtom.begin()->first.get_ptr();
4008 molecule *ph = findspecies("PHOTON");
4009 double ratev = 0.0;
4010
4011 for (mole_reaction_i p
4012 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
4013 {
4014 mole_reaction &rate = *p->second;
4015
4016 // Must have a photon in to be a dissociation rate
4017 int ipph = 0;
4018 for (int i=0;i<rate.nreactants;i++)
4019 {
4020 if (rate.reactants[i] == ph)
4021 ipph++;
4022 }
4023 if (!ipph)
4024 continue;
4025
4026 // ipsp is number of *specific* species of interest,
4027 // ipfree is number in same ionization ladder, including X-
4028 int ipspin = 0, ipfreein = 0;
4029 for (int i=0;i<rate.nreactants;i++)
4030 {
4031 if (rate.reactants[i] == sp)
4032 ++ipspin;
4033 if (rate.reactants[i]->isMonatomic() && tgt == sp->nAtom.begin()->first.get_ptr())
4034 ++ipfreein;
4035 }
4036 int ipspout = 0, ipfreeout = 0;
4037 for (int i=0;i<rate.nproducts;i++)
4038 {
4039 if (rate.products[i] == sp)
4040 ++ipspout;
4041 if (rate.products[i]->isMonatomic() && tgt == sp->nAtom.begin()->first.get_ptr())
4042 ++ipfreeout;
4043 }
4044
4045 // Must produce the species requested
4046 int newsp = ipspout-ipspin;
4047 if (newsp <= 0)
4048 continue;
4049
4050 // And must do so by breaking bonds
4051 int nbondsbroken = ipfreeout-ipfreein;
4052 if (nbondsbroken <= 0)
4053 continue;
4054 // Fraction of the generated monatomic species which were *originally* bound
4055 double fracbroken = nbondsbroken/((double)ipfreeout);
4056 ASSERT( fracbroken <= 1.0 );
4057
4058 double ratevi = reaction_rks[ rate.index ];
4059 for (int i=0;i<rate.nreactants;i++)
4060 {
4061 ratevi *= species[ rate.reactants[i]->index ].den;
4062 }
4063
4064 // Photoproduction rate is rate of production of the species
4065 // which has not come from an initially monatomic source
4066
4067 double ratesp = ratevi*newsp; // This is the total production
4068 // rate of the specific species
4069 ratesp *= fracbroken; // Scale back for any initially unbound
4070 // monatoms
4071
4072 ratev += ratesp;
4073 }
4074 return ratev;
4075}
4077{
4078 DEBUG_ENTRY("t_mole_local::source_rate_tot()");
4079
4081 double ratev = source_rate_tot(sp);
4082
4083 return ratev;
4084}
4085double t_mole_local::source_rate_tot(const molecule* const sp) const
4086{
4087 DEBUG_ENTRY("t_mole_local::source_rate_tot()");
4088 double ratev = 0;
4089
4090 for (mole_reaction_i p =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
4091 {
4092 mole_reaction &rate = *p->second;
4093 int ipthis = 0;
4094 for(int i=0;i<rate.nproducts;i++)
4095 {
4096 if( rate.products[i] == sp && rate.pvector[i]==NULL && rate.pvector_excit[i]==NULL )
4097 {
4098 ipthis++;
4099 }
4100 }
4101 if(ipthis)
4102 {
4103 double ratevi = rate.a * rate.rk();
4104 for(int i=0;i<rate.nreactants;i++)
4105 {
4106 ratevi *= species[ rate.reactants[i]->index ].den;
4107 }
4108 ratev += ipthis*ratevi;
4109 }
4110 }
4111
4112 return ratev;
4113}
4114
4115double t_mole_local::chem_heat(void) const
4116{
4117 /* >>chng 07, Feb 11 NPA. Calculate the chemical heating rate. This is defined as the net energy of the
4118 * reaction, which is:
4119 *
4120 * Energy = SUM[formation energies of reactants] - SUM[formation energies of products]
4121 *
4122 * Now take the energy, and multiply by the densities of the reactants and the rate constant, finally
4123 * you have to multiply by 1.66e-14, which is the conversion factor to go from kJ/mol to erg/atom
4124 * this gives the units in the form of erg/atom*cm3/s*cm-3*cm-3 = erg/cm-3/s/atom, which is
4125 * a heating rate
4126 */
4127
4128 DEBUG_ENTRY("t_mole_local::chem_heat()");
4129
4130 double heating = 0.;
4131 map<double,string> heatMap;
4132 molecule *ph = findspecies("PHOTON");
4133 molecule *crph = findspecies("CRPHOT");
4134 molecule *grn = findspecies("grn");
4135
4136 /* loop over all reactions */
4137 for (mole_reaction_i p
4138 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
4139 {
4140 mole_reaction &rate = *p->second;
4141
4142 // If PHOTON appears, assume it accounts for energy difference
4143 bool lgCanSkip = false;
4144 for (int i=0;i<rate.nproducts;i++)
4145 {
4146 if( rate.products[i] == ph || rate.products[i] == crph )
4147 lgCanSkip = true;
4148 }
4149 for (int i=0;i<rate.nreactants;i++)
4150 {
4151 if( rate.reactants[i] == ph || rate.reactants[i] == crph )
4152 lgCanSkip = true;
4153 }
4154 // grain catalyst reactions are handled in grain physics, don't double count here
4155 for (int i=0;i<rate.nreactants;i++)
4156 {
4157 if( rate.reactants[i] == grn && rate.rvector[i] != NULL )
4158 lgCanSkip = true;
4159 }
4160
4161 if( lgCanSkip )
4162 continue;
4163
4164 /* This loop calculates the product of the rate constant times the densities*/
4165 double rate_tot = reaction_rks[ rate.index ];
4166 for( long i=0; i < rate.nreactants; ++i )
4167 {
4168 rate_tot *= species[ rate.reactants[i]->index ].den;
4169 }
4170
4171 realnum reaction_enthalpy = 0.;
4172
4173 /* Calculate the sum of the formation energies for the reactants */
4174 for( long i=0; i < rate.nreactants; ++i )
4175 {
4176 reaction_enthalpy += rate.reactants[i]->form_enthalpy;
4177 }
4178
4179 /* Subtract from that the sum of the formation energies of the products */
4180 for( long i=0; i < rate.nproducts; ++i )
4181 {
4182 reaction_enthalpy -= rate.products[i]->form_enthalpy;
4183 }
4184
4185 /* this is the chemical heating rate. TODO. Once the H chem is merged with the C chem, then
4186 * we will have the chemical heating rate for all reactions. This is only a subset and, thusfar,
4187 * not actually used in getting the total heating. Tests with pdr_leiden_hack_f1.in show that this
4188 * heating rate can be up to 10% of the total heating */
4189
4190 double heat = reaction_enthalpy*rate_tot*(1e10/AVOGADRO); /* 1.66e-14f; */
4191 heatMap[heat] = rate.label;
4192 heating += heat;
4193 }
4194
4195 // use reverse iterator to print out biggest contributors
4196 long index = 0;
4197 // this should be a const_reverse_iterator, but pgCC 12.2-0 64-bit cannot handle this.
4198 // it appears there is no const version of heatMap.rend(); as a result the compiler
4199 // cannot find a suitable version of operator != in "it != heatMap.rend()"
4200 // The Solaris Studio compiler version 12.3 has the same problem
4201 for( map<double,string>::reverse_iterator it = heatMap.rbegin(); it != heatMap.rend(); ++it, ++index )
4202 {
4203 fprintf( ioQQQ, "DEBUGGG heat %li\t%li\t%.6e\t%s\n", index, nzone, it->first, it->second.c_str() );
4204 if( index==2 )
4205 break;
4206 }
4207 index = 0;
4208 for( map<double,string>::iterator it = heatMap.begin(); it != heatMap.end(); ++it, ++index )
4209 {
4210 if( it->first >= 0. )
4211 break;
4212 fprintf( ioQQQ, "DEBUGGG cool %li\t%li\t%.6e\t%s\n", index, nzone, it->first, it->second.c_str() );
4213 if( index==2 )
4214 break;
4215 }
4216
4217 return heating;
4218}
4219
4220/* Punch all rates involving specified species */
4221void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth)
4222{
4223 int i, ipthis;
4224 double ratevi;
4225 molecule *sp;
4226
4227 DEBUG_ENTRY("mole_punch()");
4228
4229 sp = findspecies(speciesname);
4230
4231 if (lgHeader)
4232 {
4233 fprintf (punit,"#Depth");
4234 }
4235 if (lgData)
4236 {
4237 fprintf (punit,"%.5e",depth);
4238 }
4239
4240 for (mole_reaction_i p
4241 =mole_priv::reactab.begin(); p != mole_priv::reactab.end(); ++p)
4242 {
4243 mole_reaction &rate = *p->second;
4244 ipthis = 0;
4245
4246 for (i=0;i<rate.nreactants;i++)
4247 {
4248 if( rate.reactants[i] == sp )
4249 {
4250 if( ( strcmp( args, "DEST" )==0 && rate.rvector[i]==NULL && rate.rvector_excit[i]==NULL ) ||
4251 ( strcmp( args, "DFLT" )==0 && rate.rvector[i]==NULL && rate.rvector_excit[i]==NULL ) ||
4252 ( strcmp( args, "CATA" )==0 && rate.rvector[i]!=NULL ) ||
4253 strcmp( args, "ALL " )==0 )
4254 ipthis++;
4255 }
4256 }
4257
4258 for(i=0;i<rate.nproducts;i++)
4259 {
4260 if( rate.products[i] == sp )
4261 {
4262 if( ( strcmp( args, "CREA" )==0 && rate.pvector[i]==NULL && rate.pvector_excit[i]==NULL ) ||
4263 ( strcmp( args, "DFLT" )==0 && rate.pvector[i]==NULL && rate.pvector_excit[i]==NULL ) ||
4264 ( strcmp( args, "CATA" )==0 && rate.pvector[i]!=NULL ) ||
4265 strcmp( args, "ALL " )==0 )
4266 ipthis++;
4267 }
4268 }
4269
4270 if(ipthis)
4271 {
4272 if (lgHeader)
4273 {
4274 fprintf(punit,"\t%s",rate.label.c_str());
4275 }
4276 if (lgData)
4277 {
4278 ratevi = mole.reaction_rks[ rate.index ];
4279 for(i=0;i<rate.nreactants;i++)
4280 {
4281 ratevi *= mole.species[ rate.reactants[i]->index ].den;
4282 }
4283 fprintf(punit,"\t%.3e",ratevi);
4284 }
4285 }
4286 }
4287 fprintf(punit,"\n");
4288}
4289
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 isnan
Definition cddefines.h:620
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
const int ipCARBON
Definition cddefines.h:310
@ CHARS_SPECIES
Definition cddefines.h:274
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
const int ipOXYGEN
Definition cddefines.h:312
#define POW2
Definition cddefines.h:929
T sign(T x, T y)
Definition cddefines.h:800
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
double dsexp(double x)
Definition service.cpp:953
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
double powi(double, long int)
Definition service.cpp:604
#define HMRATE(a, b, c)
Definition cddefines.h:1046
vector< GrainBin * > bin
Definition grainvar.h:583
static t_version & Inst()
Definition cddefines.h:175
bool lgEnabled
Definition h2_priv.h:345
double spon_diss_tot
Definition h2_priv.h:262
double Solomon_dissoc_rate_g
Definition h2_priv.h:264
bool lgEvaluated
Definition h2_priv.h:310
double photodissoc_BigH2_H2g
Definition h2_priv.h:258
double photoionize_rate
Definition h2_priv.h:254
double Average_A
Definition h2_priv.h:296
molecule * pvector_excit[MAXPRODUCTS]
Definition mole_priv.h:58
virtual double rk() const =0
molecule * pvector[MAXPRODUCTS]
Definition mole_priv.h:57
molecule * rvector_excit[MAXREACTANTS]
Definition mole_priv.h:55
molecule * reactants[MAXREACTANTS]
Definition mole_priv.h:53
molecule * rvector[MAXREACTANTS]
Definition mole_priv.h:54
string label
Definition mole_priv.h:51
molecule * products[MAXPRODUCTS]
Definition mole_priv.h:56
double reduced_mass
Definition mole_priv.h:59
realnum form_enthalpy
Definition mole.h:164
map< const count_ptr< chem_atom >, int, element_pointer_value_less > nAtomsMap
Definition mole.h:135
int charge
Definition mole.h:144
bool isEnabled
Definition mole.h:139
nAtomsMap nAtom
Definition mole.h:143
realnum mole_mass
Definition mole.h:165
string label
Definition mole.h:142
bool isMonatomic(void) const
Definition mole.h:157
int index
Definition mole.h:169
double den
Definition mole.h:358
double **** PhotoRate_Shell
Definition ionbal.h:111
vector< freeBound > fb
Definition iso.h:452
qList st
Definition iso.h:453
bool lgLeidenHack
Definition mole.h:286
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
double dissoc_rate(const char chSpecies[]) const
vector< double > reaction_rks
Definition mole.h:400
double sink_rate(const molecule *const sp, const mole_reaction &rate) const
double source_rate_tot(const char chSpecies[]) const
double findrate(const char buf[]) const
double findrk(const char buf[]) const
double chem_heat(void) const
valarray< class molezone > species
Definition mole.h:398
double sink_rate_tot(const char chSpecies[]) const
t_colden colden
Definition colden.cpp:5
#define ipCOL_H2s
Definition colden.h:18
#define ipCOL_H2g
Definition colden.h:16
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
double GammaBn(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double thresh, double *ainduc, double *rcool, t_phoHeat *photoHeat)
t_conv conv
Definition conv.cpp:5
bool operator<(const count_ptr< T > &a, const count_ptr< T > &b)
Definition count_ptr.h:88
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_deuterium deut
Definition deuterium.cpp:8
t_DoppVel DoppVel
Definition doppvel.cpp:5
realnum GetDopplerWidth(realnum massAMU)
realnum GetAveVelocity(realnum massAMU)
void GrainDrive(void)
Definition grains.cpp:1591
void qheat(vector< double > &, vector< double > &, long *, size_t)
GrainVar gv
Definition grainvar.cpp:5
@ MAT_SIL
Definition grainvar.h:102
@ MAT_SIL2
Definition grainvar.h:105
@ MAT_CAR2
Definition grainvar.h:104
@ MAT_CAR
Definition grainvar.h:101
const int NQGRID
Definition grainvar.h:32
diatomics hd("hd", 4100., &hmi.HD_total, Yan_H2_CS)
vector< diatomics * > diatoms
Definition h2.cpp:8
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_hextra hextra
Definition hextra.cpp:5
t_hmi hmi
Definition hmi.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH3p
Definition iso.h:31
const int ipH1s
Definition iso.h:27
const int ipHE_LIKE
Definition iso.h:63
const int ipH3s
Definition iso.h:30
const int ipH2p
Definition iso.h:29
void iso_photo(long ipISO, long nelem)
const int ipH3d
Definition iso.h:32
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
void create_isotopologues_one(ChemAtomList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, vector< string > &newLabels)
molezone * findspecieslocal(const char buf[])
mole_reaction * mole_findrate_s(const char buf[])
ChemAtomList unresolved_atom_list
molecule::nAtomsMap::iterator nAtoms_i
Definition mole.h:214
bool lgDifferByExcitation(const molecule &mol1, const molecule &mol2)
bool parse_species_label(const char label[], ChemAtomList &atomsLeftToRight, vector< int > &numAtoms, string &embellishments)
molecule * findspecies(const char buf[])
molecule * null_mole
vector< count_ptr< chem_atom > > ChemAtomList
Definition mole.h:115
@ CONFLICT
Definition mole_priv.h:69
@ ABSENT
Definition mole_priv.h:69
@ CORRECT
Definition mole_priv.h:69
#define MAXREACTANTS
Definition mole_priv.h:45
#define MAXPRODUCTS
Definition mole_priv.h:46
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
Definition mole_priv.h:37
map< string, count_ptr< mole_reaction > >::const_iterator mole_reaction_ci
Definition mole_priv.h:38
void mole_create_react(void)
@ UDFA
STATIC double mole_get_equilibrium_constant(const char buf[])
STATIC void newreact(const char label[], const char fun[], double a, double b, double c)
STATIC void compare_udfa(const count_ptr< mole_reaction > &rate)
STATIC bool lgReactBalance(const count_ptr< mole_reaction > &rate)
STATIC void read_data(const char file[], void(*parse)(char *s))
STATIC string canonicalize_reaction_label(const char label[])
STATIC long parse_reaction(count_ptr< mole_reaction > &rate, const char label[])
STATIC void mole_check_reverse_reactions(void)
void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth)
STATIC char * getcsvfield(char **s, char c)
void mole_rk_bigchange(void)
STATIC void mole_h_reactions(void)
STATIC void mole_generate_isotopologue_reactions(string atom_old, string atom_new)
STATIC void plot_sparsity(void)
STATIC void parse_base(char *s)
#define FLTEQ(a, b)
mole_reaction * mole_findrate_s(const char buf[])
STATIC double mole_partition_function(const molecule *const sp)
STATIC void mole_h2_grain_form(void)
double hmirat(double te)
double frac_H2star_hminus(void)
STATIC void canonicalize_reaction(count_ptr< mole_reaction > &rate)
static map< formula_species, struct udfa_reaction_s * > udfatab
STATIC void register_reaction_vectors(count_ptr< mole_reaction > rate)
@ BUFSIZE
static realnum albedo
void mole_update_rks(void)
void mole_cmp_num_in_out_reactions()
STATIC void parse_udfa(char *s)
map< string, count_ptr< mole_reaction > > functab
map< string, count_ptr< mole_reaction > > reactab
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double PI
Definition physconst.h:29
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double SAHA
Definition physconst.h:161
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double EN1EV
Definition physconst.h:192
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
UNUSED const double HBAReV
Definition physconst.h:204
UNUSED const double AVOGADRO
Definition physconst.h:106
UNUSED const double HION_LTE_POP
Definition physconst.h:157
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
double esc_PRD_1side(double tau, double a)
double esca0k2(double taume)
t_secondaries secondaries
static double * source
Definition species2.cpp:28
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition atmdat.h:153
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
Definition atmdat.h:152
long int nTotalIoniz
Definition conv.h:166
double eden
Definition dense.h:190
double xIonDense[LIMELM][LIMELM+1]
Definition dense.h:125
realnum AtomicWeight[LIMELM]
Definition dense.h:75
realnum cryden_ov_background
Definition hextra.h:25
bool lgH2_Chemistry_BigH2
Definition hmi.h:164
realnum UV_Cont_rel2_Draine_DB96_depth
Definition hmi.h:74
long int iheh2
Definition hmi.h:59
double rel_pop_LTE_Hmin
Definition hmi.h:194
double H2_Solomon_dissoc_rate_used_H2g
Definition hmi.h:92
double H2_H2g_to_H2s_rate_used
Definition hmi.h:89
double HMinus_induc_rec_rate
Definition hmi.h:55
double H2_Solomon_dissoc_rate_used_H2s
Definition hmi.h:98
double H2_photodissoc_used_H2g
Definition hmi.h:108
double HMinus_photo_rate
Definition hmi.h:42
double H2_photodissoc_used_H2s
Definition hmi.h:109
double rel_pop_LTE_H2s
Definition hmi.h:197
long int iheh1
Definition hmi.h:58
long int ih2pof
Definition opacity.h:230
long int ih2pnt[2]
Definition opacity.h:229
double HeatNet
Definition thermal.h:172
double te05
Definition phycon.h:57
double te
Definition phycon.h:11
double te90
Definition phycon.h:50
double sqrte
Definition phycon.h:48
double te30
Definition phycon.h:53
realnum * outlin_noplot
Definition rfield.h:200
realnum ** flux
Definition rfield.h:86
realnum ** outlin
Definition rfield.h:199
long int nflux
Definition rfield.h:43
realnum * ConInterOut
Definition rfield.h:164
realnum ** csupra
Definition secondaries.h:21
realnum x12tot
Definition secondaries.h:53
char ** chSpecies
Definition taulines.cpp:13
t_trace trace
Definition trace.cpp:5