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 /*atoms_fe2ovr compute FeII overlap with Lya */ 00004 /*fe2par evaluate FeII partition function */ 00005 #include "cddefines.h" 00006 #include "doppvel.h" 00007 #include "dense.h" 00008 #include "rfield.h" 00009 #include "iso.h" 00010 #include "phycon.h" 00011 #include "hydrogenic.h" 00012 #include "ipoint.h" 00013 #include "opacity.h" 00014 #include "radius.h" 00015 #include "atomfeii.h" 00016 #include "thirdparty.h" 00017 00018 const double WLAL = 1215.6845; 00019 00020 /* There are 373 transitions: 00021 * Wavelength (A) 00022 * absorption oscillator strength 00023 * Energy of lower level (Ryd) 00024 * Statistical weight of lower level (g) */ 00026 t_fe2ovr_la::t_fe2ovr_la() 00027 { 00028 DEBUG_ENTRY( "t_fe2ovr_la()" ); 00029 00030 const long VERSION_MAGIC = 20070717L; 00031 const static char chFile[] = "fe2ovr_la.dat"; 00032 00033 FILE *io = open_data( chFile, "r" ); 00034 00035 bool lgErr = false; 00036 long i = -1L; 00037 00038 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 ); 00039 if( lgErr || i != VERSION_MAGIC ) 00040 { 00041 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i ); 00042 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC ); 00043 cdEXIT(EXIT_FAILURE); 00044 } 00045 00046 double help=0; 00047 for( i=0; i < NFEII; i++ ) 00048 { 00049 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 ); 00050 fe2lam[i] = (realnum)help; 00051 } 00052 for( i=0; i < NFEII; i++ ) 00053 { 00054 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 ); 00055 fe2osc[i] = (realnum)help; 00056 } 00057 for( i=0; i < NFEII; i++ ) 00058 { 00059 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 ); 00060 fe2enr[i] = (realnum)help; 00061 } 00062 for( i=0; i < NFEII; i++ ) 00063 { 00064 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 ); 00065 fe2gs[i] = (realnum)help; 00066 } 00067 for( i=0; i < NFE2PR; i++ ) 00068 lgErr = lgErr || ( fscanf( io, "%le", &fe2pt[i] ) != 1 ); 00069 for( i=0; i < NFE2PR; i++ ) 00070 lgErr = lgErr || ( fscanf( io, "%le", &fe2pf[i] ) != 1 ); 00071 00072 fclose( io ); 00073 00074 ASSERT( !lgErr ); 00075 return; 00076 } 00077 00078 void t_fe2ovr_la::zero_opacity() 00079 { 00080 DEBUG_ENTRY( "zero_opacity()" ); 00081 00082 for( int i=0; i < NFEII; i++ ) 00083 { 00084 Fe2PopLte[i] = 0.f; 00085 feopc[i] = 0.f; 00086 Fe2TauLte[i] = opac.taumin; 00087 } 00088 return; 00089 } 00090 00091 void t_fe2ovr_la::init_pointers() 00092 { 00093 DEBUG_ENTRY( "init_pointers()" ); 00094 00095 for( int i=0; i < NFEII; i++ ) 00096 ipfe2[i] = ipoint(fe2enr[i]); 00097 return; 00098 } 00099 00101 void t_fe2ovr_la::tau_inc() 00102 { 00103 DEBUG_ENTRY( "tau_inc()" ); 00104 00105 for( int i=0; i < NFEII; i++ ) 00106 /* optical depths for Feii dest of lya when large feii not used */ 00107 Fe2TauLte[i] += feopc[i]*(realnum)radius.drad_x_fillfac; 00108 return; 00109 } 00110 00112 void t_fe2ovr_la::atoms_fe2ovr(void) 00113 { 00114 DEBUG_ENTRY( "atoms_fe2ovr()" ); 00115 00116 long int i; 00117 00118 static long int nZoneEval; 00119 00120 double Fe2Partn, 00121 displa, 00122 hopc, 00123 rate, 00124 weight; 00125 00126 static double BigFeWidth, 00127 BigHWidth, 00128 oldrat; 00129 00130 /* wavelength of Lya in Angstroms */ 00131 00132 /* compute efficiency of FeII emission overlapping with Ly-alpha 00133 * implemented with Fred Hamann 00134 * 00135 * make Ly-a width monotonically increasing to avoid oscillation 00136 * in deep regions of x-ray ionized clouds. 00137 * 00138 * do nothing if large FeII atom is used */ 00139 if( FeII.lgFeIILargeOn ) 00140 { 00141 return; 00142 } 00143 00144 if( nzone <= 1 ) 00145 { 00146 BigHWidth = hydro.HLineWidth; 00147 BigFeWidth = DoppVel.doppler[ipIRON]; 00148 nZoneEval = nzone; 00149 } 00150 00151 /* do not do pumping if no population,line is thin, or turned off */ 00152 if( (dense.xIonDense[ipIRON][1] <= 0. || !FeII.lgLyaPumpOn) || 00153 hydro.HLineWidth <= 0. ) 00154 { 00155 Fe2Partn = 0.; 00156 oldrat = 0.; 00157 hydro.dstfe2lya = 0.; 00158 00159 for( i=0; i < NFEII; i++ ) 00160 { 00161 Fe2PopLte[i] = 0.; 00162 } 00163 return; 00164 } 00165 00166 /* only evaluate this one time per zone to avoid oscillations 00167 * deep in x-ray ionized clouds */ 00168 if( nZoneEval == nzone && nzone > 1 ) 00169 { 00170 return; 00171 } 00172 00173 BigHWidth = MAX2(BigHWidth,(double)hydro.HLineWidth); 00174 BigFeWidth = MAX2(BigFeWidth,(double)DoppVel.doppler[ipIRON]); 00175 nZoneEval = nzone; 00176 00177 /* check that data is linked in */ 00178 ASSERT( fe2lam[0] > 0. ); 00179 00180 rate = 0.; 00181 Fe2Partn = fe2par(phycon.te); 00182 for( i=0; i < NFEII; i++ ) 00183 { 00184 /* this is displacement from line center in units of Lya width */ 00185 displa = fabs(fe2lam[i]-WLAL)/WLAL*3e10/BigHWidth; 00186 if( displa < 1.333 ) 00187 { 00188 /* have variable weighting factor depending on distance away 00189 * this comes form the Verner's fits to Adam's results */ 00190 weight = ( displa <= 0.66666 ) ? 1. : MAX2(0.,1.-(displa-0.666666)/0.66666); 00191 00192 Fe2PopLte[i] = (realnum)(fe2gs[i]/Fe2Partn*rfield.ContBoltz[ipfe2[i]-1]* 00193 dense.xIonDense[ipIRON][1]); 00194 00195 feopc[i] = (realnum)(Fe2PopLte[i]*fe2osc[i]*0.0150*(fe2lam[i]*1e-8)/BigFeWidth); 00196 00197 /* Ly-alpha line-center opacity */ 00198 if( StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1] > 0. ) 00199 { 00200 hopc = 7.60e-8*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop* 00201 dense.xIonDense[ipHYDROGEN][1]/DoppVel.doppler[ipHYDROGEN]; 00202 } 00203 else 00204 { 00205 hopc = 7.60e-8*dense.xIonDense[ipHYDROGEN][0]/DoppVel.doppler[0]; 00206 } 00207 rate += (feopc[i]/SDIV((feopc[i] + hopc)))*(BigFeWidth/DoppVel.doppler[ipHYDROGEN])* 00208 (1. - 1./(1. + 1.6*Fe2TauLte[i]))*weight; 00209 } 00210 } 00211 00212 /* dstfe2lya is total Lya deexcitation probability due to line overlap 00213 * damper is to stop feedback between here and hydrogen ionization sln */ 00216 hydro.dstfe2lya = (realnum)((rate + oldrat)/2.); 00217 oldrat = hydro.dstfe2lya; 00218 return; 00219 } 00220 00222 double t_fe2ovr_la::fe2par(double te) 00223 { 00224 DEBUG_ENTRY( "fe2par()" ); 00225 00226 double fe2par_v; 00227 00228 /* function to evaluate partition function for FeII 00229 * 00230 * Temperature (K) */ 00231 00232 if( te <= fe2pt[0] ) 00233 fe2par_v = fe2pf[0]; 00234 else if( te >= fe2pt[NFE2PR-1] ) 00235 fe2par_v = fe2pf[NFE2PR-1]; 00236 else 00237 { 00238 long i = hunt_bisect(fe2pt,NFE2PR,te); 00239 double slope = (fe2pf[i+1] - fe2pf[i])/(fe2pt[i+1] - fe2pt[i]); 00240 double du = slope*(te - fe2pt[i]); 00241 fe2par_v = fe2pf[i] + du; 00242 } 00243 return( fe2par_v ); 00244 }