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 /*iso_photo do photoionization rates for element nelem on the ipISO isoelectronic sequence */ 00004 #include "cddefines.h" 00005 #include "hydrogenic.h" 00006 #include "rfield.h" 00007 #include "opacity.h" 00008 #include "trace.h" 00009 #include "ionbal.h" 00010 #include "thermal.h" 00011 #include "gammas.h" 00012 #include "iso.h" 00013 00014 void iso_photo( 00015 /* iso sequence, hydrogen or helium for now */ 00016 long ipISO , 00017 /* the chemical element, 0 for hydrogen */ 00018 long int nelem) 00019 { 00020 long int limit , 00021 n; 00022 00023 DEBUG_ENTRY( "iso_photo()" ); 00024 00025 /* check that we were called with valid charge */ 00026 ASSERT( nelem >= 0 && nelem < LIMELM ); 00027 ASSERT( ipISO < NISO ); 00028 00029 /* do photoionization rates */ 00030 /* induced recombination; FINDUC is integral of 00031 * pho rate times EXP(-hn/kt) for induc rec 00032 * CIND is this times hnu-hnu0 to get ind rec cooling 00033 * ionbal.lgPhotoIoniz_On is 1, set to 0 with 'no photoionization' command 00034 * ipSecIon points to 7.353 Ryd, lowest energy where secondary ioniz 00035 * of hydrogen is possible */ 00036 00037 /* photoionization of ground, this is different from excited states because 00038 * there will eventually be more than one shell, (when li-like done) 00039 * and because upper limit is high-energy limit to code, so secondaries are 00040 * included */ 00041 iso.gamnc[ipISO][nelem][0] = GammaBn(iso.ipIsoLevNIonCon[ipISO][nelem][0], 00042 rfield.nflux, 00043 iso.ipOpac[ipISO][nelem][0], 00044 iso.xIsoLevNIonRyd[ipISO][nelem][0], 00045 &iso.RecomInducRate[ipISO][nelem][0], 00046 &iso.RecomInducCool_Coef[ipISO][nelem][0])* 00047 ionbal.lgPhotoIoniz_On; 00048 00049 /* heating due to photo of ground */ 00050 iso.PhotoHeat[ipISO][nelem][0] = thermal.HeatNet*ionbal.lgPhotoIoniz_On; 00051 00052 /* save these rates into ground photo rate vector */ 00053 ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][0] = iso.gamnc[ipISO][nelem][ipH1s]; 00054 ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][1] = thermal.HeatLowEnr*ionbal.lgPhotoIoniz_On; 00055 ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][2] = thermal.HeatHiEnr*ionbal.lgPhotoIoniz_On; 00056 00057 /* CompRecoilIonRate is direct photioniz rate due to 00058 * bound compton scattering of very hard x-rays+Compton scat */ 00059 /* now add in compton recoil, to ground state, save heating as high energy heat */ 00060 ASSERT( ionbal.CompRecoilIonRate[nelem][nelem-ipISO]>=0. && 00061 ionbal.CompRecoilHeatRate[nelem][nelem-ipISO]>= 0. ); 00062 iso.gamnc[ipISO][nelem][0] += ionbal.CompRecoilIonRate[nelem][nelem-ipISO]; 00063 iso.PhotoHeat[ipISO][nelem][0] += ionbal.CompRecoilHeatRate[nelem][nelem-ipISO]; 00064 00065 /* now add bound compton scattering to ground term photoionization rate */ 00066 ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][0] += ionbal.CompRecoilIonRate[nelem][nelem-ipISO]; 00067 /* add heat to high energy heating term */ 00068 ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][2] += ionbal.CompRecoilHeatRate[nelem][nelem-ipISO]; 00069 00070 /* option to print ground state photoionization rates */ 00071 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) ) 00072 { 00073 GammaPrt(iso.ipIsoLevNIonCon[ipISO][nelem][0], 00074 rfield.nflux, 00075 iso.ipOpac[ipISO][nelem][0], 00076 ioQQQ, 00077 iso.gamnc[ipISO][nelem][0], 00078 iso.gamnc[ipISO][nelem][0]*0.05); 00079 } 00080 00081 /* for excited states upper limit to integration is ground threshold */ 00082 limit = iso.ipIsoLevNIonCon[ipISO][nelem][0]-1; 00083 /* >>chng 06 aug 17, to numLevels_local instead of _max. */ 00084 for( n=1; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00085 { 00086 /* continuously update rates for n <=6, but only update 00087 * rates for higher levels when redoing static opacities */ 00088 if( StatesElem[ipISO][nelem][n].n>6 && !opac.lgRedoStatic ) 00089 break; 00094 if( hydro.lgHInducImp ) 00095 { 00096 iso.gamnc[ipISO][nelem][n] = 00097 GammaBn( 00098 iso.ipIsoLevNIonCon[ipISO][nelem][n], 00099 limit, 00100 iso.ipOpac[ipISO][nelem][n], 00101 iso.xIsoLevNIonRyd[ipISO][nelem][n], 00102 &iso.RecomInducRate[ipISO][nelem][n], 00103 &iso.RecomInducCool_Coef[ipISO][nelem][n])* 00104 ionbal.lgPhotoIoniz_On; 00105 } 00106 else 00107 { 00108 iso.gamnc[ipISO][nelem][n] = 00109 GammaK(iso.ipIsoLevNIonCon[ipISO][nelem][n], 00110 limit, 00111 iso.ipOpac[ipISO][nelem][n],1.)* 00112 ionbal.lgPhotoIoniz_On; 00113 00114 /* these are zero */ 00115 iso.RecomInducRate[ipISO][nelem][n] = 0.; 00116 iso.RecomInducCool_Coef[ipISO][nelem][n] = 0.; 00117 } 00118 iso.PhotoHeat[ipISO][nelem][n] = thermal.HeatNet*ionbal.lgPhotoIoniz_On; 00119 00120 ASSERT( iso.gamnc[ipISO][nelem][n]>= 0. ); 00121 ASSERT( iso.PhotoHeat[ipISO][nelem][n]>= 0. ); 00122 /* this loop only has excited states */ 00123 } 00124 00125 { 00126 /*@-redef@*/ 00127 enum {DEBUG_LOC=false}; 00128 /*@+redef@*/ 00129 if( DEBUG_LOC ) 00130 { 00131 if( nelem==ipHYDROGEN ) 00132 { 00133 fprintf(ioQQQ," buggbugg hphotodebugg%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00134 nzone, 00135 iso.gamnc[ipISO][nelem][0], 00136 iso.gamnc[ipISO][nelem][1], 00137 iso.gamnc[ipISO][nelem][3], 00138 iso.gamnc[ipISO][nelem][4], 00139 iso.gamnc[ipISO][nelem][5], 00140 iso.gamnc[ipISO][nelem][6]); 00141 } 00142 } 00143 } 00144 00145 /* >>chng 02 jan 19, kill excited state photoionization with case b no photo */ 00146 /* option for case b conditions, kill all excited state photoionizations */ 00147 if( opac.lgCaseB_no_photo ) 00148 { 00149 for( n=1; n < iso.numLevels_max[ipISO][nelem]; n++ ) 00150 { 00151 iso.gamnc[ipISO][nelem][n] = 0.; 00152 iso.RecomInducRate[ipISO][nelem][n] = 0.; 00153 iso.RecomInducCool_Coef[ipISO][nelem][n] = 0.; 00154 } 00155 } 00156 { 00157 /* this block turns off induced recom for some element */ 00158 /*@-redef@*/ 00159 enum {DEBUG_LOC=false}; 00160 /*@+redef@*/ 00161 if( DEBUG_LOC && ipISO==1 && nelem==5) 00162 { 00163 /* this debugging block is normally not active */ 00164 for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ ) 00165 { 00166 iso.RecomInducRate[ipISO][nelem][n] = 0.; 00167 } 00168 } 00169 } 00170 00171 if( trace.lgTrace ) 00172 { 00173 fprintf( ioQQQ, " iso_photo, ipISO%2ld nelem%2ld low, hi=",ipISO,nelem); 00174 fprintf( ioQQQ,PrintEfmt("%9.2e", iso.gamnc[ipISO][nelem][ipH1s])); 00175 ASSERT(nelem>=ipISO); 00176 fprintf( ioQQQ,PrintEfmt("%9.2e", ionbal.CompRecoilIonRate[nelem][nelem-ipISO])); 00177 fprintf( ioQQQ, " total="); 00178 fprintf( ioQQQ,PrintEfmt("%9.2e",iso.gamnc[ipISO][nelem][ipH1s] )); 00179 fprintf( ioQQQ, "\n"); 00180 } 00181 return; 00182 }