00001
00002
00003
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
00012 void atom_pop5(
00013
00014 double g[],
00015
00016
00017
00018 double ex[],
00019
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
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
00042 double p[],
00043
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
00090 tf = T1CM/phycon.te;
00091
00092
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
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
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
00127
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
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
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
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
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
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
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
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
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
00236
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 }