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 double WeightNew;
00031
00032 DEBUG_ENTRY( "RT_tau_reset()" );
00033
00034 if( trace.lgTrace )
00035 {
00036 fprintf( ioQQQ, " UPDATE estimating new optical depths\n" );
00037 if( trace.lgHBug && trace.lgIsoTraceFull[ipH_LIKE] )
00038 {
00039 fprintf( ioQQQ, " New Hydrogen outward optical depths:\n" );
00040 for( ipHi=1; ipHi < iso.numLevels_max[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
00041 {
00042 fprintf( ioQQQ, "%3ld", ipHi );
00043 for( ipLo=0; ipLo < ipHi; ipLo++ )
00044 {
00045 fprintf( ioQQQ, "%10.2e",
00046 EmisLines[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].TauIn );
00047 }
00048 fprintf( ioQQQ, "\n" );
00049 }
00050 }
00051 }
00052
00053
00054 opac.tpcah[1] = opac.tpcah[0];
00055 opac.tpcah[0] = opac.taumin;
00056
00057
00058 for( i=0; i<NCOLD; ++i )
00059 {
00060 colden.colden_old[i] = colden.colden[i];
00061 }
00062 for( i=0; i < mole.num_comole_calc; i++ )
00063 {
00064 COmole[i]->hevcol_old = COmole[i]->hevcol;
00065 }
00066
00067
00068
00069 if( iteration <= 1 )
00070 {
00071
00072 WeightNew = 1.;
00073 }
00074 else
00075 {
00076 WeightNew = 0.75;
00077 }
00078
00079
00080 opac.lgTauOutOn = true;
00081
00082 opac.telec = opac.taumin;
00083 opac.thmin = opac.taumin;
00084
00085
00086 for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00087 {
00088
00089 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00090 {
00091 if( dense.lgElmtOn[nelem] )
00092 {
00093 for( ipLo=0; ipLo < (iso.numLevels_max[ipISO][nelem] - 1); ipLo++ )
00094 {
00095 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00096 {
00097
00098
00099
00100
00101
00102
00103
00104
00105 RT_line_one_tau_reset(&EmisLines[ipISO][nelem][ipHi][ipLo],WeightNew);
00106 }
00107 }
00108
00109 for( ipHi=2; ipHi <iso.nLyman[ipISO]; ipHi++ )
00110 {
00111
00112
00113
00114 RT_line_one_tau_reset(&iso.ExtraLymanLines[ipISO][nelem][ipHi],WeightNew);
00115 }
00116 }
00117 }
00118 }
00119
00120
00121
00122
00123 if( opac.lgCaseB )
00124 {
00125 for( nelem=0; nelem < LIMELM; nelem++ )
00126 {
00127 if( dense.lgElmtOn[nelem] )
00128 {
00129 float f;
00130
00131 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn = opac.tlamin;
00132
00133 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauCon = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn;
00134 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot =
00135 2.f*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn;
00136 f = opac.tlamin/EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity;
00137
00138 for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00139 {
00140 EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn =
00141 f*EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].opacity;
00142
00143 EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauCon = EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn;
00144 EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauTot =
00145 2.f*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn;
00146 }
00147 }
00148 }
00149
00150
00151
00152 for( nelem=1; nelem < LIMELM; nelem++ )
00153 {
00154 if( dense.lgElmtOn[nelem] )
00155 {
00156 double Aprev;
00157 float ratio;
00158
00159 EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn = opac.tlamin;
00160
00161 EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauCon = EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn;
00162
00163 EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauTot =
00164 2.f*EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn;
00165
00166 ratio = opac.tlamin/EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].opacity;
00167
00168
00169
00170 Aprev = EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Aul;
00171
00172 i = ipHe2p1P+1;
00173
00174
00175 while( i < iso.numLevels_max[ipHE_LIKE][nelem] )
00176
00177 {
00178
00179
00180 if( EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul> Aprev/10. ||
00181 iso.quant_desig[ipHE_LIKE][nelem][i].s < 0 )
00182 {
00183 Aprev = EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul;
00184 EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn =
00185 ratio*EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].opacity;
00186
00187 EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauCon = EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn;
00188 EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauTot =
00189 2.f*EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn;
00190
00191
00192 }
00193 ++ i;
00194 }
00195 }
00196 }
00197 }
00198
00199
00200 for( i=1; i <= nLevel1; i++ )
00201 {
00202 RT_line_one_tau_reset(&TauLines[i],WeightNew);
00203
00204 ASSERT( fabs(TauLines[i].TauIn) <= fabs(TauLines[i].TauTot) );
00205 }
00206
00207
00208 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine*sizeof(float) );
00209
00210
00211 for( i=0; i < nWindLine; i++ )
00212 {
00213 if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
00214 {
00215 RT_line_one_tau_reset(&TauLine2[i],WeightNew);
00216 }
00217 }
00218
00219
00220 for( i=0; i < nHFLines; i++ )
00221 {
00222 RT_line_one_tau_reset(&HFLines[i],WeightNew);
00223 }
00224
00225
00226 for( i=0; i < nUTA; i++ )
00227 {
00228 if( UTALines[i].Aul > 0. )
00229 {
00230
00231
00232
00233 double hsave = UTALines[i].heat;
00234 RT_line_one_tau_reset(&UTALines[i],WeightNew);
00235 UTALines[i].heat = hsave;
00236 }
00237 }
00238
00239
00240 for( i=0; i < nCORotate; i++ )
00241 {
00242 RT_line_one_tau_reset(&C12O16Rotate[i],WeightNew);
00243 RT_line_one_tau_reset(&C13O16Rotate[i],WeightNew);
00244 }
00245
00246
00247 H2_RT_tau_reset();
00248
00249
00250 FeII_RT_tau_reset();
00251
00252 if( opac.lgCaseB )
00253 {
00254 for( i=0; i < rfield.nupper; i++ )
00255 {
00256
00257
00258
00259 opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]/2.f;
00260
00261 opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]/2.f;
00262 }
00263 }
00264 else if( geometry.lgSphere )
00265 {
00266
00267 for( i=0; i < rfield.nupper; i++ )
00268 {
00269
00270
00271 opac.TauAbsGeo[1][i] = 2.f*opac.TauAbsFace[i];
00272 opac.TauAbsGeo[0][i] = opac.TauAbsFace[i];
00273 opac.TauScatGeo[1][i] = 2.f*opac.TauScatFace[i];
00274 opac.TauScatGeo[0][i] = opac.TauScatFace[i];
00275 opac.TauTotalGeo[1][i] = opac.TauScatGeo[1][i] + opac.TauAbsGeo[1][i];
00276 opac.TauTotalGeo[0][i] = opac.TauScatGeo[0][i] + opac.TauAbsGeo[0][i];
00277 }
00278 }
00279 else
00280 {
00281
00282 for( i=0; i < rfield.nupper; i++ )
00283 {
00284 opac.TauTotalGeo[1][i] = opac.TauTotalGeo[0][i];
00285 opac.TauTotalGeo[0][i] = opac.taumin;
00286 opac.TauAbsGeo[1][i] = opac.TauAbsGeo[0][i];
00287 opac.TauAbsGeo[0][i] = opac.taumin;
00288 opac.TauScatGeo[1][i] = opac.TauScatGeo[0][i];
00289 opac.TauScatGeo[0][i] = opac.taumin;
00290 }
00291 }
00292
00293
00294 for( i=0; i < rfield.nupper; i++ )
00295 {
00296
00297 opac.TauAbsTotal[i] = opac.TauAbsFace[i];
00298
00299
00300 opac.E2TauAbsTotal[i] = (float)e2( opac.TauAbsTotal[i] );
00301
00302 opac.TauScatFace[i] = opac.taumin;
00303 opac.TauAbsFace[i] = opac.taumin;
00304 }
00305
00306
00307 rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];
00308
00309 DEBUG_EXIT( "RT_tau_reset()" );
00310 return;
00311 }
00312