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