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 /*IonSilic determine ionization balance of Silicon */ 00004 #include "cddefines.h" 00005 #include "dense.h" 00006 #include "ionbal.h" 00007 00008 void IonSilic(void) 00009 { 00010 const int NDIM = ipSILICON+1; 00011 00012 static const double dicoef[2][NDIM] = { 00013 {1.10e-3,5.87e-3,5.03e-3,5.43e-3,8.86e-3,1.68e-2,2.49e-2,3.13e-2, 00014 4.25e-2,6.18e-2,1.38e-2,.327,.189,0.}, 00015 {0.,.753,.188,.450,0.,1.80,1.88,2.01,1.22,.303,1.42,.306,.286,0.} 00016 }; 00017 static const double dite[2][NDIM] = { 00018 {7.70e4,9.63e4,8.75e4,1.05e6,1.14e6,4.85e5,4.15e5,3.66e5,3.63e5, 00019 3.88e5,2.51e5,1.88e7,1.99e7,0.}, 00020 {0.,6.46e4,4.71e4,7.98e5,0.,1.03e6,1.91e6,2.11e6,2.14e6,1.12e6, 00021 3.93e6,3.60e6,4.14e6,0.} 00022 }; 00023 static const double ditcrt[NDIM] = {1.1e4,1.1e4,1.1e4,1.7e5,9.5e4,8.0e4, 00024 7.4e4,6.8e4,6.6e4,6.5e4,4.5e4,3.7e6,6.3e6,1e20}; 00025 static const double aa[NDIM] = {-0.0219,3.2163,0.1203,0.,0.,0.,0.,0., 00026 0.,0.,0.,0.,0.,0.}; 00027 static const double bb[NDIM] = {0.4364,-12.0571,-2.6900,0.,0.,0.,0., 00028 0.,0.,0.,0.,0.,0.,0.}; 00029 static const double cc[NDIM] = {0.0684,16.2118,19.1943,0.,0.,0.,0., 00030 0.,0.,0.,0.,0.,0.,0.}; 00031 static const double dd[NDIM] = {-0.0032,-0.5886,-0.1479,0.,0.,0.,0., 00032 0.,0.,0.,0.,0.,0.,0.}; 00033 static const double ff[NDIM] = {0.1342,0.5613,0.1118,0.1,0.,0.,0.,0., 00034 0.,0.,0.,0.,0.,0.}; 00035 00036 bool lgPrtDebug=false; 00037 00038 DEBUG_ENTRY( "IonSilic()" ); 00039 00040 /* silicon, atomic number 14 */ 00041 if( !dense.lgElmtOn[ipSILICON] ) 00042 return; 00043 00044 ion_zero(ipSILICON); 00045 00046 ion_photo(ipSILICON,false); 00047 00048 # if 0 00049 static long nzUsed = -1; 00050 static double OldRate = 0.; 00051 00052 /* 2008 apr 18, this appears to try to damp out changes in the valence 00053 * shell photo rate - perhaps due to Lya oscillations in dust free 00054 * environments? */ 00055 if( nzone > 1 && OldRate > 0. ) 00056 { 00057 if( nzone != nzUsed ) 00058 { 00059 ionbal.PhotoRate_Shell[ipSILICON][0][4][0] = 00060 (ionbal.PhotoRate_Shell[ipSILICON][0][4][0] + OldRate)/2.; 00061 nzUsed = nzone; 00062 } 00063 else 00064 { 00065 ionbal.PhotoRate_Shell[ipSILICON][0][4][0] = OldRate; 00066 } 00067 } 00068 OldRate = ionbal.PhotoRate_Shell[ipSILICON][0][4][0]; 00069 # endif 00070 00071 /* find collisional ionization rates */ 00072 ion_collis(ipSILICON); 00073 00074 /* get recombination coefficients */ 00075 /*lint -e64 type mismatch */ 00076 ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipSILICON); 00077 /*lint +e64 type mismatch */ 00078 00079 /* solve for ionization balance */ 00080 /*if( nzone>700 ) lgPrtDebug = true;*/ 00081 ion_solver(ipSILICON,lgPrtDebug); 00082 00083 # if 0 00084 /* 2008 apr 18, this is all dead code from when co network fed back from 00085 * the ionization */ 00086 if( trace.lgTrace && trace.lgHeavyBug ) 00087 { 00088 fprintf( ioQQQ, " silicon\t%.2f",fnzone ); 00089 for( int i=0; i <= 10; i++ ) 00090 { 00091 fprintf( ioQQQ, "\t%.3e", dense.xIonDense[ipSILICON][i]); 00092 } 00093 fprintf( ioQQQ, "\n" ); 00094 } 00095 fprintf( ioQQQ, "DEBUG silicon\t%.2f",fnzone ); 00096 for( i=0; i < 2; i++ ) 00097 { 00098 fprintf( ioQQQ, "\t%.3e", dense.xIonDense[ipSILICON][i]/dense.gas_phase[ipSILICON]); 00099 } 00100 fprintf( ioQQQ, "\t%.3e", co.hevmol[ipATSI]/SDIV(co.hevmol[ipSIP])); 00101 fprintf( ioQQQ, "\t%.3e", ionbal.RateIonizTot[ipSILICON][0]); 00102 fprintf( ioQQQ, "\t%.3e", dense.gas_phase[ipHYDROGEN] ); 00103 { 00104 # include "grainvar.h" 00105 fprintf( ioQQQ, "\t%.3e", gv.GrainChTrRate[ipSILICON][1][0] ); 00106 } 00107 fprintf( ioQQQ, "\n" ); 00108 # endif 00109 return; 00110 }