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 float 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
00101 DEBUG_EXIT( "atom_pop5()" );
00102 return;
00103 }
00104
00105
00106 e12 = sexp(ex[0]*tf);
00107 e23 = sexp(ex[1]*tf);
00108 e34 = sexp(ex[2]*tf);
00109 e45 = sexp(ex[3]*tf);
00110 e13 = e12*e23;
00111 e14 = e13*e34;
00112 e15 = e14*e45;
00113 e24 = e23*e34;
00114 e25 = e24*e45;
00115 e35 = e34*e45;
00116
00117
00118 if( e15 == 0. )
00119 {
00120 p[0] = 0.;
00121 p[1] = 0.;
00122 p[2] = 0.;
00123 p[3] = 0.;
00124 p[4] = 0.;
00125
00126 DEBUG_EXIT( "atom_pop5()" );
00127 return;
00128 }
00129
00130
00131
00132 c21 = dense.cdsqte*cs12/g[1];
00133 c12 = c21*g[1]/g[0]*e12;
00134
00135 c31 = dense.cdsqte*cs13/g[2];
00136 c13 = c31*g[2]/g[0]*e13;
00137
00138 c41 = dense.cdsqte*cs14/g[3];
00139 c14 = c41*g[3]/g[0]*e14;
00140
00141 c51 = dense.cdsqte*cs15/g[4];
00142 c15 = c51*g[4]/g[0]*e15;
00143
00144 c32 = dense.cdsqte*cs23/g[2];
00145 c23 = c32*g[2]/g[1]*e23;
00146
00147 c42 = dense.cdsqte*cs24/g[3];
00148 c24 = c42*g[3]/g[1]*e24;
00149
00150 c52 = dense.cdsqte*cs25/g[4];
00151 c25 = c52*g[4]/g[1]*e25;
00152
00153 c43 = dense.cdsqte*cs34/g[3];
00154 c34 = c43*g[3]/g[2]*e34;
00155
00156 c53 = dense.cdsqte*cs35/g[4];
00157 c35 = c53*g[4]/g[2]*e35;
00158
00159 c54 = dense.cdsqte*cs45/g[4];
00160 c45 = c54*g[4]/g[3]*e45;
00161
00162
00163 for( i=0; i < 5; i++ )
00164 {
00165 zz[i][4] = 1.0;
00166 zz[5][i] = 0.;
00167 }
00168 zz[5][4] = abund;
00169
00170
00171 zz[0][0] = c12 + c13 + c14 + c15;
00172 zz[1][0] = -a21 - c21;
00173 zz[2][0] = -a31 - c31;
00174 zz[3][0] = -a41 - c41;
00175 zz[4][0] = -a51 - c51;
00176
00177
00178 zz[0][1] = -c12;
00179 zz[1][1] = c21 + a21 + c23 + c24 + c25;
00180 zz[2][1] = -c32 - a32;
00181 zz[3][1] = -c42 - a42;
00182 zz[4][1] = -c52 - a52;
00183
00184
00185 zz[0][2] = -c13;
00186 zz[1][2] = -c23;
00187 zz[2][2] = a31 + a32 + c31 + c32 + c34 + c35;
00188 zz[3][2] = -c43 - a43;
00189 zz[4][2] = -c53 - a53;
00190
00191
00192 zz[0][3] = -c14;
00193 zz[1][3] = -c24;
00194 zz[2][3] = -c34;
00195 zz[3][3] = a41 + c41 + a42 + c42 + a43 + c43 + c45;
00196 zz[4][3] = -c54 - a54;
00197
00198
00199 dmax = -1e0;
00200
00201 for( i=0; i < 6; i++ )
00202 {
00203 for( j=0; j < 5; j++ )
00204 {
00205 dmax = MAX2(zz[i][j],dmax);
00206 }
00207 }
00208
00209 for( i=0; i < 6; i++ )
00210 {
00211 for( j=0; j < 5; j++ )
00212 {
00213 zz[i][j] /= dmax;
00214 }
00215 }
00216
00217
00218 for( j=0; j < 5; j++ )
00219 {
00220 for( i=0; i < 5; i++ )
00221 {
00222 amat[i][j] = zz[i][j];
00223 }
00224 bvec[j] = zz[5][j];
00225 }
00226
00227 ner = 0;
00228
00229
00230
00231 getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner);
00232 getrs_wrapper('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner);
00233
00234
00235 if( ner != 0 )
00236 {
00237 fprintf( ioQQQ, " atom_pop5: dgetrs finds singular or ill-conditioned matrix\n" );
00238 puts( "[Stop in atom_pop5]" );
00239 cdEXIT(EXIT_FAILURE);
00240 }
00241
00242
00243
00244 for( i=0; i < 5; i++ )
00245 {
00246 zz[5][i] = bvec[i];
00247 }
00248
00250 p[1] = MAX2(0.e0,zz[5][1]);
00251 p[2] = MAX2(0.e0,zz[5][2]);
00252 p[3] = MAX2(0.e0,zz[5][3]);
00253 p[4] = MAX2(0.e0,zz[5][4]);
00254 p[0] = abund - p[1] - p[2] - p[3] - p[4];
00255
00256 DEBUG_EXIT( "atom_pop5()" );
00257 return;
00258 }
00259