00001
00002
00003
00004
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 #include "taulines.h"
00018
00019 const double WLAL = 1215.6845;
00020
00021
00022
00023
00024
00025
00027 t_fe2ovr_la::t_fe2ovr_la()
00028 {
00029 DEBUG_ENTRY( "t_fe2ovr_la()" );
00030
00031 const long VERSION_MAGIC = 20070717L;
00032 static const char chFile[] = "fe2ovr_la.dat";
00033
00034 FILE *io = open_data( chFile, "r" );
00035
00036 bool lgErr = false;
00037 long i = -1L;
00038
00039 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
00040 if( lgErr || i != VERSION_MAGIC )
00041 {
00042 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i );
00043 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00044 cdEXIT(EXIT_FAILURE);
00045 }
00046
00047 double help=0;
00048 for( i=0; i < NFEII; i++ )
00049 {
00050 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
00051 fe2lam[i] = (realnum)help;
00052 }
00053 for( i=0; i < NFEII; i++ )
00054 {
00055 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
00056 fe2osc[i] = (realnum)help;
00057 }
00058 for( i=0; i < NFEII; i++ )
00059 {
00060 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
00061 fe2enr[i] = (realnum)help;
00062 }
00063 for( i=0; i < NFEII; i++ )
00064 {
00065 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
00066 fe2gs[i] = (realnum)help;
00067 }
00068 for( i=0; i < NFE2PR; i++ )
00069 lgErr = lgErr || ( fscanf( io, "%le", &fe2pt[i] ) != 1 );
00070 for( i=0; i < NFE2PR; i++ )
00071 lgErr = lgErr || ( fscanf( io, "%le", &fe2pf[i] ) != 1 );
00072
00073 fclose( io );
00074
00075 ASSERT( !lgErr );
00076 return;
00077 }
00078
00079 void t_fe2ovr_la::zero_opacity()
00080 {
00081 DEBUG_ENTRY( "zero_opacity()" );
00082
00083 for( int i=0; i < NFEII; i++ )
00084 {
00085 Fe2PopLte[i] = 0.f;
00086 feopc[i] = 0.f;
00087 Fe2TauLte[i] = opac.taumin;
00088 }
00089 return;
00090 }
00091
00092 void t_fe2ovr_la::init_pointers()
00093 {
00094 DEBUG_ENTRY( "init_pointers()" );
00095
00096 for( int i=0; i < NFEII; i++ )
00097 ipfe2[i] = ipoint(fe2enr[i]);
00098 return;
00099 }
00100
00102 void t_fe2ovr_la::tau_inc()
00103 {
00104 DEBUG_ENTRY( "tau_inc()" );
00105
00106 for( int i=0; i < NFEII; i++ )
00107
00108 Fe2TauLte[i] += feopc[i]*(realnum)radius.drad_x_fillfac;
00109 return;
00110 }
00111
00113 void t_fe2ovr_la::atoms_fe2ovr(void)
00114 {
00115 DEBUG_ENTRY( "atoms_fe2ovr()" );
00116
00117 long int i;
00118
00119 static long int nZoneEval;
00120
00121 double Fe2Partn,
00122 displa,
00123 hopc,
00124 rate,
00125 weight;
00126
00127 static double BigFeWidth,
00128 BigHWidth;
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139 if( FeII.lgFeIILargeOn )
00140 {
00141 return;
00142 }
00143
00144 if( nzone <= 1 )
00145 {
00146 BigHWidth = hydro.HLineWidth;
00147 BigFeWidth = GetDopplerWidth(dense.AtomicWeight[ipIRON]);
00148 nZoneEval = nzone;
00149 }
00150
00151
00152 if( (dense.xIonDense[ipIRON][1] <= 0. || !FeII.lgLyaPumpOn) ||
00153 hydro.HLineWidth <= 0. )
00154 {
00155 Fe2Partn = 0.;
00156 hydro.dstfe2lya = 0.;
00157
00158 for( i=0; i < NFEII; i++ )
00159 {
00160 Fe2PopLte[i] = 0.;
00161 }
00162 return;
00163 }
00164
00165
00166
00167 if( nZoneEval == nzone && nzone > 1 )
00168 {
00169 return;
00170 }
00171
00172 BigHWidth = MAX2(BigHWidth,(double)hydro.HLineWidth);
00173 BigFeWidth = MAX2(BigFeWidth,(double)GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
00174 nZoneEval = nzone;
00175
00176
00177 ASSERT( fe2lam[0] > 0. );
00178
00179 rate = 0.;
00180 Fe2Partn = fe2par(phycon.te);
00181 for( i=0; i < NFEII; i++ )
00182 {
00183
00184 displa = fabs(fe2lam[i]-WLAL)/WLAL*3e10/BigHWidth;
00185 if( displa < 1.333 )
00186 {
00187
00188
00189 weight = ( displa <= 0.66666 ) ? 1. : MAX2(0.,1.-(displa-0.666666)/0.66666);
00190
00191 Fe2PopLte[i] = (realnum)(fe2gs[i]/Fe2Partn*rfield.ContBoltz[ipfe2[i]-1]*
00192 dense.xIonDense[ipIRON][1]);
00193
00194 feopc[i] = (realnum)(Fe2PopLte[i]*fe2osc[i]*0.0150*(fe2lam[i]*1e-8)/BigFeWidth);
00195
00196
00197 if( StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop > 0. )
00198 {
00199 hopc = 7.60e-8*StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop/
00200 GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]);
00201 }
00202 else
00203 {
00204 hopc = 7.60e-8*dense.xIonDense[ipHYDROGEN][0]/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]);
00205 }
00206 rate += (feopc[i]/SDIV((feopc[i] + hopc)))*(BigFeWidth/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]))*
00207 (1. - 1./(1. + 1.6*Fe2TauLte[i]))*weight;
00208 }
00209 }
00210
00211
00212 hydro.dstfe2lya = (realnum)rate;
00213 return;
00214 }
00215
00217 double t_fe2ovr_la::fe2par(double te)
00218 {
00219 DEBUG_ENTRY( "fe2par()" );
00220
00221 double fe2par_v;
00222
00223
00224
00225
00226
00227 if( te <= fe2pt[0] )
00228 fe2par_v = fe2pf[0];
00229 else if( te >= fe2pt[NFE2PR-1] )
00230 fe2par_v = fe2pf[NFE2PR-1];
00231 else
00232 {
00233 long i = hunt_bisect(fe2pt,NFE2PR,te);
00234 double slope = (fe2pf[i+1] - fe2pf[i])/(fe2pt[i+1] - fe2pt[i]);
00235 double du = slope*(te - fe2pt[i]);
00236 fe2par_v = fe2pf[i] + du;
00237 }
00238 return( fe2par_v );
00239 }