-
Notifications
You must be signed in to change notification settings - Fork 17
/
h2collisions.c
115 lines (77 loc) · 3.29 KB
/
h2collisions.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
/* ------- file: -------------------------- h2collisions.c ----------
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Sat Sep 19 15:27:09 2009 --
-------------------------- ----------RH-- */
/* --- Calculate collisional excitation rates for the vibrational
transitions in the CO ground electronic state X^1\Sigma^+.
Calculated are the collisional coefficients Omega_ul [m^3 s^-1],
where C_ul = Omega_ul * n_pert, and n_pert is the number density
of perturbers [m^-3].
See: Ayres, T.R., Wiedemann, G., 1989, ApJ 338, 1033
-- -------------- */
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "rh.h"
#include "atom.h"
#include "atmos.h"
#include "error.h"
#include "constant.h"
/* --- Collisional constants in Landau - Teller format -- ----------- */
#define A_H 3.0
#define B_H 18.1
#define A_HE 87.0
#define B_HE 19.1
#define A_H2 64.0
#define B_H2 19.1
/* --- Harmonic oscillator constant for H2 -- -------------- */
#define OMEGA0 2.103E+05
/* --- Function prototypes -- -------------- */
/* --- Global variables -- -------------- */
extern Atmosphere atmos;
extern char messageStr[];
/* ------- begin -------------------------- H2collisions.c ---------- */
void H2collisions(struct Molecule *molecule)
{
const char routineName[] = "H2collisions";
register int k;
double hcomega_k, *beta, C_0;
if (!strstr(molecule->ID, "H2")) {
sprintf(messageStr, "Molecule is not H2: %s\n", molecule->ID);
Error(ERROR_LEVEL_2, routineName, messageStr);
}
C_0 = KBOLTZMANN / ATM_TO_PA;
hcomega_k = HPLANCK * CLIGHT * OMEGA0 / KBOLTZMANN;
beta = (double *) malloc(atmos.Nspace * sizeof(double));
for (k = 0; k < atmos.Nspace; k++)
beta[k] = hcomega_k / atmos.T[k];
/* --- For now only calculate a v-independent collisional
de-excitation rate C_ul. The upward rate is then given by:
C_lu = nv^*_u / nv^*_l * C_ul
-- -------------- */
molecule->C_ul = (double *) malloc(atmos.Nspace * sizeof(double));
for (k = 0; k < atmos.Nspace; k++) {
if (molecule->n[k]) {
/* --- Neutral hydrogen contribution -- -------------- */
molecule->C_ul[k] = C_0 * atmos.T[k] *
exp(B_H - A_H * pow(atmos.T[k], -0.33333333)) /
(1.0 - exp(-beta[k])) * atmos.H->n[0][k];
/* --- Neutral helium contribution
(assuming all helium is neutral) -- -------------- */
molecule->C_ul[k] += C_0 * atmos.T[k] *
exp(B_HE - A_HE * pow(atmos.T[k], -0.33333333)) /
(1.0 - exp(-beta[k])) * atmos.elements[1].abund * atmos.nHtot[k];
/* --- H2 contribution -- -------------- */
molecule->C_ul[k] += C_0 * atmos.T[k] *
exp(B_H2 - A_H2 * pow(atmos.T[k], -0.33333333)) /
(1.0 - exp(-beta[k])) * atmos.H2->n[k];
/* --- Electronic contribution -- -------------- */
molecule->C_ul[k] += (1.4E-09 * CUBE(CM_TO_M)) / sqrt(beta[k]) *
((1.0 + beta[k]) + 19.0*exp(-3.22*beta[k]) *
(1.0 + 4.22*beta[k])) * atmos.ne[k];
}
}
free(beta);
}
/* ------- end ---------------------------- H2collisions.c ---------- */