00001
00002
00003
00004 #include "cddefines.h"
00005 #include "taulines.h"
00006 #include "atomfeii.h"
00007 #include "dense.h"
00008 #include "conv.h"
00009 #include "atoms.h"
00010 #include "rfield.h"
00011 #include "wind.h"
00012 #include "iso.h"
00013 #include "h2.h"
00014 #include "opacity.h"
00015 #include "trace.h"
00016 #include "lines_service.h"
00017 #include "atmdat.h"
00018 #include "hydrogenic.h"
00019 #include "rt.h"
00020 #include "cosmology.h"
00021 #include "physconst.h"
00022 #include "doppvel.h"
00023 #include "mole.h"
00024
00025
00026 void RT_line_all( void )
00027 {
00028 long int i,
00029 ion,
00030 ipISO,
00031 nelem;
00032 long ipHi , ipLo;
00033
00034
00035 bool lgPescUpdate = conv.lgFirstSweepThisZone || conv.lgIonStageTrimed;
00036
00037 DEBUG_ENTRY( "RT_line_all()" );
00038
00039 if( trace.lgTrace )
00040 fprintf( ioQQQ, " RT_line_all called\n" );
00041
00042
00043
00044 if( !rfield.lgDoLineTrans && conv.nPres2Ioniz )
00045 {
00046 return;
00047 }
00048
00049
00050
00051 if( conv.lgLastSweepThisZone )
00052 {
00053
00054 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
00055
00056 if( 0 && cosmology.lgDo )
00057 {
00058 realnum dVel = (realnum)SPEEDLIGHT * ( 1.f - 1.f/POW2(1.f+cosmology.redshift_start-cosmology.redshift_current) );
00059 rfield.ipFineConVelShift = -(long int)( dVel/ rfield.fine_opac_velocity_width + 0.5 );
00060 }
00061 else
00062 {
00063
00064
00065
00066
00067
00068 realnum dWind = wind.windv - wind.windv0;
00069
00070
00071
00072 rfield.ipFineConVelShift = -(long int)(dWind / rfield.fine_opac_velocity_width+0.5);
00073 }
00074 }
00075
00076 # if 0
00077
00078
00079
00080 {
00081 #include "doppvel.h"
00082 double dTau , aa;
00083 ipISO = 0; nelem = 0;ipLo = 0;
00084 ipHi = 23;
00085 dTau = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() *
00086 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() /
00087 GetDopplerWidth(dense.AtomicWeight[nelem]) + opac.opacity_abs[iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1];
00088 dTau *= radius.drad_x_fillfac;
00089 aa = log(1. + dTau ) / SDIV(dTau);
00090 fprintf(ioQQQ,"DEBUG dTau\t%.2f\t%.5e\t%.5e\t%.5e\n",fnzone,dTau,
00091 radius.drad_x_fillfac,
00092 aa);
00093 }
00094 # endif
00095
00096
00097 if( lgPescUpdate )
00098 RT_stark();
00099
00100
00101
00102 for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
00103 {
00104
00105 for( nelem=ipISO; nelem < LIMELM; ++nelem )
00106 {
00107
00108
00109 ion = nelem+1-ipISO;
00110
00111
00112 if( !dense.lgElmtOn[nelem] )
00113 continue;
00114
00115 if( ion <= dense.IonHigh[nelem] )
00116 {
00117
00118 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ++ipHi )
00119 {
00120 for( ipLo=0; ipLo < ipHi; ++ipLo )
00121 {
00122
00123
00124 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() < 0 )
00125 continue;
00126
00127
00128
00129 fixit();
00130 RT_line_one( iso_sp[ipISO][nelem].trans(ipHi,ipLo),
00131 true,(realnum)iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up,
00132 GetDopplerWidth(dense.AtomicWeight[nelem]));
00133
00134
00135 enum {DEBUG_LOC=false};
00136 if( DEBUG_LOC && nelem==1&& ipLo==0 )
00137 {
00138 fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n",
00139 ipHi ,
00140 fnzone,
00141 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest());
00142 }
00143 }
00144 }
00145
00146
00147
00148 if( opac.lgCaseB_no_pdest )
00149 {
00150 ipLo = 0;
00151
00152 for( ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
00153 {
00154 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
00155 continue;
00156
00157
00158 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest() = SMALLFLOAT;
00159 }
00160 }
00161
00162 ipLo = 0;
00163
00164
00165
00166 if( dense.xIonDense[nelem][ion] > 1e-30 && conv.lgLastSweepThisZone )
00167 {
00168 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
00169 {
00170 TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
00171
00172 (*tr).Emis().PopOpc() = iso_sp[ipISO][nelem].st[0].Pop();
00173 (*(*tr).Lo()).Pop() =
00174 iso_sp[ipISO][nelem].st[ipLo].Pop();
00175
00176
00177 RT_line_one( *tr, true, 0.f,
00178 GetDopplerWidth(dense.AtomicWeight[nelem]));
00179 }
00180 }
00181
00182 }
00183 }
00184 }
00185
00186
00187
00188
00189 if( lgTauGood( iso_sp[ipH_LIKE][ipHYDROGEN].trans(iso_ctrl.nLyaLevel[ipH_LIKE],0) ) )
00190 {
00191
00192
00193
00194
00195
00196
00197
00198 t_fe2ovr_la::Inst().atoms_fe2ovr();
00199
00200
00201
00202
00203
00204
00205 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest() += hydro.dstfe2lya;
00206 }
00207
00208
00209
00210 if( !hydro.lgLymanPumping )
00211 {
00212 ipISO = ipH_LIKE;
00213 nelem = ipHYDROGEN;
00214 ipLo = 0;
00215 for( ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
00216 {
00217 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump() = 0.;
00218 }
00219 }
00220
00221 {
00222
00223 enum {DEBUG_LOC=false};
00224 if( DEBUG_LOC && nzone>433 )
00225 {
00226
00227 DumpLine(iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0) );
00228 }
00229 }
00230
00231
00232 for( i=1; i <= nLevel1; i++ )
00233 {
00234 RT_line_one( TauLines[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
00235 }
00236
00237 {
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 static realnum P_elec_esc_ipTFe56;
00249 static long nSave;
00250 if( nzone <= 1 || nzone!=nSave )
00251 {
00252 P_elec_esc_ipTFe56 = TauLines[ipTFe56].Emis().Pelec_esc();
00253 nSave = nzone;
00254 }
00255 else
00256 TauLines[ipTFe56].Emis().Pelec_esc() = P_elec_esc_ipTFe56;
00257 }
00258
00259
00260 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
00261 {
00262 if( dBaseSpecies[ipSpecies].lgActive )
00263 {
00264 realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
00265 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
00266 tr != dBaseTrans[ipSpecies].end(); ++tr)
00267 {
00268 int ipHi = (*tr).ipHi();
00269 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
00270 continue;
00271 if( (*tr).ipCont() > 0 )
00272 {
00273 RT_line_one( *tr, true, 0.f, DopplerWidth );
00274 }
00275 }
00276 }
00277 }
00278
00279
00280
00281
00282 if( conv.lgFirstSweepThisZone || conv.lgLastSweepThisZone )
00283 {
00284 for( i=0; i < nUTA; i++ )
00285 {
00286
00287 UTALines[i].Emis().PopOpc() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
00288 (*UTALines[i].Lo()).Pop() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
00289 (*UTALines[i].Hi()).Pop() = 0.;
00290 RT_line_one( UTALines[i], true,0.f,
00291 GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
00292 }
00293 }
00294
00295
00296
00297 for( ipISO=ipHE_LIKE; ipISO<NISO; ++ipISO )
00298 {
00299
00300 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00301 {
00302 if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
00303 {
00304 for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ++ipLo )
00305 {
00306 RT_line_one( SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]], true,0.f,
00307 GetDopplerWidth(dense.AtomicWeight[nelem]) );
00308 }
00309 }
00310 }
00311 }
00312
00313
00314 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00315 (*diatom)->H2_RTMake( );
00316
00317 return;
00318 }