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 /*ion_photo fill array PhotoRate with photoionization rates for heavy elements */ 00004 #include "cddefines.h" 00005 #include "yield.h" 00006 #include "heavy.h" 00007 #include "opacity.h" 00008 #include "dense.h" 00009 #include "thermal.h" 00010 #include "conv.h" 00011 #include "grainvar.h" 00012 #include "elementnames.h" 00013 #include "gammas.h" 00014 #include "ionbal.h" 00015 00016 void ion_photo( 00017 /* nlem is atomic number on C scale, 0 for H */ 00018 long int nelem , 00019 /* debugging flag to turn on print */ 00020 bool lgPrintIt ) 00021 { 00022 long int ion, 00023 iphi, 00024 iplow, 00025 ipop, 00026 limit_hi, 00027 limit_lo, 00028 ns; 00029 00030 DEBUG_ENTRY( "ion_photo()" ); 00031 00032 /* IonLow(nelem) and IonHigh(nelem) are bounds for evaluation*/ 00033 00034 /*begin sanity checks */ 00035 ASSERT( nelem < LIMELM ); 00036 ASSERT( dense.IonLow[nelem] >= 0 ); 00037 ASSERT( dense.IonLow[nelem] <= nelem); 00038 ASSERT( dense.IonHigh[nelem] <= nelem + 1); 00039 /*end sanity checks */ 00040 00041 /* NB - in following, nelem is on c scale, so is 29 for Zn */ 00042 00043 /* min since iso-like atom rates produced in iso_photo. 00044 * IonHigh and IonLow range from 0 to nelem+1, but for photo rates 00045 * we want 0 to nelem since cannot destroy fully stripped ion. 00046 * since iso-seq done elsewhere, want to actually do IonHigh-xxx.*/ 00047 /* >>chng 00 dec 07, logic on limit_hi now precisely identical to ion_solver */ 00048 /* >>chng 02 mar 31, change limit_hi to < in loop (had been <=) and 00049 * also coded for arbitrary number of iso sequences */ 00050 /*limit_hi = MIN2(nelem,dense.IonHigh[nelem]); 00051 limit_hi = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/ 00052 limit_hi = MIN2( dense.IonHigh[nelem] , nelem+1-NISO ); 00053 00054 /* >>chng 03 sep 26, always do atom itself since may be needed for molecules */ 00055 limit_hi = MAX2( 1 , limit_hi ); 00056 00057 /* when grains are present want to use atoms as lower bound to number of stages of ionization, 00058 * since atomic rates needed for species within grains */ 00059 if( !conv.nPres2Ioniz && gv.lgDustOn ) 00060 { 00061 limit_lo = 0; 00062 } 00063 else 00064 { 00065 limit_lo = dense.IonLow[nelem]; 00066 } 00067 00068 /* >>chng 01 dec 11, lower bound now limit_lo */ 00069 /* loop over all ions for this element */ 00070 /* >>chng 02 mar 31, now ion < limit_hi not <= */ 00071 for( ion=limit_lo; ion < limit_hi; ion++ ) 00072 { 00073 /* loop over all shells for this ion */ 00074 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 00075 { 00076 /* always reevaluate the outer shell, and all shells if lgRedoStatic is set */ 00077 if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) ) 00078 { 00079 /* option to redo the rates only on occasion */ 00080 iplow = opac.ipElement[nelem][ion][ns][0]; 00081 iphi = opac.ipElement[nelem][ion][ns][1]; 00082 ipop = opac.ipElement[nelem][ion][ns][2]; 00083 00084 /* compute the photoionization rate, ionbal.lgPhotoIoniz_On is 1, set 0 00085 * with "no photoionization" command */ 00086 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 00087 GammaK(iplow,iphi, 00088 ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0))*ionbal.lgPhotoIoniz_On; 00089 00090 /* these three lines must be kept parallel with the lines 00091 * in GammaK ion*/ 00092 00093 /* the heating rate */ 00094 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = thermal.HeatLowEnr*ionbal.lgPhotoIoniz_On; 00095 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = thermal.HeatHiEnr*ionbal.lgPhotoIoniz_On; 00096 } 00097 } 00098 00099 /* add on compton recoil ionization for atoms to outer shell */ 00100 /* >>chng 02 mar 24, moved here from ion_solver */ 00101 /* this is the outer shell */ 00102 ns = (Heavy.nsShells[nelem][ion]-1); 00103 /* this must be moved to photoionize and have code parallel to iso_photo code */ 00104 ionbal.PhotoRate_Shell[nelem][ion][ns][0] += ionbal.CompRecoilIonRate[nelem][ion]; 00105 /* add the heat as secondary-ionization capable heating */ 00106 ionbal.PhotoRate_Shell[nelem][ion][ns][2] += ionbal.CompRecoilHeatRate[nelem][ion]; 00107 } 00108 00109 /* option to print information about these rates for this element */ 00110 if( lgPrintIt ) 00111 { 00112 /* option to print rates for particular shell */ 00113 ns = 5; 00114 ion = 1; 00115 GammaPrt( 00116 opac.ipElement[nelem][ion][ns][0], 00117 opac.ipElement[nelem][ion][ns][1], 00118 opac.ipElement[nelem][ion][ns][2], 00119 ioQQQ, /* io unit we will write to */ 00120 ionbal.PhotoRate_Shell[nelem][ion][ns][0], 00121 0.05); 00122 00123 /* outer loop is from K to most number of shells present in atom */ 00124 for( ns=0; ns < Heavy.nsShells[nelem][0]; ns++ ) 00125 { 00126 fprintf( ioQQQ, "\n %s", elementnames.chElementNameShort[nelem] ); 00127 fprintf( ioQQQ, " %s" , Heavy.chShell[ns]); 00128 /* MB hydrogenic photo rate may not be included in beow */ 00129 for( ion=0; ion < dense.IonHigh[nelem]; ion++ ) 00130 { 00131 if( Heavy.nsShells[nelem][ion] > ns ) 00132 { 00133 fprintf( ioQQQ, " %8.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] ); 00134 } 00135 else 00136 { 00137 break; 00138 } 00139 } 00140 } 00141 fprintf(ioQQQ,"\n"); 00142 } 00143 return; 00144 }