cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atom_pop5.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 "dense.h"
8 #include "thirdparty.h"
9 #include "atoms.h"
10 
11 /*atom_pop5 do five level atom population and cooling */
12 void atom_pop5(
13  /* vector giving statistical weights on the five levels */
14  double g[],
15  /* vector giving the excitation energy differences of the 5 levels. The energies
16  * are the energy in wavenumbers between adjacent levels. So ex[0] is the energy
17  * 1-2, ex[1] is the energy between 2-3, etc */
18  double ex[],
19  /* the collision strengths for the levels */
20  double cs12,
21  double cs13,
22  double cs14,
23  double cs15,
24  double cs23,
25  double cs24,
26  double cs25,
27  double cs34,
28  double cs35,
29  double cs45,
30  /* the transition probabilities between the various levels */
31  double a21,
32  double a31,
33  double a41,
34  double a51,
35  double a32,
36  double a42,
37  double a52,
38  double a43,
39  double a53,
40  double a54,
41  /* the destroyed level populations (units cm-3) for the five levels */
42  double p[],
43  /* the total density of this ion */
44  realnum abund)
45 {
46  int32 ipiv[5], ner;
47 
48  long int i, j;
49 
50  double c12,
51  c13,
52  c14,
53  c15,
54  c21,
55  c23,
56  c24,
57  c25,
58  c31,
59  c32,
60  c34,
61  c35,
62  c41,
63  c42,
64  c43,
65  c45,
66  c51,
67  c52,
68  c53,
69  c54,
70  e12,
71  e13,
72  e14,
73  e15,
74  e23,
75  e24,
76  e25,
77  e34,
78  e35,
79  e45,
80  tf;
81 
82  double amat[5][5],
83  bvec[5],
84  dmax,
85  zz[6][6];
86 
87  DEBUG_ENTRY( "atom_pop5()" );
88 
89  /* tf = 1.438 / te */
90  tf = T1CM/phycon.te;
91 
92  /* quit if no species present */
93  if( abund <= 0. )
94  {
95  p[0] = 0.;
96  p[1] = 0.;
97  p[2] = 0.;
98  p[3] = 0.;
99  p[4] = 0.;
100  return;
101  }
102 
103  /* get some Boltzmann factors */
104  e12 = sexp(ex[0]*tf);
105  e23 = sexp(ex[1]*tf);
106  e34 = sexp(ex[2]*tf);
107  e45 = sexp(ex[3]*tf);
108  e13 = e12*e23;
109  e14 = e13*e34;
110  e15 = e14*e45;
111  e24 = e23*e34;
112  e25 = e24*e45;
113  e35 = e34*e45;
114 
115  /* quit it highest level Boltzmann factor too large */
116  if( e15 == 0. )
117  {
118  p[0] = 0.;
119  p[1] = 0.;
120  p[2] = 0.;
121  p[3] = 0.;
122  p[4] = 0.;
123  return;
124  }
125 
126  /* get collision rates,
127  * dense.cdsqte is 8.629e-6 / sqrte * eden */
128  c21 = dense.cdsqte*cs12/g[1];
129  c12 = c21*g[1]/g[0]*e12;
130 
131  c31 = dense.cdsqte*cs13/g[2];
132  c13 = c31*g[2]/g[0]*e13;
133 
134  c41 = dense.cdsqte*cs14/g[3];
135  c14 = c41*g[3]/g[0]*e14;
136 
137  c51 = dense.cdsqte*cs15/g[4];
138  c15 = c51*g[4]/g[0]*e15;
139 
140  c32 = dense.cdsqte*cs23/g[2];
141  c23 = c32*g[2]/g[1]*e23;
142 
143  c42 = dense.cdsqte*cs24/g[3];
144  c24 = c42*g[3]/g[1]*e24;
145 
146  c52 = dense.cdsqte*cs25/g[4];
147  c25 = c52*g[4]/g[1]*e25;
148 
149  c43 = dense.cdsqte*cs34/g[3];
150  c34 = c43*g[3]/g[2]*e34;
151 
152  c53 = dense.cdsqte*cs35/g[4];
153  c35 = c53*g[4]/g[2]*e35;
154 
155  c54 = dense.cdsqte*cs45/g[4];
156  c45 = c54*g[4]/g[3]*e45;
157 
158  /* rate equations equal zero */
159  for( i=0; i < 5; i++ )
160  {
161  zz[i][4] = 1.0;
162  zz[5][i] = 0.;
163  }
164  zz[5][4] = abund;
165 
166  /* level one balance equation */
167  zz[0][0] = c12 + c13 + c14 + c15;
168  zz[1][0] = -a21 - c21;
169  zz[2][0] = -a31 - c31;
170  zz[3][0] = -a41 - c41;
171  zz[4][0] = -a51 - c51;
172 
173  /* level two balance equation */
174  zz[0][1] = -c12;
175  zz[1][1] = c21 + a21 + c23 + c24 + c25;
176  zz[2][1] = -c32 - a32;
177  zz[3][1] = -c42 - a42;
178  zz[4][1] = -c52 - a52;
179 
180  /* level three balance equation */
181  zz[0][2] = -c13;
182  zz[1][2] = -c23;
183  zz[2][2] = a31 + a32 + c31 + c32 + c34 + c35;
184  zz[3][2] = -c43 - a43;
185  zz[4][2] = -c53 - a53;
186 
187  /* level four balance equation */
188  zz[0][3] = -c14;
189  zz[1][3] = -c24;
190  zz[2][3] = -c34;
191  zz[3][3] = a41 + c41 + a42 + c42 + a43 + c43 + c45;
192  zz[4][3] = -c54 - a54;
193 
194  /* divide both sides of equation by largest number to stop overflow */
195  dmax = -1e0;
196 
197  for( i=0; i < 6; i++ )
198  {
199  for( j=0; j < 5; j++ )
200  {
201  dmax = MAX2(zz[i][j],dmax);
202  }
203  }
204 
205  for( i=0; i < 6; i++ )
206  {
207  for( j=0; j < 5; j++ )
208  {
209  zz[i][j] /= dmax;
210  }
211  }
212 
213  /* this one may be more robust */
214  for( j=0; j < 5; j++ )
215  {
216  for( i=0; i < 5; i++ )
217  {
218  amat[i][j] = zz[i][j];
219  }
220  bvec[j] = zz[5][j];
221  }
222 
223  ner = 0;
224 
225  /* solve matrix */
226  getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner);
227  getrs_wrapper('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner);
228 
229  if( ner != 0 )
230  {
231  fprintf( ioQQQ, " atom_pop5: dgetrs finds singular or ill-conditioned matrix\n" );
232  cdEXIT(EXIT_FAILURE);
233  }
234 
235  /* now put results back into z so rest of code treats only
236  * one case - as if matin1 had been used */
237  for( i=0; i < 5; i++ )
238  {
239  zz[5][i] = bvec[i];
240  }
241 
243  p[1] = MAX2(0.e0,zz[5][1]);
244  p[2] = MAX2(0.e0,zz[5][2]);
245  p[3] = MAX2(0.e0,zz[5][3]);
246  p[4] = MAX2(0.e0,zz[5][4]);
247  p[0] = abund - p[1] - p[2] - p[3] - p[4];
248  return;
249 }

Generated for cloudy by doxygen 1.8.3.1