cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*RT_tau_reset after first iteration, updates the optical depths, mirroring this 00004 * routine but with the previous iteration's variables */ 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 /*RT_tau_reset update total optical depth scale, 00021 * called by cloudy after iteration is complete */ 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 if( Transitions[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].ipCont <= 0 ) 00046 fprintf( ioQQQ, "%10.2e", 1e-30 ); 00047 else 00048 fprintf( ioQQQ, "%10.2e", 00049 Transitions[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].Emis->TauIn ); 00050 } 00051 fprintf( ioQQQ, "\n" ); 00052 } 00053 } 00054 } 00055 00056 /* pumping of CaH */ 00057 opac.tpcah[1] = opac.tpcah[0]; 00058 opac.tpcah[0] = opac.taumin; 00059 00060 /* save column densities of H species */ 00061 for( i=0; i<NCOLD; ++i ) 00062 { 00063 colden.colden_old[i] = colden.colden[i]; 00064 } 00065 for( i=0; i < mole.num_comole_calc; i++ ) 00066 { 00067 COmole[i]->hevcol_old = COmole[i]->hevcol; 00068 } 00069 00070 /* ======================================================================== */ 00071 /* must take average of old and new optical depths - were old in place? */ 00072 if( iteration <= 1 ) 00073 { 00074 /* this is first pass */ 00075 WeightNew = 1.; 00076 } 00077 else 00078 { 00079 WeightNew = 0.75; 00080 } 00081 00082 opac.telec = opac.taumin; 00083 opac.thmin = opac.taumin; 00084 00085 /* all iso sequences */ 00086 for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00087 { 00088 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00089 { 00090 if( dense.lgElmtOn[nelem] ) 00091 { 00092 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00093 { 00094 if( iso.lgDielRecom[ipISO] ) 00095 RT_line_one_tau_reset(&SatelliteLines[ipISO][nelem][ipHi],WeightNew); 00096 00097 /* the bound-bound transitions within the model atoms */ 00098 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00099 { 00100 enum {DEBUG_LOC=false}; 00101 if( DEBUG_LOC ) 00102 { 00103 if( ipISO==1 && nelem==1 && ipHi==iso.nLyaLevel[ipISO] && ipLo==0 ) 00104 fprintf(ioQQQ,"DEBUG rt before loop %li %li %li %li tot %.3e in %.3e\n", 00105 ipISO, nelem, ipHi , ipLo , 00106 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot , 00107 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn); 00108 } 00109 /*RT_line_one_tau_reset computes average of old and new optical depths 00110 * for new scale at end of iter */ 00111 RT_line_one_tau_reset(&Transitions[ipISO][nelem][ipHi][ipLo],WeightNew); 00112 if( DEBUG_LOC ) 00113 { 00114 if( ipISO==1 && nelem==1 && ipHi==iso.nLyaLevel[ipISO] && ipLo==0 ) 00115 fprintf(ioQQQ,"DEBUG rt after loop %li %li %li %li tot %.3e in %.3e\n", 00116 ipISO, nelem, ipHi , ipLo , 00117 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot , 00118 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn); 00119 } 00120 } 00121 } 00122 /* the extra Lyman lines */ 00123 for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ ) 00124 { 00125 /* fully transfer all of the extra lines even though 00126 * have not solved for their upper level populations */ 00127 RT_line_one_tau_reset(&ExtraLymanLines[ipISO][nelem][ipHi],WeightNew); 00128 } 00129 } 00130 } 00131 } 00132 00133 /* >>>chng 99 nov 11 did not have case b for hydrogenic species on second and 00134 * higher iterations */ 00135 /* option to clobber these taus for Lyman lines, if case b is set */ 00136 if( opac.lgCaseB ) 00137 { 00138 for( nelem=0; nelem < LIMELM; nelem++ ) 00139 { 00140 if( dense.lgElmtOn[nelem] ) 00141 { 00142 realnum f; 00143 /* La may be case B, tlamin set to 1e9 by default with case b command */ 00144 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn = opac.tlamin; 00145 /* >>>chng 99 nov 22, did not reset TauCon */ 00146 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauCon = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn; 00147 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = 00148 2.f*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn; 00149 f = opac.tlamin/Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity; 00150 00151 for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00152 { 00153 if( Transitions[ipH_LIKE][nelem][ipHi][ipH1s].ipCont <= 0 ) 00154 continue; 00155 00156 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = 00157 f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity; 00158 /* reset line optical depth to continuum source */ 00159 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauCon = Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn; 00160 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot = 00161 2.f*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn; 00162 } 00163 } 00164 } 00165 00166 /* now do helium like sequence - different since collapsed levels 00167 * all go to ground */ 00168 for( nelem=1; nelem < LIMELM; nelem++ ) 00169 { 00170 if( dense.lgElmtOn[nelem] ) 00171 { 00172 double Aprev; 00173 realnum ratio; 00174 /* La may be case B, tlamin set to 1e9 by default with case b command */ 00175 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn = opac.tlamin; 00176 00177 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauCon = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn; 00178 00179 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauTot = 00180 2.f*Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn; 00181 00182 ratio = opac.tlamin/Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->opacity; 00183 00184 /* this will be the trans prob of the previous lyman line, will use this to 00185 * find the next one up in the series */ 00186 Aprev = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->Aul; 00187 00188 i = ipHe2p1P+1; 00189 /* >>chng 02 jan 05, remove explicit list of lyman lines, use As to guess 00190 * which are which - this will work for any number of levels */ 00191 for( i = ipHe2p1P+1; i < iso.numLevels_max[ipHE_LIKE][nelem]; i++ ) 00192 { 00193 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 ) 00194 continue; 00195 00196 /* >>chng 02 mar 19 use proper test for resonance collapsed lines */ 00197 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul> Aprev/10. || 00198 StatesElem[ipHE_LIKE][nelem][i].S < 0 ) 00199 { 00200 Aprev = Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul; 00201 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn = 00202 ratio*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->opacity; 00203 /* reset line optical depth to continuum source */ 00204 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauCon = 00205 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn; 00206 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauTot = 00207 2.f*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn; 00208 /*fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",nelem, i, 00209 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Aul, Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn );*/ 00210 } 00211 } 00212 } 00213 } 00214 } 00215 00216 /* all heavy element lines */ 00217 for( i=1; i <= nLevel1; i++ ) 00218 { 00219 RT_line_one_tau_reset(&TauLines[i],WeightNew); 00220 /* >>chng 96 jul 06 following sanity check added */ 00221 ASSERT( fabs(TauLines[i].Emis->TauIn) <= fabs(TauLines[i].Emis->TauTot) ); 00222 } 00223 00224 /* zero out fine opacity fine grid fine mesh array */ 00225 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine*sizeof(realnum) ); 00226 00227 /* all level 2 heavy element lines */ 00228 for( i=0; i < nWindLine; i++ ) 00229 { 00230 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00231 { 00232 RT_line_one_tau_reset(&TauLine2[i],WeightNew); 00233 } 00234 } 00235 00236 /* all hyperfine structure lines */ 00237 for( i=0; i < nHFLines; i++ ) 00238 { 00239 RT_line_one_tau_reset(&HFLines[i],WeightNew); 00240 } 00241 00242 /*Atoms & Molecules*/ 00243 for( i=0; i < linesAdded2; i++) 00244 { 00245 RT_line_one_tau_reset( atmolEmis[i].tran,WeightNew ); 00246 } 00247 00248 /* inner shell lines */ 00249 for( i=0; i < nUTA; i++ ) 00250 { 00251 if( UTALines[i].Emis->Aul > 0. ) 00252 { 00253 /* these are line optical depth arrays 00254 * inward optical depth */ 00255 /* heat is special for this array - it is heat per pump */ 00256 double hsave = UTALines[i].Coll.heat; 00257 RT_line_one_tau_reset(&UTALines[i],WeightNew); 00258 UTALines[i].Coll.heat = hsave; 00259 } 00260 } 00261 00262 /* co carbon monoxide lines */ 00263 for( i=0; i < nCORotate; i++ ) 00264 { 00265 RT_line_one_tau_reset(&C12O16Rotate[i],WeightNew); 00266 RT_line_one_tau_reset(&C13O16Rotate[i],WeightNew); 00267 } 00268 00269 /* the large H2 molecule */ 00270 H2_RT_tau_reset(); 00271 00272 /* large FeII atom */ 00273 FeII_RT_tau_reset(); 00274 00275 if( opac.lgCaseB ) 00276 { 00277 for( i=0; i < rfield.nupper; i++ ) 00278 { 00279 /* DEPABS and SCT are abs and sct optical depth for depth only 00280 * we will not change total optical depths, just reset inner to half 00281 * TauAbsGeo(i,2) = 2.*TauAbsFace(i) */ 00282 opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]/2.f; 00283 /* TauScatGeo(i,2) = 2.*TauScatFace(i) */ 00284 opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]/2.f; 00285 } 00286 } 00287 else if( geometry.lgSphere ) 00288 { 00289 /* closed or spherical case */ 00290 for( i=0; i < rfield.nupper; i++ ) 00291 { 00292 /* [1] is total optical depth from previous iteration, 00293 * [0] is optical depth at current position */ 00294 opac.TauAbsGeo[1][i] = 2.f*opac.TauAbsFace[i]; 00295 opac.TauAbsGeo[0][i] = opac.TauAbsFace[i]; 00296 opac.TauScatGeo[1][i] = 2.f*opac.TauScatFace[i]; 00297 opac.TauScatGeo[0][i] = opac.TauScatFace[i]; 00298 opac.TauTotalGeo[1][i] = opac.TauScatGeo[1][i] + opac.TauAbsGeo[1][i]; 00299 opac.TauTotalGeo[0][i] = opac.TauScatGeo[0][i] + opac.TauAbsGeo[0][i]; 00300 } 00301 } 00302 else 00303 { 00304 /* open geometry */ 00305 for( i=0; i < rfield.nupper; i++ ) 00306 { 00307 opac.TauTotalGeo[1][i] = opac.TauTotalGeo[0][i]; 00308 opac.TauTotalGeo[0][i] = opac.taumin; 00309 opac.TauAbsGeo[1][i] = opac.TauAbsGeo[0][i]; 00310 opac.TauAbsGeo[0][i] = opac.taumin; 00311 opac.TauScatGeo[1][i] = opac.TauScatGeo[0][i]; 00312 opac.TauScatGeo[0][i] = opac.taumin; 00313 } 00314 } 00315 00316 /* same for open and closed geometries */ 00317 for( i=0; i < rfield.nupper; i++ ) 00318 { 00319 /* total optical depth across computed shell */ 00320 opac.TauAbsTotal[i] = opac.TauAbsFace[i]; 00321 /* e2( tau across shell), optical depth from ill face to shielded face of cloud 00322 * not that opac.TauAbsFace is reset to small number just after this */ 00323 opac.E2TauAbsTotal[i] = (realnum)e2( opac.TauAbsTotal[i] ); 00324 /* TauAbsFace and TauScatFace are abs and sct optical depth to ill face */ 00325 opac.TauScatFace[i] = opac.taumin; 00326 opac.TauAbsFace[i] = opac.taumin; 00327 } 00328 00329 /* this is optical depth at x-ray point defining effective optical depth */ 00330 rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1]; 00331 return; 00332 }