00001
00002
00003
00004 #include "cddefines.h"
00005 #include "atoms.h"
00006 #include "trace.h"
00007 #include "iso.h"
00008 #include "dense.h"
00009 #include "mole.h"
00010 #include "opacity.h"
00011 #include "thermal.h"
00012 #include "gammas.h"
00013 #include "ionbal.h"
00014
00015 void IonMagne(void)
00016 {
00017 const int NDIM = ipMAGNESIUM+1;
00018
00019 static const double dicoef[2][NDIM] = {
00020 {4.49e-4,1.95e-3,5.12e-3,7.74e-3,1.17e-2,3.69e-2,3.63e-2,4.15e-2,8.86e-3,.252,.144,0.},
00021 {.021,.074,.323,.636,.807,.351,.548,.233,.318,.315,.291,0.}
00022 };
00023 static const double dite[2][NDIM] = {
00024 {5.01e4,6.06e5,4.69e5,3.74e5,3.28e5,4.80e5,3.88e5,3.39e5,2.11e5,1.40e7,1.50e7,0.},
00025 {2.81e4,1.44e6,7.55e5,7.88e5,1.02e6,9.73e5,7.38e5,3.82e5,1.54e6,2.64e6,3.09e6,0.}
00026 };
00027 static const double ditcrt[NDIM] = {4.0e3,7.4e4,6.6e4,5.5e4,4.4e4,4.5e4,
00028 4.5e4,5.0e5,3.4e4,2.4e6,4.0e6,1e20};
00029 static const double aa[NDIM] = {1.2044,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00030 static const double bb[NDIM] = {-4.6836,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00031 static const double cc[NDIM] = {7.6620,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00032 static const double dd[NDIM] = {-0.5930,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00033 static const double ff[NDIM] = {1.6260,0.1,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00034
00035 DEBUG_ENTRY( "IonMagne()" );
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 if( !dense.lgElmtOn[ipMAGNESIUM] )
00047 {
00048 atoms.xMg2Max = 0.;
00049 return;
00050 }
00051
00052 ion_zero(ipMAGNESIUM);
00053
00054 ion_photo(ipMAGNESIUM,false);
00055
00056
00057
00058
00059 ion_collis(ipMAGNESIUM);
00060
00061
00062 ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipMAGNESIUM);
00063
00064
00065 if( dense.IonLow[ipMAGNESIUM] <= 1 )
00066 {
00067 t_phoHeat dummy;
00068
00069 realnum rmg2l = (realnum)GammaK(opac.ipmgex,
00070 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0],opac.ipOpMgEx,1., &dummy );
00071 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] += rmg2l*atoms.popmg2;
00072
00073 if( nzone <= 1 )
00074 {
00075 atoms.xMg2Max = 0.;
00076 }
00077 else if( ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] > 1e-30 )
00078 {
00079
00080 atoms.xMg2Max = (realnum)(MAX2(atoms.xMg2Max,rmg2l*atoms.popmg2/
00081 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0]));
00082 }
00083 }
00084 else
00085 {
00086 atoms.xMg2Max = 0.;
00087 }
00088
00089
00090 if( dense.IonLow[ipMAGNESIUM] <= 0 )
00091 {
00092
00093
00094
00095 if( !co.lgNoCOMole )
00096 {
00097
00098 ionbal.PhotoRate_Shell[ipMAGNESIUM][0][3][0] +=
00099 CO_findrk("S+,Mg=>S,Mg+")*dense.xIonDense[ipSULPHUR][1] +
00100 CO_findrk("Si+,Mg=>Si,Mg+")*dense.xIonDense[ipSILICON][1] +
00101 CO_findrk("C+,Mg=>C,Mg+")*dense.xIonDense[ipCARBON][1];
00102
00103 }
00104 }
00105
00106
00107 ion_solver(ipMAGNESIUM,false);
00108
00109 if( trace.lgTrace && trace.lgHeavyBug )
00110 {
00111 fprintf( ioQQQ, " IonMagne returns; frac=" );
00112 for( int i=0; i < 10; i++ )
00113 {
00114 fprintf( ioQQQ, "%10.3e", dense.xIonDense[ipMAGNESIUM][i]/
00115 dense.gas_phase[ipMAGNESIUM] );
00116 }
00117 fprintf( ioQQQ, "\n" );
00118 }
00119 return;
00120 }