cloudy: /home66/gary/public_html/cloudy/c08_branch/source/rt_tau

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 }
/home66/gary/public_html/cloudy/c08_branch/source/rt_tau_reset.cpp
