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 /*IonCarbo compute ionization balance for carbon */ 00004 #include "cddefines.h" 00005 #include "thermal.h" 00006 #include "carb.h" 00007 #include "trace.h" 00008 #include "dense.h" 00009 #include "phycon.h" 00010 #include "hmi.h" 00011 #include "mole.h" 00012 #include "rfield.h" 00013 #include "punch.h" 00014 #include "ionbal.h" 00015 00016 void IonCarbo(void) 00017 { 00018 const int NDIM = ipCARBON+1; 00019 00020 static const double dicoef[2][NDIM] = { 00021 {2.54e-3,6.15e-3,1.62e-3,4.78e-2,3.22e-2,0.}, {4.42e-2,5.88e-2,0.343,0.362,0.315,0.} 00022 }; 00023 static const double dite[2][NDIM] = { 00024 {1.57e5,1.41e5,8.19e4,3.44e6,4.06e6,0.}, {3.74e5,1.41e5,1.59e5,5.87e5,8.31e5,0.} 00025 }; 00026 static const double ditcrt[NDIM] = {1.2e4,1.2e4,1.1e4,4.4e5,7.0e5,1e20}; 00027 /* >>refer C DR Nussbaumer H. & Storey, P 1983, A&A, 126, 75 */ 00028 static const double aa[NDIM] = {.0108,1.8267,2.3196,0.,0.,0.}; 00029 static const double bb[NDIM] = {-0.1075,4.1012,10.7328,0.,0.,0.}; 00030 static const double cc[NDIM] = {.2810,4.8443,6.8830,0.,0.,0.}; 00031 static const double dd[NDIM] = {-0.0193,.2261,-0.1824,0.,0.,0.}; 00032 /* NB C+ - C0 recom numerically unstable below 2000K, use est instead */ 00033 static const double ff[NDIM] = {0.000,0.5960,0.4101,0.1,0.1,0.}; 00034 00035 double save_rec; 00036 00037 DEBUG_ENTRY( "IonCarbo()" ); 00038 00039 if( trace.lgTrace && trace.lgCarBug ) 00040 { 00041 fprintf( ioQQQ, " IonCarbo called.\n" ); 00042 } 00043 00044 if( !dense.lgElmtOn[ipCARBON] ) 00045 { 00046 carb.p1909 = 0.; 00047 carb.p2326 = 0.; 00048 thermal.heating[ipCARBON][9] = 0.; 00049 return; 00050 } 00051 00052 /* zero out ionization balance arrays */ 00053 ion_zero(ipCARBON); 00054 00055 ion_photo(ipCARBON,false); 00056 00057 /* >>chng 05 aug 10, add Leiden hack here to get same C0 photo rate 00058 * as UMIST - negates difference in grain opacities */ 00059 if(!co.lgUMISTrates) 00060 { 00061 int nelem=ipCARBON , ion=0 , ns=2; 00062 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 00063 (HMRATE((1e-10)*3.0,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face* 00064 exp(-(3.0*rfield.extin_mag_V_point))/1.66)); 00065 /* heating rates */ 00066 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.; 00067 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.; 00068 /*fprintf(ioQQQ,"DEBUG hack %li\tav\t%.2e\tgo\t%.2e\trate\t%.2e\tabun\t%.2e\n", 00069 nzone , 00070 rfield.extin_mag_V_point, 00071 hmi.UV_Cont_rel2_Habing_TH85_face, 00072 ionbal.PhotoRate_Shell[nelem][ion][ns][0] , 00073 dense.xIonDense[ipCARBON][1]/SDIV(dense.xIonDense[ipCARBON][0]) );*/ 00074 } 00075 00076 /* find collisional ionization rates */ 00077 ion_collis(ipCARBON); 00078 00079 /* get recombination coefficients */ 00080 /*lint -e64 double = pointer */ 00081 ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipCARBON); 00082 /*lint +e64 double = pointer */ 00083 00084 /* Photoproduction of 3P */ 00085 carb.p1909 = ionbal.PhotoRate_Shell[ipCARBON][1][1][0]; 00086 00087 /* photoproduction of C II] 2326 */ 00088 carb.p2326 = ionbal.PhotoRate_Shell[ipCARBON][0][1][0]; 00089 00090 /* >>chng 03 sep 29, synch up ion and co solvers. 00091 * the co solver will have a different C/C+ balance that 00092 * is strongly affected by the chemistry if we are in neutral gas 00093 * in this case use co solver's C/C+ balance */ 00094 save_rec = ionbal.RateRecomTot[ipCARBON][0]; 00095 00096 bool lgDEBUG=false; 00097 /*if( iteration>1 )*/ 00098 /*if( nzone > 258 ) 00099 lgDEBUG = true; 00100 else 00101 lgDEBUG = false;*/ 00102 ion_solver(ipCARBON,lgDEBUG); 00103 00104 /* reset the var we just hosed */ 00105 if( save_rec > 0. ) 00106 ionbal.RateRecomTot[ipCARBON][0] = save_rec; 00107 00108 /* rate will be cm-3 s-1, into triplets only 00109 * >>chng 96 nov 22, rid of carb() */ 00110 carb.p1909 *= dense.xIonDense[ipCARBON][1]*0.62; 00111 carb.p2326 *= dense.xIonDense[ipCARBON][0]*0.1; 00112 00113 if( trace.lgTrace ) 00114 { 00115 fprintf( ioQQQ, " IonCarbo returns; fracs=" ); 00116 for( int i=0; i < 7; i++ ) 00117 { 00118 fprintf( ioQQQ, " %10.3e", dense.xIonDense[ipCARBON][i]/ 00119 dense.gas_phase[ipCARBON] ); 00120 } 00121 fprintf( ioQQQ, "\n" ); 00122 } 00123 00124 /* We have already checked the validity of cont_gaunt_calc in sanitycheck.c. 00125 * Now we check to see if the InterpolateGff routine also works correctly. */ 00126 enum {AGN=false}; 00127 /* if set true there will be two problems at very low temps */ 00128 if( AGN ) 00129 { 00130 phycon.te=10.; 00131 /* this tells ion_recomb to punch data */ 00132 punch.lgioRecom = true; 00133 /* this is where it will be punched */ 00134 punch.ioRecom = ioQQQ; 00135 while( phycon.te<1e7 ) 00136 { 00137 /* get recombination coefficients */ 00138 /*lint -e64 double = pointer */ 00139 ion_recomb(false,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc, 00140 dd,ff,ipCARBON); 00141 /*lint +e64 double = pointer */ 00142 TempChange(phycon.te *2.f , true); 00143 } 00144 /* bail out */ 00145 cdEXIT(EXIT_SUCCESS); 00146 } 00147 return; 00148 }