48 STATIC void newreact(
const char label[],
const char fun[],
double a,
double b,
double c);
88 template<
class T>
void newfunc()
94 ratefunc findfunc(
const char name[]);
97 class mole_reaction_hmrate;
99 class mole_reaction_th85rate;
101 class mole_reaction_crnurate_noalbedo;
105 class mole_reaction_grn_react;
107 class mole_reaction_h2_collh_excit;
109 class mole_reaction_h2_collh2_excit;
111 class mole_reaction_h2_collh_deexc;
113 class mole_reaction_h2_collh2_deexc;
115 class mole_reaction_rh2g_dis_h;
117 class mole_reaction_rh2s_dis_h;
119 class mole_reaction_rh2g_dis_h2;
121 class mole_reaction_rh2s_dis_h2;
123 class mole_reaction_rh2s_dis_h2_nodeexcit;
125 class mole_reaction_bh2g_dis_h;
127 class mole_reaction_bh2s_dis_h;
129 class mole_reaction_bh2g_dis_h2;
131 class mole_reaction_hneut;
133 class mole_reaction_cionhm;
150 typedef mole_reaction_null T;
152 virtual T*
Create()
const {
return new T;}
153 virtual const char*
name() {
return "null";}
191 for(n=0;n<nreact;n++)
220 r *= pow(te/300.,rate->
b);
222 r *= exp(-rate->
c/te);
228 typedef mole_reaction_hmrate_exo T;
230 virtual T*
Create()
const {
return new T;}
232 virtual const char*
name() {
return "hmrate_exo";}
243 te =
MIN2( te, -10. * this->c );
245 return pow(te/300.,this->b)*exp(-te/this->c);
251 typedef mole_reaction_hmrate T;
253 virtual T*
Create()
const {
return new T;}
254 virtual const char*
name() {
return "hmrate";}
263 typedef mole_reaction_constrate T;
265 virtual T*
Create()
const {
return new T;}
266 virtual const char*
name() {
return "constrate";}
307 typedef mole_reaction_th85rate T;
309 virtual T*
Create()
const {
return new T;}
310 virtual const char*
name() {
return "th85rate";}
313 return th85rate(
this);
339 typedef mole_reaction_crnurate_noalbedo T;
341 virtual T*
Create()
const {
return new T;}
343 virtual const char*
name() {
return "crnurate_noalbedo";}
346 return crnurate_noalbedo(
this);
352 typedef mole_reaction_crnurate T;
354 virtual T*
Create()
const {
return new T;}
355 virtual const char*
name() {
return "crnurate";}
358 return crnurate_noalbedo(
this)/(1.0-
albedo);
367 typedef mole_reaction_cryden_ov_bg T;
369 virtual T*
Create()
const {
return new T;}
370 virtual const char*
name() {
return "cryden_ov_bg";}
379 typedef mole_reaction_co_lnu_c_o_lnu T;
381 virtual T*
Create()
const {
return new T;}
382 virtual const char*
name() {
return "co_lnu_c_o_lnu";}
394 for( ns=0; ns<2; ++ns )
440 typedef mole_reaction_vib_evap T;
442 virtual T*
Create()
const {
return new T;}
443 virtual const char*
name() {
return "vib_evap";}
446 double binding_energy, exponent, vib_freq;
452 binding_energy = this->
b;
453 double bin_total=0.0;
454 for(
size_t nd=0; nd <
gv.
bin.size() ; nd++ )
456 double bin_area =
gv.
bin[nd]->IntArea*
gv.
bin[nd]->cnv_H_pCM3;
457 exponent += exp(-binding_energy/
gv.
bin[nd]->tedust)*bin_area;
458 bin_total += bin_area;
460 exponent /= bin_total;
461 const double surface_density_of_sites = 1.5e15;
464 vib_freq = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*binding_energy/(PI*PI*this->reactants[0]->mole_mass));
484 typedef mole_reaction_grn_abs T;
486 virtual T*
Create()
const {
return new T;}
487 virtual const char*
name() {
return "grn_abs";}
495 if (this->reactants[0]->n_nuclei() != 0)
496 mass = this->reactants[0]->mole_mass;
498 mass = this->reactants[1]->mole_mass;
500 return sqrt(8.*BOLTZMANN*
phycon.
te/(PI*mass));
506 typedef mole_reaction_grn_react T;
508 virtual T*
Create()
const {
return new T;}
509 virtual const char*
name() {
return "grn_react";}
512 return grn_react(
this);
536 double quant_barrier = 1e-8;
537 double surface_density_of_sites = 1.5e15;
538 fixit(
"rate->a is _always_ overall scaling factor");
542 double activ_barrier = rate->
c;
545 fixit(
"Ordering of reactants may have switched");
546 double vib_freq_i = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_i/(PI*PI*rate->
reactants[0]->
mole_mass));
547 double vib_freq_j = sqrt(2*surface_density_of_sites*(0.3*BOLTZMANN)*E_j/(PI*PI*rate->
reactants[1]->
mole_mass));
551 double dust_density = 0.0;
554 fixit(
"should there be an area weighting to exponential terms?");
556 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
558 double bin_density =
gv.
bin[nd]->IntArea*
gv.
bin[nd]->cnv_H_pCM3;
559 Exp_i += exp(-E_i/
gv.
bin[nd]->tedust)*bin_density;
560 Exp_j += exp(-E_j/
gv.
bin[nd]->tedust)*bin_density;
561 dust_density += bin_density/(4*1e-10);
569 double diff_i = vib_freq_i*Exp_i/total_sites;
570 double diff_j = vib_freq_j*Exp_j/total_sites;
574 double scale = exp(-2*(quant_barrier/(HBAReV*EN1EV))*sqrt(2.0*rate->
reduced_mass*0.3*BOLTZMANN*activ_barrier));
577 return scale*(diff_i + diff_j)/
SDIV(dust_density);
582 typedef mole_reaction_grn_photo T;
584 virtual T*
Create()
const {
return new T;}
585 virtual const char*
name() {
return "grn_photo";}
600 fixit(
"grainn photoionization");
613 typedef mole_reaction_th85rate_co T;
615 virtual T*
Create()
const {
return new T;}
617 virtual const char*
name() {
return "th85rate_co";}
641 if (this->reactants[0]->n_nuclei() != 0)
642 column =
mole.
species[ this->reactants[0]->index ].column;
644 column =
mole.
species[ this->reactants[1]->index ].column;
646 esc_co = 4.4e-15 * column /
651 return esca0k2(esc_co)*th85rate(
this);
657 typedef mole_reaction_oh_c2h2_co_ch3 T;
659 virtual T*
Create()
const {
return new T;}
660 virtual const char*
name() {
return "oh_c2h2_co_ch3";}
679 typedef mole_reaction_h_hnc_hcn_h T;
681 virtual T*
Create()
const {
return new T;}
683 virtual const char*
name() {
return "h_hnc_hcn_h";}
719 return 1. - 4.938e-6;
724 double frac_H2star_grains(
void)
750 typedef mole_reaction_h2ph3p T;
752 virtual T*
Create()
const {
return new T;}
753 virtual const char*
name() {
return "h2ph3p";}
772 typedef mole_reaction_hopexch T;
774 virtual T*
Create()
const {
return new T;}
775 virtual const char*
name() {
return "hopexch";}
784 typedef mole_reaction_hpoexch T;
786 virtual T*
Create()
const {
return new T;}
787 virtual const char*
name() {
return "hpoexch";}
796 typedef mole_reaction_hmattach T;
798 virtual T*
Create()
const {
return new T;}
799 virtual const char*
name() {
return "hmattach";}
808 typedef mole_reaction_hmihph2p T;
810 virtual T*
Create()
const {
return new T;}
811 virtual const char*
name() {
return "hmihph2p";}
835 typedef mole_reaction_hmphoto T;
837 virtual T*
Create()
const {
return new T;}
838 virtual const char*
name() {
return "hmphoto";}
870 typedef mole_reaction_cionhm T;
872 virtual T*
Create()
const {
return new T;}
873 virtual const char*
name() {
return "cionhm";}
880 double assoc_detach(
void)
895 y=545969508.1323510+x*71239.23653059864;
901 typedef mole_reaction_c3bod T;
903 virtual T*
Create()
const {
return new T;}
904 virtual const char*
name() {
return "c3bod";}
913 typedef mole_reaction_asdfg T;
915 virtual T*
Create()
const {
return new T;}
916 virtual const char*
name() {
return "asdfg";}
925 typedef mole_reaction_asdfs T;
927 virtual T*
Create()
const {
return new T;}
928 virtual const char*
name() {
return "asdfs";}
937 typedef mole_reaction_asdbg T;
939 virtual T*
Create()
const {
return new T;}
940 virtual const char*
name() {
return "asdbg";}
950 typedef mole_reaction_asdbs T;
952 virtual T*
Create()
const {
return new T;}
953 virtual const char*
name() {
return "asdbs";}
963 typedef mole_reaction_bhneut T;
965 virtual T*
Create()
const {
return new T;}
966 virtual const char*
name() {
return "bhneut";}
977 return (hneut(
this)*ratio*
1000 return 3.4738192887404660e-008;
1005 typedef mole_reaction_hneut T;
1007 virtual T*
Create()
const {
return new T;}
1008 virtual const char*
name() {
return "hneut";}
1018 typedef mole_reaction_h2_spon_diss T;
1020 virtual T*
Create()
const {
return new T;}
1021 virtual const char*
name() {
return "h2_spon_diss";}
1030 typedef mole_reaction_h2_ind_rad_diss T;
1032 virtual T*
Create()
const {
return new T;}
1033 virtual const char*
name() {
return "h2_ind_rad_diss";}
1043 fixit(
"Remove factor of dense.gas_phase[ipHYDROGEN] factor from"
1044 "derivation of rate_h2_form_grains_used to avoid"
1054 virtual const char*
name() {
return "grnh2tot";}
1057 return grnh2tot(
this);
1063 typedef mole_reaction_grnh2 T;
1065 virtual T*
Create()
const {
return new T;}
1066 virtual const char*
name() {
return "grnh2";}
1069 return grnh2tot(
this)*(1.-frac_H2star_grains());
1075 typedef mole_reaction_grnh2s T;
1077 virtual T*
Create()
const {
return new T;}
1078 virtual const char*
name() {
return "grnh2s";}
1081 return grnh2tot(
this)*frac_H2star_grains();
1087 typedef mole_reaction_radasc T;
1089 virtual T*
Create()
const {
return new T;}
1090 virtual const char*
name() {
return "radasc";}
1098 return hmrate(
this) *
1110 typedef mole_reaction_assoc_ion T;
1112 virtual T*
Create()
const {
return new T;}
1113 virtual const char*
name() {
return "assoc_ion";}
1120 return hmrate(
this) *
1151 typedef mole_reaction_rh2g_dis_h T;
1153 virtual T*
Create()
const {
return new T;}
1154 virtual const char*
name() {
return "rh2g_dis_h";}
1157 return rh2g_dis_h(
this);
1172 fprintf(
ioQQQ,
"invalid parameter for rh2s_dis_h\n" );
1180 typedef mole_reaction_rh2s_dis_h T;
1182 virtual T*
Create()
const {
return new T;}
1183 virtual const char*
name() {
return "rh2s_dis_h";}
1186 return rh2s_dis_h(
this);
1202 fprintf(
ioQQQ,
"invalid parameter for rh2g_dis_h2\n" );
1205 return hmrate4(5.5e-29*0.5/(SAHA*3.634e-5)*sqrt(300.),0.5,5.195e4,
phycon.
te);
1210 typedef mole_reaction_rh2g_dis_h2 T;
1212 virtual T*
Create()
const {
return new T;}
1213 virtual const char*
name() {
return "rh2g_dis_h2";}
1216 return rh2g_dis_h2(
this);
1231 fprintf(
ioQQQ,
"invalid parameter for rh2s_dis_h2\n" );
1239 typedef mole_reaction_rh2s_dis_h2 T;
1241 virtual T*
Create()
const {
return new T;}
1242 virtual const char*
name() {
return "rh2s_dis_h2";}
1245 return rh2s_dis_h2(
this);
1260 fprintf(
ioQQQ,
"invalid parameter for rh2s_dis_h2_nodeexcit\n" );
1266 class mole_reaction_rh2s_dis_h2_nodeexcit :
public mole_reaction
1268 typedef mole_reaction_rh2s_dis_h2_nodeexcit T;
1270 virtual T*
Create()
const {
return new T;}
1271 virtual const char*
name() {
return "rh2s_dis_h2_nodeexcit";}
1274 return rh2s_dis_h2_nodeexcit(
this);
1284 typedef mole_reaction_bh2g_dis_h T;
1286 virtual T*
Create()
const {
return new T;}
1287 virtual const char*
name() {
return "bh2g_dis_h";}
1290 return bh2g_dis_h(
this);
1300 typedef mole_reaction_bh2s_dis_h T;
1302 virtual T*
Create()
const {
return new T;}
1303 virtual const char*
name() {
return "bh2s_dis_h";}
1306 return bh2s_dis_h(
this);
1316 typedef mole_reaction_bh2g_dis_h2 T;
1318 virtual T*
Create()
const {
return new T;}
1319 virtual const char*
name() {
return "bh2g_dis_h2";}
1322 return bh2g_dis_h2(
this);
1328 typedef mole_reaction_bh2s_dis_h2 T;
1330 virtual T*
Create()
const {
return new T;}
1331 virtual const char*
name() {
return "bh2s_dis_h2";}
1340 typedef mole_reaction_h2photon T;
1342 virtual T*
Create()
const {
return new T;}
1343 virtual const char*
name() {
return "h2photon";}
1354 typedef mole_reaction_h2crphot T;
1356 virtual T*
Create()
const {
return new T;}
1357 virtual const char*
name() {
return "h2crphot";}
1366 typedef mole_reaction_h2crphh T;
1368 virtual T*
Create()
const {
return new T;}
1369 virtual const char*
name() {
return "h2crphh";}
1394 typedef mole_reaction_h2scrphh T;
1396 virtual T*
Create()
const {
return new T;}
1397 virtual const char*
name() {
return "h2scrphh";}
1418 typedef mole_reaction_radath T;
1420 virtual T*
Create()
const {
return new T;}
1421 virtual const char*
name() {
return "radath";}
1430 typedef mole_reaction_gamtwo T;
1432 virtual T*
Create()
const {
return new T;}
1433 virtual const char*
name() {
return "gamtwo";}
1449 typedef mole_reaction_hlike_phot T;
1451 virtual T*
Create()
const {
return new T;}
1452 virtual const char*
name() {
return "hlike_phot";}
1468 typedef mole_reaction_h2s_sp_decay T;
1470 virtual T*
Create()
const {
return new T;}
1471 virtual const char*
name() {
return "h2s_sp_decay";}
1500 typedef mole_reaction_h2_collh2_deexc T;
1502 virtual T*
Create()
const {
return new T;}
1503 virtual const char*
name() {
return "h2_collh2_deexc";}
1506 return h2_collh2_deexc(
this);
1523 typedef mole_reaction_h2_collh_deexc T;
1525 virtual T*
Create()
const {
return new T;}
1526 virtual const char*
name() {
return "h2_collh_deexc";}
1529 return h2_collh_deexc(
this);
1547 typedef mole_reaction_h2_collh2_excit T;
1549 virtual T*
Create()
const {
return new T;}
1550 virtual const char*
name() {
return "h2_collh2_excit";}
1553 return h2_collh2_excit(
this);
1571 typedef mole_reaction_h2_collh_excit T;
1573 virtual T*
Create()
const {
return new T;}
1574 virtual const char*
name() {
return "h2_collh_excit";}
1577 return h2_collh_excit(
this);
1583 typedef mole_reaction_h2gexcit T;
1585 virtual T*
Create()
const {
return new T;}
1586 virtual const char*
name() {
return "h2gexcit";}
1595 typedef mole_reaction_h2sdissoc T;
1597 virtual T*
Create()
const {
return new T;}
1598 virtual const char*
name() {
return "h2sdissoc";};
1613 typedef mole_reaction_h2gdissoc T;
1615 virtual T*
Create()
const {
return new T;}
1616 virtual const char*
name() {
return "h2gdissoc";}
1631 typedef mole_reaction_hd_photodissoc T;
1633 virtual T*
Create()
const {
return new T;}
1634 virtual const char*
name() {
return "hd_photodissoc";}
1665 else if( te < 1200. )
1669 else if( te < 3800. )
1674 else if( te <= 1.4e4 )
1698 typedef mole_reaction_gamheh T;
1700 virtual T*
Create()
const {
return new T;}
1701 virtual const char*
name() {
return "gamheh";}
1713 for( i=
hmi.
iheh1-1; i < limit; i++ )
1726 enum {exclude, base, umisthack, federman, lithium, deuterium, ti, misc, in_code, generated};
1730 static bool lgReactInitialized =
false;
1740 if( lgReactInitialized )
1743 lgReactInitialized =
true;
1753 newfunc<mole_reaction_null>();
1755 map<string,bool> canonical;
1764 newfunc<mole_reaction_hmrate>();
1765 newfunc<mole_reaction_hmrate_exo>();
1766 newfunc<mole_reaction_constrate>();
1767 newfunc<mole_reaction_th85rate>();
1768 newfunc<mole_reaction_crnurate>();
1769 newfunc<mole_reaction_crnurate_noalbedo>();
1770 newfunc<mole_reaction_co_lnu_c_o_lnu>();
1771 newfunc<mole_reaction_vib_evap>();
1772 newfunc<mole_reaction_th85rate_co>();
1773 newfunc<mole_reaction_grn_abs>();
1775 newfunc<mole_reaction_grn_react>();
1777 newfunc<mole_reaction_grn_photo>();
1778 newfunc<mole_reaction_oh_c2h2_co_ch3>();
1779 newfunc<mole_reaction_h_hnc_hcn_h>();
1781 newfunc<mole_reaction_gamheh>();
1782 newfunc<mole_reaction_hd_photodissoc>();
1783 newfunc<mole_reaction_h2gdissoc>();
1784 newfunc<mole_reaction_h2sdissoc>();
1785 newfunc<mole_reaction_h2gexcit>();
1786 newfunc<mole_reaction_h2_collh_excit>();
1787 newfunc<mole_reaction_h2_collh2_excit>();
1788 newfunc<mole_reaction_h2_collh_deexc>();
1789 newfunc<mole_reaction_h2_collh2_deexc>();
1790 newfunc<mole_reaction_h2s_sp_decay>();
1791 newfunc<mole_reaction_hlike_phot>();
1792 newfunc<mole_reaction_gamtwo>();
1793 newfunc<mole_reaction_h2ph3p>();
1794 newfunc<mole_reaction_radath>();
1795 newfunc<mole_reaction_cryden_ov_bg>();
1796 newfunc<mole_reaction_h2scrphh>();
1797 newfunc<mole_reaction_h2crphh>();
1798 newfunc<mole_reaction_h2photon>();
1799 newfunc<mole_reaction_h2crphot>();
1800 newfunc<mole_reaction_rh2g_dis_h>();
1801 newfunc<mole_reaction_rh2s_dis_h>();
1802 newfunc<mole_reaction_rh2g_dis_h2>();
1803 newfunc<mole_reaction_rh2s_dis_h2>();
1804 newfunc<mole_reaction_rh2s_dis_h2_nodeexcit>();
1805 newfunc<mole_reaction_bh2g_dis_h>();
1806 newfunc<mole_reaction_bh2s_dis_h>();
1807 newfunc<mole_reaction_bh2g_dis_h2>();
1808 newfunc<mole_reaction_bh2s_dis_h2>();
1809 newfunc<mole_reaction_radasc>();
1810 newfunc<mole_reaction_assoc_ion>();
1811 newfunc<mole_reaction_grnh2>();
1812 newfunc<mole_reaction_grnh2s>();
1813 newfunc<mole_reaction_h2_spon_diss>();
1814 newfunc<mole_reaction_h2_ind_rad_diss>();
1815 newfunc<mole_reaction_bhneut>();
1816 newfunc<mole_reaction_hneut>();
1817 newfunc<mole_reaction_asdbs>();
1818 newfunc<mole_reaction_asdbg>();
1819 newfunc<mole_reaction_asdfs>();
1820 newfunc<mole_reaction_asdfg>();
1821 newfunc<mole_reaction_c3bod>();
1822 newfunc<mole_reaction_cionhm>();
1823 newfunc<mole_reaction_hmphoto>();
1824 newfunc<mole_reaction_hmihph2p>();
1825 newfunc<mole_reaction_hmattach>();
1860 newreact(
"C+,OH=>CO,H+",
"hmrate",0.,0.,0.);
1865 newreact(
"H,H,grn=>H2,grn",
"grnh2",1.,0.,0.);
1866 newreact(
"H,H,grn=>H2*,grn",
"grnh2s",1.,0.,0.);
1867 newreact(
"H-,PHOTON=>H,e-",
"hmphoto",1.,0.,0.);
1872 fixit(
"this should be atom_list instead of unresolved_element_list, but we have not defined ionized species of all isotopes yet!!!");
1876 if( !(*atom)->lgHasLinkedIon() )
1878 long nelem = (*atom)->el->Z-1;
1882 sprintf(react,
"H-,%s+=>H,%s", (*atom)->label().c_str(), (*atom)->label().c_str() );
1883 newreact(react,
"hmrate",4e-6/sqrt(300.),-0.5,0.);
1887 newreact(
"H,e-=>H-,PHOTON",
"hmattach",1.,0.,0.);
1888 newreact(
"H-,H+=>H2+,e-",
"hmihph2p",1.,0.,0.);
1889 newreact(
"H-,e-=>H,e-,e-",
"cionhm",1.,0.,0.);
1890 newreact(
"H,e-,e-=>H-,e-",
"c3bod",1.,0.,0.);
1891 newreact(
"H,H-=>H2,e-",
"asdfg",1.,0.,0.);
1892 newreact(
"H,H-=>H2*,e-",
"asdfs",1.,0.,0.);
1893 newreact(
"H2,e-=>H,H-",
"asdbg",1.,0.,0.);
1894 newreact(
"H2*,e-=>H,H-",
"asdbs",1.,0.,0.);
1895 newreact(
"H-,H+=>H,H",
"hneut",1.,0.,0.);
1896 newreact(
"H,H=>H-,H+",
"bhneut",1.,0.,0.);
1901 newreact(
"H2*,PHOTON=>H,H,PHOTON",
"h2_ind_rad_diss",1.,0.,0.);
1902 newreact(
"H,H2+=>H+,H2",
"hmrate",0.,0.,0.);
1907 newreact(
"H2,PHOTON=>H,H,PHOTON",
"h2_ind_rad_diss",1.,0.,0.);
1908 newreact(
"H,H2+=>H+,H2",
"hmrate",6.4e-10f,0.,0.);
1913 newreact(
"H,H=>H2,PHOTON",
"radasc",3e-14,0.,0.);
1914 newreact(
"H,H=>H2+,e-",
"assoc_ion",3.27e-10,-0.35,17829.);
1915 newreact(
"H2,H=>H,H,H",
"rh2g_dis_h",1.,0.,0.);
1916 newreact(
"H2,H2=>H,H,H2",
"rh2g_dis_h2",1.,0.,0.);
1917 newreact(
"H,H,H=>H2,H",
"bh2g_dis_h",1.,0.,0.);
1918 newreact(
"H,H,H2=>H2,H2",
"bh2g_dis_h2",1.,0.,0.);
1920 newreact(
"H2,PHOTON=>H2+,e-",
"h2photon",1.,0.,0.);
1921 newreact(
"H2,CRPHOT=>H2+,e-",
"h2crphot",1.,0.,0.);
1922 newreact(
"H2*,PHOTON=>H2+,e-",
"h2photon",1.,0.,0.);
1923 newreact(
"H2*,CRPHOT=>H2+,e-",
"h2crphot",1.,0.,0.);
1924 newreact(
"H2,CRPHOT=>H,H",
"h2crphh",1.,0.,0.);
1925 newreact(
"H2,CRPHOT=>H+,H-",
"cryden_ov_bg",3.9e-21,0.,0.);
1926 newreact(
"H2*,CRPHOT=>H+,H-",
"cryden_ov_bg",3.9e-21,0.,0.);
1927 newreact(
"H2,CRPHOT=>H+,H,e-",
"cryden_ov_bg",2.2e-19,0.,0.);
1928 newreact(
"H2*,CRPHOT=>H+,H,e-",
"cryden_ov_bg",2.2e-19,0.,0.);
1929 newreact(
"H+,H=>H2+,PHOTON",
"radath",1.,0.,0.);
1930 newreact(
"H2+,PHOTON=>H,H+",
"gamtwo",1.,0.,0.);
1931 newreact(
"H2+,CRPHOT=>H,H+",
"hlike_phot",1.,0.,0.);
1934 newreact(
"H3+,CRPHOT=>H2+,H+,e-",
"hlike_phot",2.,0.,0.);
1935 newreact(
"H,H+,e-=>H2+,e-",
"hmrate",2e-7 * SAHA * (4./(2.*1.)) * 3.634e-5 * pow(300.,-1.50),-1.50,0.);
1939 newreact(
"H2,CRPHOT=>H2+,e-",
"hmrate",4.4e-17,0.,0.);
1940 newreact(
"H2,CRPHOT=>H,H+,e-",
"hmrate",1e-19,0.,0.);
1941 newreact(
"H2*,CRPHOT=>H,H+,e-",
"hmrate",1e-19,0.,0.);
1942 newreact(
"H2*,CRPHOT=>H2+,e-",
"hmrate",4.4e-17,0.,0.);
1943 newreact(
"H2,CRPHOT=>H,H",
"hmrate",5e-18,0.,0.);
1944 newreact(
"e-,H3+=>H2,H",
"hmrate",2.5e-8,-0.3,0.);
1945 newreact(
"e-,H3+=>H,H,H",
"hmrate",7.5e-8,-0.3,0.);
1946 newreact(
"H+,H=>H2+,PHOTON",
"hmrate",5.3e-19,1.85,0);
1950 newreact(
"H2+,e-=>H,H",
"hmrate",1.6e-8,-0.43,0.);
1951 newreact(
"H2+,PHOTON=>H,H+",
"th85rate",5.7e-10,0.,1.9);
1954 newreact(
"H2*,CRPHOT=>H,H",
"h2scrphh",1.,0.,0.);
1955 newreact(
"H2,H2+=>H,H3+",
"h2ph3p",1.,0.,0.);
1956 newreact(
"H2*=>H2,PHOTON",
"h2s_sp_decay",1.,0.,0.);
1957 newreact(
"H2*,H2=>H2,H2",
"h2_collh2_deexc",1.,0.,0.);
1958 newreact(
"H2*,H=>H2,H",
"h2_collh_deexc",1.,0.,0.);
1959 newreact(
"H2,H2=>H2*,H2",
"h2_collh2_excit",1.,0.,0.);
1960 newreact(
"H2,H=>H2*,H",
"h2_collh_excit",1.,0.,0.);
1963 newreact(
"H2*,H=>H,H,H",
"rh2s_dis_h",1.,0.,0.);
1964 newreact(
"H2*,H2=>H2,H,H",
"rh2s_dis_h2",1.,0.,0.);
1965 newreact(
"H2*,H2*=>H2,H,H",
"rh2s_dis_h2",1.,0.,0.);
1966 newreact(
"H2*,H2*=>H2*,H,H",
"rh2s_dis_h2_nodeexcit",1.,0.,0.);
1970 newreact(
"H2,He=>He,H,H",
"rh2g_dis_h",1.,0.,0.);
1971 newreact(
"H2,H+=>H+,H,H",
"rh2g_dis_h",1.,0.,0.);
1972 newreact(
"H2,H3+=>H3+,H,H",
"rh2g_dis_h",1.,0.,0.);
1973 newreact(
"H2,e-=>e-,H,H",
"rh2g_dis_h",1.,0.,0.);
1974 newreact(
"H2*,He=>He,H,H",
"rh2s_dis_h",1.,0.,0.);
1975 newreact(
"H2*,H+=>H+,H,H",
"rh2s_dis_h",1.,0.,0.);
1976 newreact(
"H2*,H3+=>H3+,H,H",
"rh2s_dis_h",1.,0.,0.);
1977 newreact(
"H2*,e-=>e-,H,H",
"rh2s_dis_h",1.,0.,0.);
1980 newreact(
"He,H,H=>H2,He",
"bh2g_dis_h",1.,0.,0.);
1981 newreact(
"H+,H,H=>H2,H+",
"bh2g_dis_h",1.,0.,0.);
1982 newreact(
"H3+,H,H=>H2,H3+",
"bh2g_dis_h",1.,0.,0.);
1983 newreact(
"e-,H,H=>H2,e-",
"bh2g_dis_h",1.,0.,0.);
1984 newreact(
"H,H,H=>H2*,H",
"bh2s_dis_h",1.,0.,0.);
1985 newreact(
"H2,H,H=>H2*,H2",
"bh2s_dis_h2",1.,0.,0.);
1986 newreact(
"H2,H,H=>H2*,H2*",
"bh2s_dis_h2",1.,0.,0.);
1987 newreact(
"H2*,H,H=>H2*,H2*",
"bh2s_dis_h2",1.,0.,0.);
1988 newreact(
"He,H,H=>H2*,He",
"bh2s_dis_h",1.,0.,0.);
1989 newreact(
"H+,H,H=>H2*,H+",
"bh2s_dis_h",1.,0.,0.);
1990 newreact(
"H3+,H,H=>H2*,H3+",
"bh2s_dis_h",1.,0.,0.);
1991 newreact(
"e-,H,H=>H2*,e-",
"bh2s_dis_h",1.,0.,0.);
1994 newreact(
"H2,PHOTON=>H2*",
"h2gexcit",1.,0.,0.);
1995 newreact(
"H2*,PHOTON=>H,H",
"h2sdissoc",1.,0.,0.);
1996 newreact(
"H2,PHOTON=>H,H",
"h2gdissoc",1.,0.,0.);
1997 newreact(
"HeH+,PHOTON=>H,He+",
"gamheh",1.,0.,0.);
2002 newreact(
"OHgrn,Hgrn=>H2Ogrn",
"grn_react",1.,0.,0.);
2007 fprintf(stderr,
"Finished testing against UDFA database\n");
2023 it->second->index = index;
2037 DEBUG_ENTRY(
"mole_generate_isotopologue_reactions()" );
2039 bool lgDebug =
false;
2043 vector<string> newReactionStrings;
2044 vector<mole_reaction*> oldReactions;
2045 vector<realnum> branchingRatios;
2050 bool lgFound =
false;
2054 for(
long j=0; j<it->second->nproducts; ++j )
2057 if( it->second->products[j]->nNuclide.find( atomOld ) != it->second->products[j]->nNuclide.end() )
2058 numSites += it->second->products[j]->nNuclide[ atomOld ];
2061 for(
long i=0; i<it->second->nreactants; ++i )
2064 for(
nNucs_i atom=it->second->reactants[i]->nNuclide.begin(); atom != it->second->reactants[i]->nNuclide.end(); ++atom )
2066 if( atom->first->label() == atom_old )
2080 fprintf(
ioQQQ,
"DEBUGGG mole_generate_isotopologue_reactions %s ..........\n", it->first.c_str() );
2082 for(
long i=0; i<it->second->nreactants; ++i )
2085 if( it->second->reactants[i]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
2088 vector<string> react_iso_labels;
2090 vector< int > numAtoms;
2091 string embellishments;
2093 bool lgParseOK =
parse_species_label( it->second->reactants[i]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
2104 for(
unsigned j=0; j<react_iso_labels.size(); ++j )
2106 int numAtomsTot = 0;
2107 for(
long k=0; k<it->second->nproducts; ++k )
2110 if( it->second->products[k]->mole_mass < 0.99 * ATOMIC_MASS_UNIT )
2113 atomsLeftToRight.clear();
2115 embellishments.clear();
2117 lgParseOK =
parse_species_label( it->second->products[k]->label.c_str(), atomsLeftToRight, numAtoms, embellishments );
2118 ASSERT( lgParseOK ==
true );
2120 for(
unsigned position = 0; position < atomsLeftToRight.size(); ++position )
2122 string prod_iso_label;
2133 if( prod_iso_label.empty() )
2137 string react_string;
2139 for(
long i1=0; i1<i; ++i1 )
2141 react_string += it->second->reactants[i1]->label;
2142 if( i1 != it->second->nreactants-1 )
2143 react_string +=
",";
2145 react_string += react_iso_labels[j];
2146 if( i != it->second->nreactants-1 )
2147 react_string +=
",";
2148 for(
long i2=i+1; i2<it->second->nreactants; ++i2 )
2150 react_string += it->second->reactants[i2]->label;
2151 if( i2 != it->second->nreactants-1 )
2152 react_string +=
",";
2155 react_string +=
"=>";
2157 for(
long k1=0; k1<k; ++k1 )
2159 react_string += it->second->products[k1]->label;
2160 if( k1 != it->second->nproducts-1 )
2161 react_string +=
",";
2163 react_string += prod_iso_label;
2164 if( k != it->second->nproducts-1 )
2165 react_string +=
",";
2166 for(
long k2=k+1; k2<it->second->nproducts; ++k2 )
2168 react_string += it->second->products[k2]->label;
2169 if( k2 != it->second->nproducts-1 )
2170 react_string +=
",";
2175 newReactionStrings.push_back( canon_react_string );
2176 oldReactions.push_back( it->second.get_ptr() );
2178 branchingRatios.push_back( numAtoms[position]/numSites );
2181 fprintf(
ioQQQ,
"DEBUGGG mole_generate_isotopologue_reactions .................... %s\t\t(%2i/%2i).\n",
2182 canon_react_string.c_str(), numAtoms[position], numSites );
2184 numAtomsTot += numAtoms[position];
2187 ASSERT( numAtomsTot == numSites );
2192 ASSERT( oldReactions.size() == newReactionStrings.size() );
2195 vector<mole_reaction*>::const_iterator it2 = oldReactions.begin();
2196 vector<realnum>::const_iterator it3 = branchingRatios.begin();
2197 for( vector<string>::const_iterator it1 = newReactionStrings.begin(); it1 != newReactionStrings.end(); ++it1, ++it2, ++it3 )
2199 fixit(
"make adjustments to a for mass?");
2204 ASSERT( *it3 <= 1. + FLT_EPSILON );
2205 fixit(
"multiply by *it3 here. ASSERT at mole_reactions.cpp:1190 will blow currently.");
2206 newreact( it1->c_str(), (*it2)->name(), (*it2)->a , (*it2)->b, (*it2)->c );
2217 char chLabel[50], chLabelSave[50];
2224 strcpy( chLabel, p->second->label.c_str() );
2225 strcpy( chLabelSave, p->second->label.c_str() );
2226 char *str = chLabel;
2227 const char *delim =
"=>";
2228 char *chNewProducts = strtok( str, delim );
2229 char *chNewReactants = strtok( NULL, delim );
2230 char chNewReaction[50];
2232 strcpy( chNewReaction, chNewReactants );
2233 strcat( chNewReaction,
"=>" );
2234 strcat( chNewReaction, chNewProducts );
2244 fprintf(
ioQQQ,
"Warning! No reverse reaction for %30s. ", p->second->label.c_str() );
2248 fixit(
"NB reverse reactions should be generated here");
2257 DEBUG_ENTRY(
"mole_get_equilibrium_condition()" );
2266 DEBUG_ENTRY(
"mole_get_equilibrium_condition()" );
2272 double ln_result = 0.;
2278 ln_result += log(fac);
2285 ln_result -= log(fac);
2298 if( sp->
label ==
"PHOTON" || sp->
label ==
"CRPHOT" )
2300 fixit(
"How can we adapt existing structures to have a photon energy or range available here?");
2301 fixit(
"include 2hnu^3/c^2. Understand HNU3C2 macro!");
2304 else if( sp->
label ==
"CRP" || sp->
label ==
"grn" )
2307 fixit(
"need to figure out stat weight for any given particle");
2310 double deltaH = sp->
form_enthalpy * (1e10/AVOGADRO/BOLTZMANN);
2314 ASSERT( part_fun >= 0. );
2323 vector<long> numForm, numDest;
2341 for(
unsigned i=0; i<numForm.size(); ++i )
2343 if( numForm[i]==0 && numDest[i]==0 )
2345 if( numForm[i]>1 && numDest[i]>1 )
2358 char *label, *reactstr, *f;
2373 STATIC void newreact(
const char label[],
const char fun[],
double a,
double b,
double c)
2377 ratefunc rate = findfunc(fun);
2378 if (rate.get_ptr() == NULL)
2380 fprintf(stderr,
"Rate function %s not known for reaction %s. Aborting. Sorry.\n",fun,label);
2386 fprintf(
ioQQQ,
"\n PROBLEM Reaction %s has negative pre-coefficient. Aborting. Sorry.\n", label);
2390 rate->label = label;
2394 rate->source = source;
2409 fprintf(
ioQQQ,
" W-reaction %s disabled\n",rate->label.c_str());
2414 const char *rateLabelPtr = rate->label.c_str();
2418 rate->udfastate =
ABSENT;
2422 if (rate->udfastate ==
ABSENT)
2424 fprintf(stderr,
"Reaction %s not in UDFA\n",rateLabelPtr);
2429 if (rate->nreactants == 2 && rate->reactants[0]->mole_mass!=0.0 && rate->reactants[1]->mole_mass!=0.0 )
2431 rate->reduced_mass = 1./(1./rate->reactants[0]->mole_mass+1./rate->reactants[1]->mole_mass);
2435 rate->reduced_mass = 0.;
2445 fprintf(
ioQQQ,
"Attention: duplicate reaction %s -- using new version\n",rateLabelPtr);
2466 bool lgProd =
false;
2468 for(
int i=0;!i || label[i-1]!=
'\0';i++)
2470 if(label[i] ==
',' || label[i] ==
'=' || label[i] ==
'\0')
2476 fprintf(
ioQQQ,
"Mole_reactions: ignoring reaction %s (species %s not active)\n",label,buf.c_str());
2484 fprintf(stderr,
"Mole_reactions: Too many reactants in %s, only %d allowed\n",label,MAXREACTANTS);
2494 fprintf(stderr,
"Mole_reactions: Too many products in %s, only %d allowed\n",label,MAXPRODUCTS);
2503 if (label[i] !=
'>')
2535 ratefunc rate = findfunc(
"null");
2536 rate->label = label;
2573 rate->
label = newLabel;
2583 bool lgTrivial =
false;
2728 fprintf(sparsefp,
"P4\n%d %d\n",
2736 ch = (ch << 1) | (c[i][j] != 0.0);
2759 int dcharge, n,
sign;
2766 nel[it->first] += it->second;
2772 nel[it->first] -= it->second;
2777 fprintf(stderr,
"Reaction %s charge out of balance by %d\n",
2778 rate->
label.c_str(),dcharge);
2782 for(
nNucs_i it = nel.begin(); it != nel.end(); ++it )
2790 fprintf(stderr,
"Error: reaction %s %s %d of element %s\n",
2791 rate->
label.c_str(),sign==1?
"destroys":
"creates",
2793 it->first->label().c_str() );
2805 class formula_species {
2810 bool operator< (
const formula_species &a,
const formula_species &b)
2815 if (a.reactants[i]<b.reactants[i])
2817 else if (a.reactants[i]>b.reactants[i])
2822 if (a.products[i]<b.products[i])
2824 else if (a.products[i]>b.products[i])
2830 class udfa_reaction {
2835 double a, b, c, trange[2];
2839 static map <formula_species,count_ptr<udfa_reaction> >
udfatab;
2849 fprintf(stderr,
"Error, could not read %s\n",file);
2853 fixit(
"this seg-faults if file ends in blank line!");
2862 #define FLTEQ(a,b) (fabs((a)-(b)) <= 1e-6*fabs((a)+(b)))
2866 unsigned int havespecies = 1, i, n;
2878 r->l.reactants[n] = NULL;
2890 if (!strncmp(f,
"CRPHOT",6))
2899 r->l.products[n] = NULL;
2937 r->source = f[0]?f[0]:
'?';
2940 r->trange[n] = atof(f);
2945 if (udfatab.find(r->l) != udfatab.end())
2947 fprintf(stderr,
"Duplicate reaction\n");
2976 map<formula_species, count_ptr<udfa_reaction> >::iterator p
2978 if (p == udfatab.end() )
3002 ratefunc findfunc (
const char name[])
3010 enum { DEBUG_MOLE =
false };
3022 long index = rate.
index;
3027 if (fabs(newrk-oldrk) > 0.1*newrk)
3029 rate.
label.c_str(),oldrk,newrk);
3037 enum { DEBUG_MOLE =
false };
3048 double bigchange = 0.;
3049 unsigned long bigindex = ULONG_MAX;
3054 sum = oldrk+newrk, diff = newrk-oldrk;
3057 double change = fabs(diff)/sum;
3058 if (change > bigchange)
3070 if (bigindex == (
unsigned) rate.
index)
3076 pct = 100.*(newrk-oldrk)/oldrk;
3077 fprintf(
ioQQQ,
"Zone %ld, big chemistry rate change %s:"
3078 " %15.8g => %15.8g (%6.2g%%)\n",
3096 double T_ortho_para_crit,xi_ELRD, beta_alpha_ELRD, recombination_efficiency_ELRD, Td;
3108 fixit(
"Is this still necessary?");
3121 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3131 vector<double> qtemp, qprob;
3138 qheat(qtemp,qprob,&qnbin,nd);
3140 if(
gv.
bin[nd]->lgUseQHeat )
3148 qtemp[0] =
gv.
bin[nd]->tedust;
3151 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.;
3153 for( k=0; k < qnbin; k++ )
3163 double conversion_efficiency_HM79 = 1/(1. + 1e4*
sexp(920./qtemp[k]));
3166 gv.
bin[nd]->rate_h2_form_grains_HM79 += qprob[k] * sticking_probability_H *
3167 conversion_efficiency_HM79;
3173 gv.
bin[nd]->rate_h2_form_grains_HM79 *= 0.5 * AveVelH *
3174 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3176 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_HM79 > 0. );
3188 double conversion_efficiency_HM79 = 1/(1. + 1e4*
sexp(920./
gv.
bin[nd]->tedust));
3193 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3195 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * conversion_efficiency_HM79;
3196 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_HM79 > 0. );
3207 double sqrt_term = 35.399494936611667;
3209 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.;
3211 for( k=0; k < qnbin; k++ )
3213 double beta_alpha = 0.25 * sqrt_term *
sexp(200./qtemp[k] );
3215 double xi = 1./ (1. + 1.3e13*
sexp(1.5*1e4/qtemp[k])*sqrt_term/(2.*f) );
3217 double beta = 3e12 *
sexp( 320. / qtemp[k] );
3220 double recombination_efficiency_CT02 = xi / (1. + 0.005*f/2./
SDIV(beta) + beta_alpha );
3226 gv.
bin[nd]->rate_h2_form_grains_CT02 += qprob[k] * sticking_probability_H *
3227 recombination_efficiency_CT02;
3234 gv.
bin[nd]->rate_h2_form_grains_CT02 *= 0.5 * AveVelH *
3235 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3237 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_CT02 > 0. );
3247 double sqrt_term = 35.399494936611667;
3248 double beta_alpha = 0.25 * sqrt_term * exp(-200./
gv.
bin[nd]->tedust);
3250 double xi = 1./ (1. + 1.3e13*exp(-1.5e4/
gv.
bin[nd]->tedust)*sqrt_term/(2.*f) );
3252 double beta = 3e12 * exp(-320./
gv.
bin[nd]->tedust);
3255 double recombination_efficiency_CT02 = beta*xi / (beta + 0.0025*f + beta*beta_alpha);
3261 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3262 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_CT02;
3263 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_CT02 >= 0. );
3276 gv.
bin[nd]->rate_h2_form_grains_ELRD= 0.;
3282 for( k=0; k < qnbin; k++ )
3286 beta_alpha_ELRD = exp(-800./Td) / (0.5389970511202651 * exp(-540./Td) + 5.6333909478365e-14*sqrt(Td) ) ;
3287 xi_ELRD= 1./(1. + 4.61e24*
sexp(45000./Td));
3288 recombination_efficiency_ELRD= (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3289 sticking_probability_H = 1./(1. + 0.04*sqrt(Td+
phycon.
te) +
3292 gv.
bin[nd]->rate_h2_form_grains_ELRD+= qprob[k] * sticking_probability_H *
3293 recombination_efficiency_ELRD;
3301 gv.
bin[nd]->rate_h2_form_grains_ELRD*= 0.5 * AveVelH *
3302 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3304 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD> 0. );
3309 for( k=0; k < qnbin; k++ )
3313 beta_alpha_ELRD = exp(-450./Td) / (0.4266153643741715*exp(-340./Td) + 2.5335919594255e-14*sqrt(Td) ) ;
3314 xi_ELRD= 1./(1. + 7.00e24*
sexp(45000./Td));
3315 recombination_efficiency_ELRD= (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3316 sticking_probability_H = 1./(1. + 0.04*sqrt(Td+
phycon.
te) +
3319 gv.
bin[nd]->rate_h2_form_grains_ELRD+= qprob[k] * sticking_probability_H *
3320 recombination_efficiency_ELRD;
3327 gv.
bin[nd]->rate_h2_form_grains_ELRD*= 0.5 * AveVelH *
3328 gv.
bin[nd]->IntArea/4. *
gv.
bin[nd]->cnv_H_pCM3;
3330 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD> 0. );
3340 gv.
bin[nd]->rate_h2_form_grains_ELRD= 0.;
3346 Td =
gv.
bin[nd]->tedust;
3348 beta_alpha_ELRD = exp(-800./Td) / (0.5389970511202651 * exp(-540./Td) + 5.6333909478365e-14*sqrt(Td) ) ;
3349 xi_ELRD= 1./(1. + 4.61e24*
sexp(45000./Td));
3350 recombination_efficiency_ELRD = (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3356 gv.
bin[nd]->rate_h2_form_grains_ELRD = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3357 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_ELRD;
3359 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD > 0. );
3364 Td =
gv.
bin[nd]->tedust;
3366 beta_alpha_ELRD = exp(-450./Td) / (0.4266153643741715*exp(-340./Td) + 2.5335919594255e-14*sqrt(Td) ) ;
3367 xi_ELRD= 1./(1. + 7.00e24*
sexp(45000./Td));
3368 recombination_efficiency_ELRD = (1. / (1. + beta_alpha_ELRD))*xi_ELRD;
3374 gv.
bin[nd]->rate_h2_form_grains_ELRD = 0.5 * AveVelH *
gv.
bin[nd]->IntArea/4. *
3375 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * recombination_efficiency_ELRD;
3377 ASSERT(
gv.
bin[nd]->rate_h2_form_grains_ELRD > 0. );
3398 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H;
3424 if(
gv.
bin[nd]->tedust < T_ortho_para_crit )
3430 gv.
bin[nd]->cnv_H_pCM3 * sticking_probability_H * efficiency_opr;
3449 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3451 gv.
bin[nd]->rate_h2_form_grains_CT02 = 0.;
3452 gv.
bin[nd]->rate_h2_form_grains_HM79 = 0.;
3465 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
3472 gv.
bin[nd]->rate_h2_form_grains_used =
3482 gv.
bin[nd]->rate_h2_form_grains_used =
3490 gv.
bin[nd]->rate_h2_form_grains_used =
3498 gv.
bin[nd]->rate_h2_form_grains_used =
3528 static double teused=-1;
3538 cr_H2dis_ELWERT_H2g,
3539 cr_H2dis_ELWERT_H2s;
3646 enum {DEBUG_LOC=
false};
3661 fixit(
"Wasted cycles if we don't use Stancil's rates below? Why not put this down there/if used?");
3665 (*diatom)->Mol_Photo_Diss_Rates();
3697 enum {DEBUG_LOC=
false};
3715 enum {DEBUG_LOC=
false};
3717 if( DEBUG_LOC &&
nzone>400)
3826 double sqrtx = sqrt(1.+x);
3830 double fshield = 0.965/
POW2(1.+x/b5) + 0.035/sqrtx *
3882 double x_H2g, x_H2s,
3883 fshield_H2g, fshield_H2s,
3885 static double a_H2g, a_H2s,
3918 a_H2g = 0.06 * log_G0_face + 1.32;
3919 a_H2s = 0.375 * log_G0_face + 2.125;
3921 e1_H2g = -0.05 * log_G0_face + 2.25;
3922 e1_H2s = -0.125 * log_G0_face + 2.625;
3924 e2_H2g = -0.005 * log_G0_face + 0.625;
3926 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
3933 k_f_H2s =
MAX2(0.1,2.375 * log_G0_face - 1.875 );
3936 k_H2g_to_H2s =
MAX2(1.,-1.75 * log_G0_face + 11.25);
3951 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
3952 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
4116 enum {DEBUG_LOC=
false};
4120 fprintf(
ioQQQ,
" Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",
4143 return &(*p->second);
4219 for(
int i=0;i<rate.
nreactants && ipthis == -1;i++)
4228 double ratevi = rate.
a * rate.
rk();
4272 int ipspin = 0, ipfreein = 0;
4280 int ipspout = 0, ipfreeout = 0;
4290 int newsp = ipspout-ipspin;
4295 int nbondsbroken = ipfreeout-ipfreein;
4296 if (nbondsbroken <= 0)
4299 double fracbroken = nbondsbroken/((double)ipfreeout);
4300 ASSERT( fracbroken <= 1.0 );
4311 double ratesp = ratevi*newsp;
4313 ratesp *= fracbroken;
4347 double ratevi = rate.
a * rate.
rk();
4352 ratev += ipthis*ratevi;
4374 double heating = 0.;
4375 map<double,string> heatMap;
4387 bool lgCanSkip =
false;
4415 realnum reaction_enthalpy = 0.;
4424 for(
long i=0; i < rate.
nproducts; ++i )
4434 double heat = reaction_enthalpy*rate_tot*(1e10/AVOGADRO);
4435 heatMap[heat] = rate.
label;
4445 for( map<double,string>::reverse_iterator it = heatMap.rbegin(); it != heatMap.rend(); ++it, ++index )
4447 fprintf(
ioQQQ,
"DEBUGGG heat %li\t%li\t%.6e\t%s\n", index,
nzone, it->first, it->second.c_str() );
4452 for( map<double,string>::iterator it = heatMap.begin(); it != heatMap.end(); ++it, ++index )
4454 if( it->first >= 0. )
4456 fprintf(
ioQQQ,
"DEBUGGG cool %li\t%li\t%.6e\t%s\n", index,
nzone, it->first, it->second.c_str() );
4479 double T2 = T_gas / 100.;
4480 double Tg2 = T_grain / 100.;
4481 double S = 1./(1. + 0.4*sqrt(Tg2 + T2) + 0.2*T2 + 0.08*T2*T2);
molecule * reactants[MAXREACTANTS]
double sink_rate_tot(const char chSpecies[]) const
t_mole_global mole_global
STATIC void register_reaction_vectors(count_ptr< mole_reaction > rate)
double H2_Solomon_dissoc_rate_used_H2g
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
vector< double > reaction_rks
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
STATIC double sticking_probability_H_func(double T_gas, double T_grain)
molecule::nNucsMap::iterator nNucs_i
void mole_create_react(void)
double H2_photodissoc_ELWERT_H2s
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
NORETURN void TotalInsanity(void)
bool lgH2_Chemistry_BigH2
double Average_collH2_excit
STATIC void mole_check_reverse_reactions(void)
double H2_H2g_to_H2s_rate_used
realnum UV_Cont_rel2_Draine_DB96_face
virtual const char * name()=0
bool lgLeiden_Keep_ipMH2s
double H2star_forms_grains
bool operator<(const count_ptr< T > &a, const count_ptr< T > &b)
double H2_H2g_to_H2s_rate_BHT90
double Cont_Dissoc_Rate_H2g
double H2_photodissoc_used_H2s
double H2_Solomon_dissoc_rate_TH85_H2s
double H2_photodissoc_BHT90
STATIC void mole_h2_grain_form(void)
double findrk(const char buf[]) const
map< string, count_ptr< mole_reaction > > reactab
double rate_grain_op_conserve
double H2_photodissoc_TH85
double H2_Solomon_dissoc_rate_BD96_H2g
sys_float sexp(sys_float x)
double H2_Solomon_dissoc_rate_ELWERT_H2s
double H2_H2g_to_H2s_rate_ELWERT
molezone * findspecieslocal(const char buf[])
double source_rate_tot(const char chSpecies[]) const
STATIC char * getcsvfield(char **s, char c)
double rate_h2_form_grains_set
void iso_photo(long ipISO, long nelem)
ChemNuclideList nuclide_list
double H2_Solomon_dissoc_rate_BD96_H2s
bool parse_species_label(const char label[], ChemNuclideList &atomsLeftToRight, vector< int > &numAtoms, string &embellishments)
double anu(size_t i) const
double Average_collH2_dissoc_g
double Average_collH2_deexcit
map< string, bool > offReactions
double H2_Solomon_dissoc_rate_TH85_H2g
STATIC void newreact(const char label[], const char fun[], double a, double b, double c)
static t_version & Inst()
double Solomon_dissoc_rate_g
t_iso_sp iso_sp[NISO][LIMELM]
double H2_Solomon_dissoc_rate_BHT90_H2g
double xIonDense[LIMELM][LIMELM+1]
double Average_collH2_dissoc_s
bool fp_equal(sys_float x, sys_float y, int n=3)
double H2_Solomon_dissoc_rate_used_H2s
STATIC long parse_reaction(count_ptr< mole_reaction > &rate, const char label[])
molecule * products[MAXPRODUCTS]
void mole_cmp_num_in_out_reactions(void)
STATIC void mole_h_reactions(void)
double chem_heat(void) const
double H2_photodissoc_ELWERT_H2g
STATIC bool lgReactionTrivial(count_ptr< mole_reaction > &rate)
vector< double > old_reaction_rks
molecule * rvector[MAXREACTANTS]
const bool ENABLE_QUANTUM_HEATING
void qheat(vector< double > &, vector< double > &, long *, size_t)
double photodissoc_BigH2_H2s
double H2_H2g_to_H2s_rate_BD96
map< string, count_ptr< mole_reaction > >::iterator mole_reaction_i
molecule * findspecies(const char buf[])
double Average_collH_deexcit
void create_isotopologues_one_position(unsigned position, ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, string &newLabel)
valarray< class molezone > species
virtual double rk() const =0
double sink_rate(const molecule *const sp, const mole_reaction &rate) const
double photodissoc_BigH2_H2g
void mole_rk_bigchange(void)
count_ptr< chem_nuclide > findnuclide(const char buf[])
vector< diatomics * > diatoms
double powi(double, long int)
STATIC bool lgReactBalance(const count_ptr< mole_reaction > &rate)
map< string, count_ptr< mole_reaction > > functab
double H2_Solomon_dissoc_rate_BHT90_H2s
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
realnum GetDopplerWidth(realnum massAMU)
bool exists(const molecule *m)
double column(const genericState &gs)
map< const count_ptr< chem_nuclide >, int, element_pointer_value_less > nNucsMap
realnum GetAveVelocity(realnum massAMU)
realnum UV_Cont_rel2_Habing_TH85_face
double H2_photodissoc_used_H2g
double esc_PRD_1side(double tau, double a)
double **** PhotoRate_Shell
double rate_grain_J1_to_J0
realnum gas_phase[LIMELM]
realnum UV_Cont_rel2_Draine_DB96_depth
vector< count_ptr< chem_nuclide > > ChemNuclideList
realnum AtomicWeight[LIMELM]
STATIC void mole_generate_isotopologue_reactions(string atom_old, string atom_new)
double frac_H2star_hminus()
double findrate(const char buf[]) const
molecule * pvector[MAXPRODUCTS]
map< string, count_ptr< mole_reaction > >::const_iterator mole_reaction_ci
STATIC void compare_udfa(const count_ptr< mole_reaction > &rate)
bool lgH2_grain_deexcitation
char chH2_small_model_type
bool isMonatomic(void) const
double Average_collH_excit
#define DEBUG_ENTRY(funcname)
realnum UV_Cont_rel2_Habing_spec_depth
double powpq(double x, int p, int q)
double Solomon_dissoc_rate_s
STATIC double sticking_probability_H_HM79(double T_gas, double T_grain)
bool lgDifferByExcitation(const molecule &mol1, const molecule &mol2)
STATIC void plot_sparsity(void)
double Cont_Dissoc_Rate_H2s
double GammaBn(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double thresh, double *ainduc, double *rcool, t_phoHeat *photoHeat)
void mole_update_rks(void)
diatomics hd("hd", 4100.,&hmi.HD_total, Yan_H2_CS)
STATIC double mole_get_equilibrium_condition(const char buf[])
STATIC void read_data(const char file[], void(*parse)(char *s))
STATIC double mole_partition_function(const molecule *const sp)
realnum UV_Cont_rel2_Habing_TH85_depth
int fprintf(const Output &stream, const char *format,...)
double dissoc_rate(const char chSpecies[]) const
mole_reaction * mole_findrate_s(const char buf[])
sys_float SDIV(sys_float x)
molecule * rvector_excit[MAXREACTANTS]
double H2_H2g_to_H2s_rate_TH85
STATIC void parse_udfa(char *s)
virtual mole_reaction * Create() const =0
double Average_collH_dissoc_g
double hmrate4(double a, double b, double c, double te)
double esca0k2(double taume)
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
STATIC void canonicalize_reaction(count_ptr< mole_reaction > &rate)
STATIC string canonicalize_reaction_label(const char label[])
double Average_collH_dissoc_s
t_secondaries secondaries
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
double H2_Solomon_dissoc_rate_ELWERT_H2g
double HMinus_induc_rec_rate
STATIC void parse_base(char *s)
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
void create_isotopologues(ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, vector< string > &newLabels)
vector< diatomics * >::iterator diatom_iter
double rate_h2_form_grains_used_total
static map< formula_species, count_ptr< udfa_reaction > > udfatab
molecule * pvector_excit[MAXPRODUCTS]
double HMinus_induc_rec_cooling
double H2star_forms_hminus
1.8.5