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.numLevels_max[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
00039 {
00040 fprintf( ioQQQ, "%3ld", ipHi );
00041 for( ipLo=0; ipLo < ipHi; ipLo++ )
00042 {
00043 if( Transitions[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].ipCont <= 0 )
00044 fprintf( ioQQQ, "%10.2e", 1e-30 );
00045 else
00046 fprintf( ioQQQ, "%10.2e",
00047 Transitions[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][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.num_comole_calc; i++ )
00060 {
00061 COmole[i]->hevcol_old = COmole[i]->hevcol;
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.numLevels_max[ipISO][nelem]; ipHi++ )
00075 {
00076 if( iso.lgDielRecom[ipISO] )
00077 RT_line_one_tau_reset(&SatelliteLines[ipISO][nelem][ipHi]);
00078 }
00079
00080 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; 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.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 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot ,
00092 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn);
00093 }
00094
00095
00096 RT_line_one_tau_reset(&Transitions[ipISO][nelem][ipHi][ipLo]);
00097 if( DEBUG_LOC )
00098 {
00099 if( ipISO==1 && nelem==1 && ipHi==iso.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 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot ,
00103 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn);
00104 }
00105 }
00106 }
00107
00108 for( ipHi=StatesElemNEW[nelem][nelem-ipISO][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ )
00109 {
00110
00111
00112 RT_line_one_tau_reset(&ExtraLymanLines[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 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn = opac.tlamin;
00130
00131 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauCon = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn;
00132 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot =
00133 2.f*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn;
00134 f = opac.tlamin/Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
00135
00136 for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00137 {
00138 if( Transitions[ipH_LIKE][nelem][ipHi][ipH1s].ipCont <= 0 )
00139 continue;
00140
00141 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn =
00142 f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity;
00143
00144 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauCon = Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn;
00145 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot =
00146 2.f*Transitions[ipH_LIKE][nelem][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 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn = opac.tlamin;
00161
00162 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauCon = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
00163
00164 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauTot =
00165 2.f*Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
00166
00167 ratio = opac.tlamin/Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->opacity;
00168
00169
00170
00171 Aprev = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->Aul;
00172
00173 i = ipHe2p1P+1;
00174
00175
00176 for( i = ipHe2p1P+1; i < iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
00177 {
00178 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 )
00179 continue;
00180
00181
00182 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul> Aprev/10. ||
00183 StatesElemNEW[nelem][nelem-ipHE_LIKE][i].S < 0 )
00184 {
00185 Aprev = Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul;
00186 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn =
00187 ratio*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->opacity;
00188
00189 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauCon =
00190 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
00191 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauTot =
00192 2.f*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
00193
00194
00195 }
00196 }
00197 }
00198 }
00199 }
00200
00201
00202 for( i=1; i <= nLevel1; i++ )
00203 {
00204 RT_line_one_tau_reset(&TauLines[i]);
00205
00206 ASSERT( fabs(TauLines[i].Emis->TauIn) <= fabs(TauLines[i].Emis->TauTot) );
00207 }
00208
00209
00210 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
00211
00212
00213 for( i=0; i < nWindLine; i++ )
00214 {
00215 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
00216 {
00217 RT_line_one_tau_reset(&TauLine2[i]);
00218 }
00219 }
00220
00221
00222 for( i=0; i < nHFLines; i++ )
00223 {
00224 RT_line_one_tau_reset(&HFLines[i]);
00225 }
00226
00227
00228 for( i=0; i < linesAdded2; i++)
00229 {
00230 RT_line_one_tau_reset(dBaseLines[i].tran);
00231 }
00232
00233
00234 for( i=0; i < nUTA; i++ )
00235 {
00236
00237
00238
00239 double hsave = UTALines[i].Coll.heat;
00240 RT_line_one_tau_reset(&UTALines[i]);
00241 UTALines[i].Coll.heat = hsave;
00242 }
00243
00244
00245 H2_RT_tau_reset();
00246
00247
00248 FeII_RT_tau_reset();
00249
00250 if( opac.lgCaseB )
00251 {
00252 for( i=0; i < rfield.nupper; i++ )
00253 {
00254
00255
00256
00257 opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]/2.f;
00258
00259 opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]/2.f;
00260 }
00261 }
00262 else if( geometry.lgSphere )
00263 {
00264
00265 for( i=0; i < rfield.nupper; i++ )
00266 {
00267
00268
00269 opac.TauAbsGeo[1][i] = 2.f*opac.TauAbsFace[i];
00270 opac.TauAbsGeo[0][i] = opac.TauAbsFace[i];
00271 opac.TauScatGeo[1][i] = 2.f*opac.TauScatFace[i];
00272 opac.TauScatGeo[0][i] = opac.TauScatFace[i];
00273 opac.TauTotalGeo[1][i] = opac.TauScatGeo[1][i] + opac.TauAbsGeo[1][i];
00274 opac.TauTotalGeo[0][i] = opac.TauScatGeo[0][i] + opac.TauAbsGeo[0][i];
00275 }
00276 }
00277 else
00278 {
00279
00280 for( i=0; i < rfield.nupper; i++ )
00281 {
00282 opac.TauTotalGeo[1][i] = opac.TauTotalGeo[0][i];
00283 opac.TauTotalGeo[0][i] = opac.taumin;
00284 opac.TauAbsGeo[1][i] = opac.TauAbsGeo[0][i];
00285 opac.TauAbsGeo[0][i] = opac.taumin;
00286 opac.TauScatGeo[1][i] = opac.TauScatGeo[0][i];
00287 opac.TauScatGeo[0][i] = opac.taumin;
00288 }
00289 }
00290
00291
00292 for( i=0; i < rfield.nupper; i++ )
00293 {
00294
00295 opac.TauAbsTotal[i] = opac.TauAbsFace[i];
00296
00297
00298 opac.E2TauAbsTotal[i] = (realnum)e2( opac.TauAbsTotal[i] );
00299
00300 opac.TauScatFace[i] = opac.taumin;
00301 opac.TauAbsFace[i] = opac.taumin;
00302 }
00303
00304
00305 rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];
00306 return;
00307 }