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 /*IonNitro ionization balance for nitrogen */ 00004 #include "cddefines.h" 00005 #include "dense.h" 00006 #include "thermal.h" 00007 #include "atoms.h" 00008 #include "opacity.h" 00009 #include "trace.h" 00010 #include "conv.h" 00011 #include "gammas.h" 00012 #include "ionbal.h" 00013 00014 void IonNitro(void) 00015 { 00016 const int NDIM = ipNITROGEN+1; 00017 00018 static const double dicoef[2][NDIM] = { 00019 {2.98e-3,7.41e-3,1.13e-2,2.62e-3,7.5e-2,4.61e-2,0.}, 00020 {0.,.0764,.164,.243,.35,.309,0.} 00021 }; 00022 static const double dite[2][NDIM] = { 00023 {2.2e5,2.01e5,1.72e5,1.02e5,4.75e6,5.44e6,0.}, 00024 {0.,7.37e4,2.25e5,1.25e5,8.35e5,1.14e6,0.} 00025 }; 00026 static const double ditcrt[NDIM] = {1.8e4,1.8e4,2.4e4,1.5e4,6.8e5,1.0e6,1e20}; 00027 static const double aa[NDIM] = {0.0,0.0320,-0.8806,0.4134,0.,0.,0.}; 00028 static const double bb[NDIM] = {0.6310,-0.6624,11.2406,-4.6319,0.,0.,0.}; 00029 static const double cc[NDIM] = {0.1990,4.3191,30.7066,25.9172,0.,0.,0.}; 00030 static const double dd[NDIM] = {-0.0197,0.0003,-1.1721,-2.2290,0.,0.,0.}; 00031 static const double ff[NDIM] = {0.4398,0.5946,0.6127,0.2360,0.1,0.,0.}; 00032 00033 bool lgDebug; 00034 double save_rec; 00035 00036 DEBUG_ENTRY( "IonNitro()" ); 00037 00038 /* nitrogen, atomic number 7 */ 00039 if( !dense.lgElmtOn[ipNITROGEN] ) 00040 { 00041 return; 00042 } 00043 00044 ion_zero(ipNITROGEN); 00045 00046 ion_photo(ipNITROGEN,false); 00047 00048 /* find collisional ionization rates */ 00049 ion_collis(ipNITROGEN); 00050 00051 /* get recombination coefficients */ 00052 /*lint -e64 type mismatch */ 00053 ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipNITROGEN); 00054 /*lint +e64 type mismatch */ 00055 00056 /* >>chng 04 jul 09, synch up ion and co solvers. 00057 * the co solver will have a different N/N+ balance that 00058 * is strongly affected by the chemistry if we are in neutral gas 00059 * in this case use co solver's N/N+ balance */ 00060 save_rec = ionbal.RateRecomTot[ipNITROGEN][0]; 00061 /*>>chng 04 apr 24, do not test for ionhigh being 1, 00062 * no reason for test on upper stage of ionization, 00063 * if codrive called but not evaluated then hevmol is all zero */ 00064 /* >>chng 04 sep 10, rm check on search phase, no reason for it */ 00065 # if 0 00066 if( dense.IonLow[ipNITROGEN]==0 && 00067 co.hevmol[ipNP] > SMALLFLOAT && 00068 ionbal.RateIonizTot[ipNITROGEN][0]*co.hevmol[ipATN]> 0. ) 00069 { 00070 ionbal.RateRecomTot[ipNITROGEN][0] = 00071 ionbal.RateIonizTot[ipNITROGEN][0]* 00072 co.hevmol[ipATN]/co.hevmol[ipNP]; 00073 } 00074 # endif 00075 00076 /* photoionization from 2D of NI, is atomic nitrogen present? */ 00077 if( dense.xIonDense[ipNITROGEN][0] > 0. ) 00078 { 00079 atoms.d5200r = (realnum)GammaK(opac.in1[0],opac.in1[1],opac.in1[2],1.); 00080 /*aeff = atoms.d5200r + 1.28e-5;*/ 00081 /* >>chng 01 sep 02, use values from coolnitr */ 00082 /*cs5200 = 1.32e-4*phycon.te/(phycon.te10*phycon.te01);*/ 00083 /* uses 1 for total pop of atom - so is normalized to unity */ 00084 /*atoms.p2nit = (realnum)(atom_pop2(cs5200,4.,10.,aeff,2.769e4,1.)/aeff);*/ 00085 /* if very first call then set to zero */ 00086 if( !conv.nTotalIoniz ) 00087 atoms.p2nit = 0.; 00088 /* valence shell photoionization, [0][6] => atomic nitrogen */ 00089 ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0]* 00090 (1. - atoms.p2nit) + atoms.p2nit*atoms.d5200r; 00091 ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1]* 00092 (1. - atoms.p2nit) + thermal.HeatNet*atoms.p2nit; 00093 } 00094 else 00095 { 00096 atoms.p2nit = 0.; 00097 atoms.d5200r = 0.; 00098 } 00099 00100 /* solve for ionization balance */ 00101 # if 0 00102 broken();/* rm following true */ 00103 if(nzone == 1 ) 00104 lgDebug = true; 00105 else 00106 lgDebug = false; 00107 # endif 00108 lgDebug = false; 00109 ion_solver(ipNITROGEN,lgDebug); 00110 00111 /* reset the var we just hosed */ 00112 if( save_rec > 0.) 00113 ionbal.RateRecomTot[ipNITROGEN][0] = save_rec; 00114 00115 if( trace.lgTrace && trace.lgHeavyBug ) 00116 { 00117 fprintf( ioQQQ, " IonNitro retun; frac=" ); 00118 for( int i=0; i < 8; i++ ) 00119 { 00120 fprintf( ioQQQ, "%10.3e", dense.xIonDense[ipNITROGEN][i]/ 00121 dense.gas_phase[ipNITROGEN] ); 00122 } 00123 fprintf( ioQQQ, "\n" ); 00124 } 00125 return; 00126 }