cloudy trunk
Loading...
Searching...
No Matches
atom_pop5.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/*atom_pop5 do five level atom population and cooling */
4#include "cddefines.h"
5#include "physconst.h"
6#include "phycon.h"
7#include "thermal.h"
8#include "dense.h"
9#include "thirdparty.h"
10#include "atoms.h"
11
12/*atom_pop5 do five level atom population and cooling */
14 /* vector giving statistical weights on the five levels */
15 const double g[],
16 /* vector giving the excitation energy differences of the 5 levels. The energies
17 * are the energy in wavenumbers between adjacent levels. So EnerWN[0] is the energy
18 * 1-2, EnerWN[1] is the energy between 2-3, etc */
19 const double EnerWN[],
20 /* the collision strengths for the levels */
21 double cs12,
22 double cs13,
23 double cs14,
24 double cs15,
25 double cs23,
26 double cs24,
27 double cs25,
28 double cs34,
29 double cs35,
30 double cs45,
31 /* the transition probabilities between the various levels */
32 double a21,
33 double a31,
34 double a41,
35 double a51,
36 double a32,
37 double a42,
38 double a52,
39 double a43,
40 double a53,
41 double a54,
42 /* the destroyed level populations (units cm-3) for the five levels */
43 double p[],
44 /* the total density of this ion */
46 double *Cooling,
47 double *CoolingDeriv,
48 double pump01,
49 double pump02,
50 double pump03,
51 double pump04
52 )
53{
54
55 DEBUG_ENTRY( "atom_pop5()" );
56
57 /* quit if no species present */
58 ASSERT( abund>=0. );
59 if( abund == 0. )
60 {
61 p[0] = 0.;
62 p[1] = 0.;
63 p[2] = 0.;
64 p[3] = 0.;
65 p[4] = 0.;
66 *Cooling = 0.;
67 *CoolingDeriv = 0.;
68 return;
69 }
70
71 // tf = 1.438 / te, converts energy in wavenumbers into Boltzmann factor
72 double tf = T1CM/phycon.te;
73
74 // define Boltzmann factors
75 double BoltzFac[5][5];
76
77 BoltzFac[0][1] = sexp(EnerWN[0]*tf);
78 BoltzFac[1][2] = sexp(EnerWN[1]*tf);
79 BoltzFac[2][3] = sexp(EnerWN[2]*tf);
80 BoltzFac[3][4] = sexp(EnerWN[3]*tf);
81 BoltzFac[0][2] = BoltzFac[0][1]*BoltzFac[1][2];
82 BoltzFac[0][3] = BoltzFac[0][2]*BoltzFac[2][3];
83 BoltzFac[0][4] = BoltzFac[0][3]*BoltzFac[3][4];
84 BoltzFac[1][3] = BoltzFac[1][2]*BoltzFac[2][3];
85 BoltzFac[1][4] = BoltzFac[1][3]*BoltzFac[3][4];
86 BoltzFac[2][4] = BoltzFac[2][3]*BoltzFac[3][4];
87
88 /* quit it highest level Boltzmann factor too large */
89 if( (BoltzFac[0][4]+pump04) == 0. )
90 {
91 p[0] = 0.;
92 p[1] = 0.;
93 p[2] = 0.;
94 p[3] = 0.;
95 p[4] = 0.;
96 *Cooling = 0.;
97 *CoolingDeriv = 0.;
98 return;
99 }
100
101 // get collision rates, dense.cdsqte is 8.629e-6 / sqrte * eden */
102 // rates units are s^-1
103 double col[5][5];
104 col[1][0] = dense.cdsqte*cs12/g[1];
105 col[0][1] = col[1][0]*g[1]/g[0]*BoltzFac[0][1];
106
107 col[2][0] = dense.cdsqte*cs13/g[2];
108 col[0][2] = col[2][0]*g[2]/g[0]*BoltzFac[0][2];
109
110 col[3][0] = dense.cdsqte*cs14/g[3];
111 col[0][3] = col[3][0]*g[3]/g[0]*BoltzFac[0][3];
112
113 col[4][0] = dense.cdsqte*cs15/g[4];
114 col[0][4] = col[4][0]*g[4]/g[0]*BoltzFac[0][4];
115
116 col[2][1] = dense.cdsqte*cs23/g[2];
117 col[1][2] = col[2][1]*g[2]/g[1]*BoltzFac[1][2];
118
119 col[3][1] = dense.cdsqte*cs24/g[3];
120 col[1][3] = col[3][1]*g[3]/g[1]*BoltzFac[1][3];
121
122 col[4][1] = dense.cdsqte*cs25/g[4];
123 col[1][4] = col[4][1]*g[4]/g[1]*BoltzFac[1][4];
124
125 col[3][2] = dense.cdsqte*cs34/g[3];
126 col[2][3] = col[3][2]*g[3]/g[2]*BoltzFac[2][3];
127
128 col[4][2] = dense.cdsqte*cs35/g[4];
129 col[2][4] = col[4][2]*g[4]/g[2]*BoltzFac[2][4];
130
131 col[4][3] = dense.cdsqte*cs45/g[4];
132 col[3][4] = col[4][3]*g[4]/g[3]*BoltzFac[3][4];
133
134 double amat[5][5], bvec[5];
135 // homogeneous matrix - no source or sink - use conservation at 5th equation
136 for( long i=0; i<5; ++i )
137 {
138 amat[i][4] = 1.;
139 bvec[i] = 0.;
140 }
141 bvec[4] = abund;
142
143 /* level one balance equation */
144 amat[0][0] = col[0][1] + col[0][2] + col[0][3] + col[0][4] +pump01+pump02+pump03+pump04;
145 amat[1][0] = -a21 - col[1][0];
146 amat[2][0] = -a31 - col[2][0];
147 amat[3][0] = -a41 - col[3][0];
148 amat[4][0] = -a51 - col[4][0];
149
150 /* level two balance equation */
151 amat[0][1] = -col[0][1] - pump01;
152 amat[1][1] = col[1][0] + a21 + col[1][2] + col[1][3] + col[1][4];
153 amat[2][1] = -col[2][1] - a32;
154 amat[3][1] = -col[3][1] - a42;
155 amat[4][1] = -col[4][1] - a52;
156
157 /* level three balance equation */
158 amat[0][2] = -col[0][2] - pump02;
159 amat[1][2] = -col[1][2];
160 amat[2][2] = a31 + a32 + col[2][0] + col[2][1] + col[2][3] + col[2][4];
161 amat[3][2] = -col[3][2] - a43;
162 amat[4][2] = -col[4][2] - a53;
163
164 /* level four balance equation */
165 amat[0][3] = -col[0][3] - pump03;
166 amat[1][3] = -col[1][3];
167 amat[2][3] = -col[2][3];
168 amat[3][3] = a41 + col[3][0] + a42 + col[3][1] + a43 + col[3][2] + col[3][4];
169 amat[4][3] = -col[4][3] - a54;
170
171# if 0
172 // is it necessary to precondition the vars?
173 /* divide both sides of equation by largest number to stop overflow */
174 double dmax = -1e0;
175 for( i=0; i < 6; i++ )
176 for( j=0; j < 5; j++ )
177 dmax = MAX2(zz[i][j],dmax);
178
179 for( i=0; i < 6; i++ )
180 for( j=0; j < 5; j++ )
181 zz[i][j] /= dmax;
182# endif
183
184 int32 ipiv[5], ner=0;
185 /* solve matrix */
186 getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner);
187 getrs_wrapper('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner);
188
189 if( ner != 0 )
190 {
191 fprintf( ioQQQ, "DISASTER PROBLEM atom_pop5: dgetrs finds singular or ill-conditioned matrix\n" );
193 }
194
195 /* p(5) was very slightly negative (1e-40) for SII in dqher.in - highest
196 * level has smallest excitation rates may be closest to ill conditioned*/
197 p[1] = MAX2(0.e0,bvec[1]);
198 p[2] = MAX2(0.e0,bvec[2]);
199 p[3] = MAX2(0.e0,bvec[3]);
200 p[4] = MAX2(0.e0,bvec[4]);
201 // this should be majority and so numerically more
202 p[0] = abund - p[1] - p[2] - p[3] - p[4];
203
204 // level energies in ergs, needed for energy exchange
205 double Erg[5] , EnergyKelvin[5];
206 Erg[0] = 0.;
207 EnergyKelvin[0] = 0.;
208 for( long i=1; i<5; ++i )
209 {
210 Erg[i] = Erg[i-1] + EnerWN[i-1]*ERG1CM;
211 EnergyKelvin[i] = EnergyKelvin[i-1] + EnerWN[i-1]*T1CM;
212 }
213
214 *Cooling = 0.;
215 *CoolingDeriv = 0.;
216 for( long ihi=1; ihi<5; ++ihi )
217 {
218 for( long ilo=0; ilo<ihi; ++ilo )
219 {
220 double CoolOne = (p[ilo]*col[ilo][ihi] - p[ihi]*col[ihi][ilo])*
221 (Erg[ihi]-Erg[ilo]);
222
223 *Cooling += CoolOne;
224 *CoolingDeriv += CoolOne*( EnergyKelvin[ihi]*thermal.tsq1 - thermal.halfte );
225 }
226 }
227
228 return;
229}
t_abund abund
Definition abund.cpp:5
static double * amat
void atom_pop5(const double g[], const double EnerWN[], double cs12, double cs13, double cs14, double cs15, double cs23, double cs24, double cs25, double cs34, double cs35, double cs45, double a21, double a31, double a41, double a51, double a32, double a42, double a52, double a43, double a53, double a54, double p[], realnum abund, double *Cooling, double *CoolingDeriv, double pump01, double pump02, double pump03, double pump04)
Definition atom_pop5.cpp:13
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_dense dense
Definition dense.cpp:24
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double T1CM
Definition physconst.h:167
UNUSED const double ERG1CM
Definition physconst.h:164
static double * g
Definition species2.cpp:28
t_thermal thermal
Definition thermal.cpp:5
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)