cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*IonMagne ionization balance for magnesium */ 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 realnum rmg2l; 00035 00036 DEBUG_ENTRY( "IonMagne()" ); 00037 00038 /* magnesium nelem=12 00039 * 00040 * rates from Shull and van Steenberg, Ap.J. Sup 48, 95. */ 00041 00042 /* rec from +9, +10, +11 from Arnauld et al '85 */ 00043 /* Pequignot and Aldrovandi Ast Ap 161, 169. */ 00044 00045 /* mg i from nuss+storey, who say that =>Mgii is very small */ 00046 00047 if( !dense.lgElmtOn[ipMAGNESIUM] ) 00048 { 00049 atoms.xMg2Max = 0.; 00050 return; 00051 } 00052 00053 ion_zero(ipMAGNESIUM); 00054 00055 ion_photo(ipMAGNESIUM,false); 00056 /* debugging printout for shell photo rates - 0 for atom, last true, also print details */ 00057 /*GammaPrtRate( ioQQQ , 0 , ipMAGNESIUM , true );*/ 00058 00059 /* find collisional ionization rates */ 00060 ion_collis(ipMAGNESIUM); 00061 00062 /* get recombination coefficients */ 00063 ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipMAGNESIUM); 00064 00065 /* is the atom present*/ 00066 if( dense.IonLow[ipMAGNESIUM] <= 0 ) 00067 { 00068 /* can only do this one time 00069 * photoionization from excited upper state of 2798 */ 00070 rmg2l = (realnum)GammaK(opac.ipmgex,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0],opac.ipOpMgEx,1.); 00071 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] += rmg2l*atoms.popmg2; 00072 00073 /* >>chng 06 Feb 28 -- NPA. Add in charge transfer ionization of Fe with S+, Si+, and C+ */ 00074 /* only include this if molecular network is enabled - otherwise no feedback onto 00075 * Si+, S+, and C+ soln */ 00076 if( !co.lgNoCOMole ) 00077 { 00078 /* Use sink rate from last completed state of network, including _all_ present & future Mg sinks */ 00079 ionbal.PhotoRate_Shell[ipMAGNESIUM][0][3][0] += 00080 CO_findrk("S+,Mg=>S,Mg+")*dense.xIonDense[ipSULPHUR][1] + 00081 CO_findrk("Si+,Mg=>Si,Mg+")*dense.xIonDense[ipSILICON][1] + 00082 CO_findrk("C+,Mg=>C,Mg+")*dense.xIonDense[ipCARBON][1]; 00083 /* CO_sink_rate("Mg"); */ 00084 } 00085 00086 if( nzone <= 1 ) 00087 { 00088 atoms.xMg2Max = 0.; 00089 } 00090 else if( ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] > 1e-30 ) 00091 { 00092 /* remember max rel photo rate */ 00093 atoms.xMg2Max = (realnum)(MAX2(atoms.xMg2Max,rmg2l*atoms.popmg2/ 00094 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0])); 00095 } 00096 } 00097 else 00098 { 00099 atoms.xMg2Max = 0.; 00100 } 00101 00102 /* solve for ionization balance */ 00103 ion_solver(ipMAGNESIUM,false); 00104 00105 if( trace.lgTrace && trace.lgHeavyBug ) 00106 { 00107 fprintf( ioQQQ, " IonMagne returns; frac=" ); 00108 for( int i=0; i < 10; i++ ) 00109 { 00110 fprintf( ioQQQ, "%10.3e", dense.xIonDense[ipMAGNESIUM][i]/ 00111 dense.gas_phase[ipMAGNESIUM] ); 00112 } 00113 fprintf( ioQQQ, "\n" ); 00114 } 00115 return; 00116 }