cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*atom_pop5 do five level atom population and cooling */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "phycon.h" 00007 #include "dense.h" 00008 #include "thirdparty.h" 00009 #include "atoms.h" 00010 00011 /*atom_pop5 do five level atom population and cooling */ 00012 void atom_pop5( 00013 /* vector giving statistical weights on the five levels */ 00014 double g[], 00015 /* vector giving the excitation energy differences of the 5 levels. The energies 00016 * are the energy in wavenumbers between adjacent levels. So ex[0] is the energy 00017 * 1-2, ex[1] is the energy between 2-3, etc */ 00018 double ex[], 00019 /* the collision strengths for the levels */ 00020 double cs12, 00021 double cs13, 00022 double cs14, 00023 double cs15, 00024 double cs23, 00025 double cs24, 00026 double cs25, 00027 double cs34, 00028 double cs35, 00029 double cs45, 00030 /* the transition probabilities between the various levels */ 00031 double a21, 00032 double a31, 00033 double a41, 00034 double a51, 00035 double a32, 00036 double a42, 00037 double a52, 00038 double a43, 00039 double a53, 00040 double a54, 00041 /* the destroyed level populations (units cm-3) for the five levels */ 00042 double p[], 00043 /* the total density of this ion */ 00044 realnum abund) 00045 { 00046 int32 ipiv[5], ner; 00047 00048 long int i, j; 00049 00050 double c12, 00051 c13, 00052 c14, 00053 c15, 00054 c21, 00055 c23, 00056 c24, 00057 c25, 00058 c31, 00059 c32, 00060 c34, 00061 c35, 00062 c41, 00063 c42, 00064 c43, 00065 c45, 00066 c51, 00067 c52, 00068 c53, 00069 c54, 00070 e12, 00071 e13, 00072 e14, 00073 e15, 00074 e23, 00075 e24, 00076 e25, 00077 e34, 00078 e35, 00079 e45, 00080 tf; 00081 00082 double amat[5][5], 00083 bvec[5], 00084 dmax, 00085 zz[6][6]; 00086 00087 DEBUG_ENTRY( "atom_pop5()" ); 00088 00089 /* tf = 1.438 / te */ 00090 tf = T1CM/phycon.te; 00091 00092 /* quit if no species present */ 00093 if( abund <= 0. ) 00094 { 00095 p[0] = 0.; 00096 p[1] = 0.; 00097 p[2] = 0.; 00098 p[3] = 0.; 00099 p[4] = 0.; 00100 return; 00101 } 00102 00103 /* get some Boltzmann factors */ 00104 e12 = sexp(ex[0]*tf); 00105 e23 = sexp(ex[1]*tf); 00106 e34 = sexp(ex[2]*tf); 00107 e45 = sexp(ex[3]*tf); 00108 e13 = e12*e23; 00109 e14 = e13*e34; 00110 e15 = e14*e45; 00111 e24 = e23*e34; 00112 e25 = e24*e45; 00113 e35 = e34*e45; 00114 00115 /* quit it highest level Boltzmann factor too large */ 00116 if( e15 == 0. ) 00117 { 00118 p[0] = 0.; 00119 p[1] = 0.; 00120 p[2] = 0.; 00121 p[3] = 0.; 00122 p[4] = 0.; 00123 return; 00124 } 00125 00126 /* get collision rates, 00127 * dense.cdsqte is 8.629e-6 / sqrte * eden */ 00128 c21 = dense.cdsqte*cs12/g[1]; 00129 c12 = c21*g[1]/g[0]*e12; 00130 00131 c31 = dense.cdsqte*cs13/g[2]; 00132 c13 = c31*g[2]/g[0]*e13; 00133 00134 c41 = dense.cdsqte*cs14/g[3]; 00135 c14 = c41*g[3]/g[0]*e14; 00136 00137 c51 = dense.cdsqte*cs15/g[4]; 00138 c15 = c51*g[4]/g[0]*e15; 00139 00140 c32 = dense.cdsqte*cs23/g[2]; 00141 c23 = c32*g[2]/g[1]*e23; 00142 00143 c42 = dense.cdsqte*cs24/g[3]; 00144 c24 = c42*g[3]/g[1]*e24; 00145 00146 c52 = dense.cdsqte*cs25/g[4]; 00147 c25 = c52*g[4]/g[1]*e25; 00148 00149 c43 = dense.cdsqte*cs34/g[3]; 00150 c34 = c43*g[3]/g[2]*e34; 00151 00152 c53 = dense.cdsqte*cs35/g[4]; 00153 c35 = c53*g[4]/g[2]*e35; 00154 00155 c54 = dense.cdsqte*cs45/g[4]; 00156 c45 = c54*g[4]/g[3]*e45; 00157 00158 /* rate equations equal zero */ 00159 for( i=0; i < 5; i++ ) 00160 { 00161 zz[i][4] = 1.0; 00162 zz[5][i] = 0.; 00163 } 00164 zz[5][4] = abund; 00165 00166 /* level one balance equation */ 00167 zz[0][0] = c12 + c13 + c14 + c15; 00168 zz[1][0] = -a21 - c21; 00169 zz[2][0] = -a31 - c31; 00170 zz[3][0] = -a41 - c41; 00171 zz[4][0] = -a51 - c51; 00172 00173 /* level two balance equation */ 00174 zz[0][1] = -c12; 00175 zz[1][1] = c21 + a21 + c23 + c24 + c25; 00176 zz[2][1] = -c32 - a32; 00177 zz[3][1] = -c42 - a42; 00178 zz[4][1] = -c52 - a52; 00179 00180 /* level three balance equation */ 00181 zz[0][2] = -c13; 00182 zz[1][2] = -c23; 00183 zz[2][2] = a31 + a32 + c31 + c32 + c34 + c35; 00184 zz[3][2] = -c43 - a43; 00185 zz[4][2] = -c53 - a53; 00186 00187 /* level four balance equation */ 00188 zz[0][3] = -c14; 00189 zz[1][3] = -c24; 00190 zz[2][3] = -c34; 00191 zz[3][3] = a41 + c41 + a42 + c42 + a43 + c43 + c45; 00192 zz[4][3] = -c54 - a54; 00193 00194 /* divide both sides of equation by largest number to stop overflow */ 00195 dmax = -1e0; 00196 00197 for( i=0; i < 6; i++ ) 00198 { 00199 for( j=0; j < 5; j++ ) 00200 { 00201 dmax = MAX2(zz[i][j],dmax); 00202 } 00203 } 00204 00205 for( i=0; i < 6; i++ ) 00206 { 00207 for( j=0; j < 5; j++ ) 00208 { 00209 zz[i][j] /= dmax; 00210 } 00211 } 00212 00213 /* this one may be more robust */ 00214 for( j=0; j < 5; j++ ) 00215 { 00216 for( i=0; i < 5; i++ ) 00217 { 00218 amat[i][j] = zz[i][j]; 00219 } 00220 bvec[j] = zz[5][j]; 00221 } 00222 00223 ner = 0; 00224 00225 /* solve matrix */ 00226 getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner); 00227 getrs_wrapper('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner); 00228 00229 if( ner != 0 ) 00230 { 00231 fprintf( ioQQQ, " atom_pop5: dgetrs finds singular or ill-conditioned matrix\n" ); 00232 cdEXIT(EXIT_FAILURE); 00233 } 00234 00235 /* now put results back into z so rest of code treats only 00236 * one case - as if matin1 had been used */ 00237 for( i=0; i < 5; i++ ) 00238 { 00239 zz[5][i] = bvec[i]; 00240 } 00241 00243 p[1] = MAX2(0.e0,zz[5][1]); 00244 p[2] = MAX2(0.e0,zz[5][2]); 00245 p[3] = MAX2(0.e0,zz[5][3]); 00246 p[4] = MAX2(0.e0,zz[5][4]); 00247 p[0] = abund - p[1] - p[2] - p[3] - p[4]; 00248 return; 00249 }