00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "taulines.h"
00007 #include "trace.h"
00008 #include "iso.h"
00009 #include "rfield.h"
00010 #include "opacity.h"
00011 #include "h2.h"
00012 #include "mole.h"
00013 #include "geometry.h"
00014 #include "dense.h"
00015 #include "atomfeii.h"
00016 #include "colden.h"
00017 #include "rt.h"
00018
00019
00020
00021
00022 void RT_tau_reset(void)
00023 {
00024 long int i,
00025 ipHi,
00026 ipISO,
00027 nelem,
00028 ipLo;
00029
00030 DEBUG_ENTRY( "RT_tau_reset()" );
00031
00032 if( trace.lgTrace )
00033 {
00034 fprintf( ioQQQ, " UPDATE estimating new optical depths\n" );
00035 if( trace.lgHBug && trace.lgIsoTraceFull[ipH_LIKE] )
00036 {
00037 fprintf( ioQQQ, " New Hydrogen outward optical depths:\n" );
00038 for( ipHi=1; ipHi < iso_sp[ipH_LIKE][ trace.ipIsoTrace[ipH_LIKE] ].numLevels_max; ipHi++ )
00039 {
00040 fprintf( ioQQQ, "%3ld", ipHi );
00041 for( ipLo=0; ipLo < ipHi; ipLo++ )
00042 {
00043 if( iso_sp[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]].trans(ipHi,ipLo).ipCont() <= 0 )
00044 fprintf( ioQQQ, "%10.2e", 1e-30 );
00045 else
00046 fprintf( ioQQQ, "%10.2e",
00047 iso_sp[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]].trans(ipHi,ipLo).Emis().TauIn() );
00048 }
00049 fprintf( ioQQQ, "\n" );
00050 }
00051 }
00052 }
00053
00054
00055 for( i=0; i<NCOLD; ++i )
00056 {
00057 colden.colden_old[i] = colden.colden[i];
00058 }
00059 for( i=0; i < mole_global.num_calc; i++ )
00060 {
00061 mole.species[i].column_old = mole.species[i].column;
00062 }
00063
00064 opac.telec = opac.taumin;
00065 opac.thmin = opac.taumin;
00066
00067
00068 for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00069 {
00070 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00071 {
00072 if( dense.lgElmtOn[nelem] )
00073 {
00074 for( ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00075 {
00076 if( iso_ctrl.lgDielRecom[ipISO] )
00077 RT_line_one_tau_reset(SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipHi]]);
00078 }
00079
00080 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00081 {
00082
00083 for( ipLo=0; ipLo < ipHi; ipLo++ )
00084 {
00085 enum {DEBUG_LOC=false};
00086 if( DEBUG_LOC )
00087 {
00088 if( ipISO==1 && nelem==1 && ipHi==iso_ctrl.nLyaLevel[ipISO] && ipLo==0 )
00089 fprintf(ioQQQ,"DEBUG rt before loop %li %li %li %li tot %.3e in %.3e\n",
00090 ipISO, nelem, ipHi , ipLo ,
00091 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() ,
00092 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn());
00093 }
00094
00095
00096 RT_line_one_tau_reset(iso_sp[ipISO][nelem].trans(ipHi,ipLo));
00097 if( DEBUG_LOC )
00098 {
00099 if( ipISO==1 && nelem==1 && ipHi==iso_ctrl.nLyaLevel[ipISO] && ipLo==0 )
00100 fprintf(ioQQQ,"DEBUG rt after loop %li %li %li %li tot %.3e in %.3e\n",
00101 ipISO, nelem, ipHi , ipLo ,
00102 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() ,
00103 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn());
00104 }
00105 }
00106 }
00107
00108 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
00109 {
00110
00111
00112 RT_line_one_tau_reset(ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]]);
00113 }
00114 }
00115 }
00116 }
00117
00118
00119
00120
00121 if( opac.lgCaseB )
00122 {
00123 for( nelem=0; nelem < LIMELM; nelem++ )
00124 {
00125 if( dense.lgElmtOn[nelem] )
00126 {
00127 realnum f;
00128
00129 iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn() = opac.tlamin;
00130
00131 iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauCon() = iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn();
00132 iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauTot() =
00133 2.f*iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn();
00134 f = opac.tlamin/iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().opacity();
00135
00136 for( ipHi=3; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_max; ipHi++ )
00137 {
00138 if( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).ipCont() <= 0 )
00139 continue;
00140
00141 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauIn() =
00142 f*iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().opacity();
00143
00144 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauCon() = iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauIn();
00145 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauTot() =
00146 2.f*iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn();
00147 }
00148 }
00149 }
00150
00151
00152
00153 for( nelem=1; nelem < LIMELM; nelem++ )
00154 {
00155 if( dense.lgElmtOn[nelem] )
00156 {
00157 double Aprev;
00158 realnum ratio;
00159
00160 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauIn() = opac.tlamin;
00161
00162 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauCon() = iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauIn();
00163
00164 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauTot() =
00165 2.f*iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauIn();
00166
00167 ratio = opac.tlamin/iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().opacity();
00168
00169
00170
00171 Aprev = iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().Aul();
00172
00173 i = ipHe2p1P+1;
00174
00175
00176 for( i = ipHe2p1P+1; i < iso_sp[ipHE_LIKE][nelem].numLevels_max; i++ )
00177 {
00178 if( iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).ipCont() <= 0 )
00179 continue;
00180
00181
00182 if( iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().Aul()> Aprev/10. ||
00183 iso_sp[ipHE_LIKE][nelem].st[i].S() < 0 )
00184 {
00185 Aprev = iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().Aul();
00186 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauIn() =
00187 ratio*iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().opacity();
00188
00189 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauCon() =
00190 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauIn();
00191 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauTot() =
00192 2.f*iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauIn();
00193 }
00194 }
00195 }
00196 }
00197 }
00198
00199
00200 for( i=1; i <= nLevel1; i++ )
00201 {
00202 RT_line_one_tau_reset(TauLines[i]);
00203
00204 ASSERT( fabs(TauLines[i].Emis().TauIn()) <= fabs(TauLines[i].Emis().TauTot()) );
00205 }
00206
00207
00208 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
00209
00210
00211 for( i=0; i < nWindLine; i++ )
00212 {
00213 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
00214 {
00215 RT_line_one_tau_reset(TauLine2[i]);
00216 }
00217 }
00218
00219
00220 for( i=0; i < nHFLines; i++ )
00221 {
00222 RT_line_one_tau_reset(HFLines[i]);
00223 }
00224
00225
00226 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
00227 {
00228 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
00229 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
00230 {
00231 RT_line_one_tau_reset((*em).Tran());
00232 }
00233 }
00234
00235
00236 for( i=0; i < nUTA; i++ )
00237 {
00238
00239
00240
00241 double hsave = UTALines[i].Coll().heat();
00242 RT_line_one_tau_reset(UTALines[i]);
00243 UTALines[i].Coll().heat() = hsave;
00244 }
00245
00246
00247 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00248 (*diatom)->H2_RT_tau_reset();
00249
00250
00251 FeII_RT_tau_reset();
00252
00253 if( opac.lgCaseB )
00254 {
00255 for( i=0; i < rfield.nupper; i++ )
00256 {
00257
00258
00259
00260 opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]/2.f;
00261
00262 opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]/2.f;
00263 }
00264 }
00265 else if( geometry.lgSphere )
00266 {
00267
00268 for( i=0; i < rfield.nupper; i++ )
00269 {
00270
00271
00272 opac.TauAbsGeo[1][i] = 2.f*opac.TauAbsFace[i];
00273 opac.TauAbsGeo[0][i] = opac.TauAbsFace[i];
00274 opac.TauScatGeo[1][i] = 2.f*opac.TauScatFace[i];
00275 opac.TauScatGeo[0][i] = opac.TauScatFace[i];
00276 opac.TauTotalGeo[1][i] = opac.TauScatGeo[1][i] + opac.TauAbsGeo[1][i];
00277 opac.TauTotalGeo[0][i] = opac.TauScatGeo[0][i] + opac.TauAbsGeo[0][i];
00278 }
00279 }
00280 else
00281 {
00282
00283 for( i=0; i < rfield.nupper; i++ )
00284 {
00285 opac.TauTotalGeo[1][i] = opac.TauTotalGeo[0][i];
00286 opac.TauTotalGeo[0][i] = opac.taumin;
00287 opac.TauAbsGeo[1][i] = opac.TauAbsGeo[0][i];
00288 opac.TauAbsGeo[0][i] = opac.taumin;
00289 opac.TauScatGeo[1][i] = opac.TauScatGeo[0][i];
00290 opac.TauScatGeo[0][i] = opac.taumin;
00291 }
00292 }
00293
00294
00295 for( i=0; i < rfield.nupper; i++ )
00296 {
00297
00298 opac.TauAbsTotal[i] = opac.TauAbsFace[i];
00299
00300
00301 opac.E2TauAbsTotal[i] = (realnum)e2( opac.TauAbsTotal[i] );
00302
00303 opac.TauScatFace[i] = opac.taumin;
00304 opac.TauAbsFace[i] = opac.taumin;
00305 }
00306
00307
00308 rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];
00309 return;
00310 }