/*fivel do five level atom population and cooling */ #include #include "cddefines.h" #include "physconst.h" #include "edsqte.h" #include "phycon.h" #include "typmatrx.h" #include "linpack.h" #include "sexp.h" #include "matin1.h" #include "veclib.h" #include "fivel.h" void fivel(double g[], double ex[], 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[], float abund) { long int i, ipiv[5], j, job, ner; double c12, c13, c14, c15, c21, c23, c24, c25, c31, c32, c34, c35, c41, c42, c43, c45, c51, c52, c53, c54, e12, e13, e14, e15, e23, e24, e25, e34, e35, e45, tf; double amat[5][5], bvec[5], dmax, rcond, work[5], zz[6][6]; # ifdef DEBUG_FUN fputs( "<+>fivel()\n", debug_fp ); # endif /* following used for linpac matrix inversion */ /* tf = 1.438 / te */ tf = T1CM/phycon.te; /* quit if no species present */ if( abund <= 0. ) { p[0] = 0.; p[1] = 0.; p[2] = 0.; p[3] = 0.; p[4] = 0.; # ifdef DEBUG_FUN fputs( " <->fivel()\n", debug_fp ); # endif return; } /* get some Boltzmann factors */ e12 = sexp(ex[0]*tf); e23 = sexp(ex[1]*tf); e34 = sexp(ex[2]*tf); e45 = sexp(ex[3]*tf); e13 = e12*e23; e14 = e13*e34; e15 = e14*e45; e24 = e23*e34; e25 = e24*e45; e35 = e34*e45; /* quit it highest level Boltzmann factor too large */ if( e15 == 0. ) { p[0] = 0.; p[1] = 0.; p[2] = 0.; p[3] = 0.; p[4] = 0.; # ifdef DEBUG_FUN fputs( " <->fivel()\n", debug_fp ); # endif return; } /* get collision rates, * gcdsqte IS 8.629e-6 / sqrte * eden */ c21 = edsqte.cdsqte*cs12/g[1]; c12 = c21*g[1]/g[0]*e12; c31 = edsqte.cdsqte*cs13/g[2]; c13 = c31*g[2]/g[0]*e13; c41 = edsqte.cdsqte*cs14/g[3]; c14 = c41*g[3]/g[0]*e14; c51 = edsqte.cdsqte*cs15/g[4]; c15 = c51*g[4]/g[0]*e15; c32 = edsqte.cdsqte*cs23/g[2]; c23 = c32*g[2]/g[1]*e23; c42 = edsqte.cdsqte*cs24/g[3]; c24 = c42*g[3]/g[1]*e24; c52 = edsqte.cdsqte*cs25/g[4]; c25 = c52*g[4]/g[1]*e25; c43 = edsqte.cdsqte*cs34/g[3]; c34 = c43*g[3]/g[2]*e34; c53 = edsqte.cdsqte*cs35/g[4]; c35 = c53*g[4]/g[2]*e35; c54 = edsqte.cdsqte*cs45/g[4]; c45 = c54*g[4]/g[3]*e45; /* rate equations equal zero */ for( i=0; i < 5; i++ ) { zz[i][4] = 1.0; zz[5][i] = 0.; } zz[5][4] = abund; /* level one balance equation */ zz[0][0] = c12 + c13 + c14 + c15; zz[1][0] = -a21 - c21; zz[2][0] = -a31 - c31; zz[3][0] = -a41 - c41; zz[4][0] = -a51 - c51; /* level two balance equation */ zz[0][1] = -c12; zz[1][1] = c21 + a21 + c23 + c24 + c25; zz[2][1] = -c32 - a32; zz[3][1] = -c42 - a42; zz[4][1] = -c52 - a52; /* level three balance equation */ zz[0][2] = -c13; zz[1][2] = -c23; zz[2][2] = a31 + a32 + c31 + c32 + c34 + c35; zz[3][2] = -c43 - a43; zz[4][2] = -c53 - a53; /* level four balance equation */ zz[0][3] = -c14; zz[1][3] = -c24; zz[2][3] = -c34; zz[3][3] = a41 + c41 + a42 + c42 + a43 + c43 + c45; zz[4][3] = -c54 - a54; /* divide both sides of equation by largest number to stop overflow */ dmax = -1e0; for( i=0; i < 6; i++ ) { for( j=0; j < 5; j++ ) { dmax = MAX2(zz[i][j],dmax); } } for( i=0; i < 6; i++ ) { for( j=0; j < 5; j++ ) { zz[i][j] /= dmax; } } /* which matrix solver? */ if( strcmp(TypMatrx.chMatrix,"matin1 ") == 0 ) { /*matin1();*/ ner = 1; if( ner != 0 ) { fprintf( ioQQQ, " 5-level matrix error.\n" ); puts( "[Stop in fivel]" ); exit(1); } } else if( strcmp(TypMatrx.chMatrix,"linpack") == 0 ) { /* this one may be more robust */ for( j=0; j < 5; j++ ) { for( i=0; i < 5; i++ ) { amat[i][j] = zz[i][j]; } bvec[j] = zz[5][j]; } DGETRF(5,5,(double*)amat,5,ipiv,&ner); DGETRS('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner); if( ner != 0 ) { fprintf( ioQQQ, " 5-level dgetrs error\n" ); puts( "[Stop in fivel]" ); exit(1); } /* now put results back into z so rest of code treates only * one case - as if matin1 had been used */ for( i=0; i < 5; i++ ) { zz[5][i] = bvec[i]; } } else if( strcmp(TypMatrx.chMatrix,"veclib ") == 0 ) { /* Jason found this one on the Exemplar, distributed source just stops */ fprintf( ioQQQ, " this has not been checked since H atom conv\n" ); for( j=0; j < 5; j++ ) { for( i=0; i < 5; i++ ) { amat[i][j] = zz[i][j]; } bvec[j] = zz[5][j]; } job = 0; rcond = 0.; dgeco((double*)amat,5,5,ipiv,rcond,work); dgesl((double*)amat,5,5,ipiv,bvec,job); /* now put results back into z so rest of code treates only * one case - as if matin1 had been used */ for( i=0; i < 5; i++ ) { zz[5][i] = bvec[i]; } puts( "[Stop in fivel]" ); exit(1); } else { fprintf( ioQQQ, " chMatrix type insane in hydrogen, was%7.7s\n", TypMatrx.chMatrix ); puts( "[Stop in fivel]" ); exit(1); } /* 666 error! p(5) was very slightly negative (1e-40) for SII in dqher.in */ p[1] = MAX2(0.e0,zz[5][1]); p[2] = MAX2(0.e0,zz[5][2]); p[3] = MAX2(0.e0,zz[5][3]); p[4] = MAX2(0.e0,zz[5][4]); p[0] = abund - p[1] - p[2] - p[3] - p[4]; # ifdef DEBUG_FUN fputs( " <->fivel()\n", debug_fp ); # endif return; }