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
00018 const double WLAL = 1215.6845;
00019
00020
00021
00022
00023
00024
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
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
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 = DoppVel.doppler[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 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
00167
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
00178 ASSERT( fe2lam[0] > 0. );
00179
00180 rate = 0.;
00181 Fe2Partn = fe2par(phycon.te);
00182 for( i=0; i < NFEII; i++ )
00183 {
00184
00185 displa = fabs(fe2lam[i]-WLAL)/WLAL*3e10/BigHWidth;
00186 if( displa < 1.333 )
00187 {
00188
00189
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
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
00213
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
00229
00230
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 }