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

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_line_all do escape and destruction probabilities for all lines in code.  
00004  * Called with false arg by ionize, to only redo destruction probabilities,
00005  * and with true by cloudy to do both escape and destruction */
00006 #include "cddefines.h"
00007 #include "taulines.h"
00008 #include "atomfeii.h"
00009 #include "dense.h"
00010 #include "conv.h"
00011 #include "atoms.h"
00012 #include "rfield.h"
00013 #include "wind.h"
00014 #include "iso.h"
00015 #include "h2.h"
00016 #include "opacity.h"
00017 #include "trace.h"
00018 #include "lines_service.h"
00019 #include "atmdat.h"
00020 #include "hydrogenic.h"
00021 #include "rt.h"
00022 
00023 void RT_line_all(
00024         /* this is true if we want to do both escape and destruction probabilities,
00025          * and false if only destruction probabilities are needed */
00026         bool lgDoEsc ,
00027         /* flag saying whether to update fine opacities */
00028         bool lgUpdateFineOpac )
00029 {
00030         long int i,
00031                 ion,
00032                 ipISO,
00033                 nelem;
00034         long ipHi , ipLo;
00035         double factor,
00036                 coloi;
00037         double SaveLyaPesc[NISO][LIMELM] , 
00038                 SaveLyaPdest[NISO][LIMELM];
00039         /* this flag says whether to update the line RT */
00040         bool lgRT_update = lgUpdateFineOpac || lgDoEsc || conv.lgIonStageTrimed;
00041 
00042         DEBUG_ENTRY( "RT_line_all()" );
00043 
00044         if( trace.lgTrace )
00045                 fprintf( ioQQQ, "     RT_line_all called\n" );
00046 
00047         /* skip line transfer if requested, and not first call in this zone */
00048         if( !rfield.lgDoLineTrans && conv.nPres2Ioniz )
00049         {
00050                 return;
00051         }
00052 
00053         /* this array is huge and takes significant time to zero out or update, 
00054          * only do so when needed, */
00055         rfield.lgFine_opac_update = lgUpdateFineOpac;
00056         if( rfield.lgFine_opac_update )
00057         {
00058                 /* zero out fine opacity array */
00059                 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
00060 
00061                 /* this is offset within fine opacity array caused by Doppler shift
00062                  * between first zone, the standard of rest, and current position 
00063                  * in case of acceleration away from star in outflowing wind 
00064                  * dWind is positive, this means that the rest frame original
00065                  * velocity is red shifted to lower energy */
00066                 realnum dWind = wind.windv - wind.windv0;
00067                 /* will add ipVelShift to index of original energy so redshift has 
00068                  * to be negative */
00069                 rfield.ipFineConVelShift = -(long int)(dWind / rfield.fine_opac_velocity_width+0.5);
00070         }
00071 
00072 #       if 0
00073         /* this code is a copy of the code within line_one that does cloaking
00074          * within this zone.  it can be used to see how a particular line
00075          * is being treated. */
00076         {
00077 #include "doppvel.h"
00078                 double dTau , aa;
00079                 ipISO = 0; nelem = 0;ipLo = 0;
00080                 ipHi = 23;
00081                 dTau =  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc * 
00082                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity / 
00083                         DoppVel.doppler[nelem] + opac.opacity_abs[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
00084                 dTau *= radius.drad_x_fillfac_mean;
00085                 aa = log(1. + dTau ) / SDIV(dTau);
00086                 fprintf(ioQQQ,"DEBUG dTau\t%.2f\t%.5e\t%.5e\t%.5e\n",fnzone,dTau,
00087                         radius.drad_x_fillfac_mean,
00088                          aa);
00089         }
00090 #       endif
00091 
00092         /* find Stark escape probabilities for hydrogen itself */
00093         RT_stark();
00094 
00095         /*if( lgUpdateFineOpac )
00096                 fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac in rt_line_all\n");*/
00097         for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00098         {
00099                 /* loop over all iso-electronic sequences */
00100                 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00101                 {
00102                         /* parent ion stage, for H is 1, for He is 1 for He-like and 
00103                          * 2 for H-like */
00104                         ion = nelem+1-ipISO;
00105 
00106                         /* element turned off */
00107                         if( !dense.lgElmtOn[nelem] )
00108                                 continue;
00109                         /* need we consider this ion? */
00110                         if( ion <=dense.IonHigh[nelem] )
00111                         {
00112                                 /* save escape and destruction probs for each Lya since
00113                                  * these are special cases, and quantities are not damped
00114                                  * in routine that evaluates each, test is to avoid
00115                                  * evaluating iso sequence that does not exist for an element,
00116                                  * as in He-like hydrogen */
00117                                 if( ipISO<=nelem )
00118                                 {
00119                                         SaveLyaPesc[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc;
00120                                         SaveLyaPdest[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest;
00121                                 }
00122 
00123                                 /* convert pops to per unit vol rather than per ion */
00124                                 if( dense.xIonDense[nelem][ion] > 1e-30 )
00125                                 {
00126                                         factor = dense.xIonDense[nelem][ion];
00127                                 }
00128                                 else
00129                                 {
00130                                         /* case where almost no parent ion - this will make
00131                                          * very large line opacity, so background dest small */
00132                                         factor = 1.;
00133                                 }
00134 
00135                                 /* loop over all lines */
00136                                 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ++ipHi )
00137                                 {
00138                                         for( ipLo=0; ipLo < ipHi; ++ipLo )
00139                                         {
00140                                                 /* negative ipCont means this is not a real line, so do not
00141                                                  * transfer it */
00142                                                 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 0 ) 
00143                                                         continue;
00144 
00145                                                 double SavePopOpc = Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc;
00146                                                 /* must temporarily make ipLnPopOpc physical */
00147                                                 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc *= factor;
00148                                                 /*if( ipISO==0 && nelem==0 && ipHi==12 && ipLo==11 )
00149                                                         fprintf(ioQQQ,"DEBUG rt loop %li %li %li %li %.3e %.3e %.3e\n",
00150                                                         ipISO, nelem, ipHi , ipLo , 
00151                                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc , 
00152                                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn,
00153                                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot);*/
00154 
00155                                                 /* generate escape prob, pumping rate, destruction prob, 
00156                                                  * inward outward fracs */
00157                                                 RT_line_one(&Transitions[ipISO][nelem][ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,true);
00158 
00159 
00160                                                 {
00161                                                         /* set true to print pump rates*/
00162                                                         enum {DEBUG_LOC=false};
00163                                                         if( DEBUG_LOC && nelem==1&& ipLo==0  /*&& iteration==2*/ )
00164                                                         {
00165                                                                 fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n",
00166                                                                         ipHi ,
00167                                                                         fnzone,
00168                                                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest);
00169                                                         }
00170                                                 }
00171 
00172                                                 /* go back to original units so that final correction ok */
00173                                                 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = SavePopOpc;
00174                                         }
00175                                 }
00176 
00177                                 if( ipISO > 0 && lgDoEsc )
00178                                 {
00179                                         /* do not let too much Lya escape outward since so important 
00180                                          * H Lya is special case done above */
00181                                         Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc *= 
00182                                                 opac.ExpmTau[Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].ipCont-1];
00183                                 }
00184 
00185                                 /* now update two photon induced rates */
00186                                 atmdat_2phot_rate(ipISO , nelem);
00187 
00188                                 /* this is option to not do destruction probabilities
00189                                  * for case b no pdest option */
00190                                 if( opac.lgCaseB_no_pdest )
00191                                 {
00192                                         ipLo = 0;
00193                                         /* okay to let this go to numLevels_max. */
00194                                         for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem]; ++ipHi )
00195                                         {
00196                                                 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
00197                                                         continue;
00198 
00199                                                 /* hose the previously computed destruction probability */
00200                                                 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest = SMALLFLOAT; 
00201                                         }
00202                                 }
00203 
00204                                 ipLo = 0;
00205                                 /* these are the extra lyman lines for the iso sequences */
00206                                 /*for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )*/
00207                                 /* only update if significant abundance and need to update fine opac */
00208                                 if( dense.xIonDense[nelem][ion] > 1e-30 && lgUpdateFineOpac )
00209                                 {
00210                                         /* we just want the population of the ground state */
00211                                         factor = StatesElem[ipISO][nelem][0].Pop * dense.xIonDense[nelem][ion];
00212                                         for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ )
00213                                         {
00214                                                 /* must make ipLnPopOpc physical */
00215                                                 ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc = factor;
00216                                                 ExtraLymanLines[ipISO][nelem][ipHi].Lo->Pop =
00217                                                         StatesElem[ipISO][nelem][ipLo].Pop;
00218 
00219                                                 /* actually do the work */
00220                                                 RT_line_one(&ExtraLymanLines[ipISO][nelem][ipHi] , lgDoEsc , lgUpdateFineOpac,true);
00221 
00222                                                 /* reset to funny population of ground state */
00223                                                 ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc = 
00224                                                         StatesElem[ipISO][nelem][0].Pop;
00225                                         }
00226                                 }
00227 
00228                         }/* if nelem if ion <=dense.IonHigh */
00229                 }/* loop over nelem */
00230         }/* loop over ipISO */
00231 
00232         
00233         /* add on Stark escape probabilities for H itself */
00234         for( ipISO=0; ipISO<NISO; ipISO++ )
00235         {
00236                 for( nelem=ipISO; nelem<LIMELM; nelem++ )
00237                 {
00238                         if( nelem < 2 || dense.lgElmtOn[nelem] )
00239                         {
00240                                 /* note that pesc and dest are updated no matter what escprob logic we
00241                                  * specify, and that these are not updated if we have overrun the
00242                                  * optical depth scale.  only update here in that case */
00243                                 if( lgTauGood( &Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0] ) )
00244                                 {
00245                                         if( lgDoEsc )
00246                                         {
00247                                                 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00248                                                 {
00249                                                         for( ipLo=0; ipLo < ipHi; ipLo++ )
00250                                                         {
00251                                                                 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
00252                                                                         continue;
00253 
00254                                                                 /* >>chng 03 jun 07, do not clobber esp prob when line is masing -
00255                                                                  * this had effect of preventing total escape prob from getting larger than 1 */
00256                                                                 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc<1. )
00257                                                                 {
00258                                                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc = MIN2(1.f,
00259                                                                                 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc+
00260                                                                                 (realnum)iso.pestrk[ipISO][nelem][ipLo][ipHi]);
00261                                                                 }
00262                                                         }
00263                                                 }
00264                                         }
00265                                         /* always add on stark broadening term to Lya*/
00266                                         else if( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc < 1. )
00267                                         {
00268                                                 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc = MIN2(1.f,
00269                                                         Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc+
00270                                                         (realnum)iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]]);
00271                                         }
00272                                 }
00273                         }
00274                 }
00275         }
00276 
00277         /* note that pesc and dest are updated no matter what escprob logic we
00278          * specify, and that these are not updated if we have overrun the
00279          * optical depth scale.  only update here in that case */
00280         if( lgTauGood( &Transitions[ipH_LIKE][ipHYDROGEN][iso.nLyaLevel[ipH_LIKE]][0] ) )
00281         {
00282                 /* do the 8446 problem */
00283                 atom_oi_calc(&coloi);
00284 
00285                 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc = (realnum)(atoms.pmph31/
00286                         Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul);
00287 
00288                 if( trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE] )
00289                 {
00290                         fprintf( ioQQQ, "       RT_line_all calls P8446 who found pmph31=%10.2e\n", 
00291                                 atoms.pmph31 );
00292                 }
00293 
00294                 /* add on destruction of hydrogen Lya by FeII
00295                 * now add in FeII deexcitation via overlap,
00296                 * but only as loss, not photoionization, source
00297                 * dstfe2lya is Ly-alpha deexcit by FeII overlap - conversion into FeII em */
00298                 /* find FeII overlap destruction rate, 
00299                 * this does NOTHING when large FeII atom is turned on, 
00300                 * in this case evaluation is done in call to FeIILevelPops */
00301                 t_fe2ovr_la::Inst().atoms_fe2ovr();
00302                 /*fprintf(ioQQQ,"DEBUG fe2 %.2e %.2e\n", hydro.dstfe2lya ,
00303                         hydro.HLineWidth);*/
00304 
00305                 /* >>chng 00 jan 06, let dest be large than 1 to desaturate the atom */
00306                 /* >>chng 01 apr 01, add test for tout >= 0., 
00307                  * must not add to Pdest when it has not been refreshed here */
00308                 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest += hydro.dstfe2lya;
00309         }
00310 
00311         /* take mean of old and new to damp oscillations in Lya (only) */
00312         if( nzone > 1 )
00313         {
00314                 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00315                 {
00316                         /* loop over all iso-electronic sequences */
00317                         for( nelem=ipISO; nelem<LIMELM; ++nelem )
00318                         {
00319                                 ion = nelem+1-ipISO;
00320                                 if( ion <=dense.IonHigh[nelem] && ipISO<=nelem )
00321                                 {
00322                                         /* now take mean of old and new escape dest probs for Lya */
00323                                         /* only ipLY_A redist did not already take average */
00324                                         if( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->iRedisFun==ipLY_A &&
00325                                                 lgTauGood( &Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0] ) )
00326                                         {
00327                                                 /* the Lya routine is different in that escape and dest
00328                                                  * probs are always updated, unless we have overrun the
00329                                                  * optical depth scale */
00330                                                 /*lint -e644 SaveLyaPdest not init */
00331                                                 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest = 
00332                                                         ((realnum)SaveLyaPdest[ipISO][nelem] + 
00333                                                         Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest) / 2.f;
00334                                                 Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc = 
00335                                                         ((realnum)SaveLyaPesc[ipISO][nelem] + 
00336                                                         Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc) / 2.f;
00337                                                 /*lint +e644 SaveLyaPdest not init */
00338                                         }
00339                                 }
00340                         }
00341                 }
00342         }
00343         /*5986 fprintf(ioQQQ,"DEBUG he rt_all\t%li\t%li\t%.4e\t%.4e\n",
00344                 iteration, nzone, 
00345                                 Transitions[1][1][10][5].Emis->TauIn,
00346                                 Transitions[1][1][10][5].Emis->Pesc);*/
00347 
00348         /* is continuum pumping of H Lyman lines included?  yes, but turned off
00349          * with atom h-like Lyman pumping off command */
00350         if( !hydro.lgLymanPumping )
00351         {
00352                 ipISO = ipH_LIKE;
00353                 nelem = ipHYDROGEN;
00354                 ipLo = 0;
00355                 for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem]; ++ipHi )
00356                 {
00357                         /* negative ipCont means this is not a real line, so do not
00358                                 * transfer it */
00359                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->pump = 0.;
00360                 }
00361         }
00362 
00363         {
00364                 /* following should be set true to print ots contributors for he-like lines*/
00365                 /*@-redef@*/
00366                 enum {DEBUG_LOC=false};
00367                 /*@+redef@*/
00368                 if( DEBUG_LOC && nzone>433 /*&& iteration==2*/ )
00369                 {
00370                         /* option to dump a line  */
00371                         DumpLine(&Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0] );
00372                 }
00373         }
00374 
00375         /* level 1 lines */
00376         for( i=1; i <= nLevel1; i++ )
00377         {
00378                 RT_line_one(&TauLines[i] , lgDoEsc , lgUpdateFineOpac,true);
00379         }
00380         /* co carbon monoxide */
00381         for( i=0; i < nCORotate; i++ )
00382         {
00383                 RT_line_one(&C12O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,true);
00384                 RT_line_one(&C13O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,true);
00385         }
00386         for( i=0; i < nHFLines; i++ )
00387         {
00388                 RT_line_one(&HFLines[i] , lgDoEsc , lgUpdateFineOpac,true);
00389         }
00390         /* Database Atoms & Molecules*/
00391         for( i=0; i < linesAdded2; i++)
00392         {
00393                 RT_line_one(atmolEmis[i].tran , lgDoEsc , lgUpdateFineOpac,true);
00394         }
00395         /* this is a major time sink for this routine */
00396         if( lgRT_update )
00397         {
00398                 for( i=0; i < nUTA; i++ )
00399                 {
00400                         /* Aul = -1 means this line is superceded by better calculations in
00401                          * different data set */
00402                         if( UTALines[i].Emis->Aul > 0. )
00403                         {
00404                                 /* these are not defined in cooling routines so must be set here */
00405                                 UTALines[i].Emis->PopOpc = dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1];
00406                                 UTALines[i].Lo->Pop = dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1];
00407                                 UTALines[i].Hi->Pop = 0.;
00408                                 RT_line_one(&UTALines[i] , lgDoEsc , lgUpdateFineOpac,true);
00409                         }
00410                 }
00411         }
00412 
00413         /* h-like ions only have one electron, no satellite lines. */
00414         for( ipISO=ipHE_LIKE; ipISO<NISO; ++ipISO )
00415         {
00416                 /* loop over all iso-electronic sequences */
00417                 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00418                 {
00419                         if( dense.lgElmtOn[nelem] )
00420                         {
00421                                 for( ipLo=1; ipLo < (iso.numLevels_local[ipISO][nelem]-1); ++ipLo )
00422                                 {
00423                                         RT_line_one(&SatelliteLines[ipISO][nelem][ipLo], lgDoEsc, lgUpdateFineOpac, true);
00424                                 }
00425                         }
00426                 }
00427         }
00428 
00429         /* >>chng 03 aug 22, must call RT_line_one every single time, since that routine
00430          * establishes the fine opacities for the level 2 lines */
00431         /* >>chng 01 aug 11, the level 2 lines were only updated when lgDoEsc was true,
00432          * this meant that dest probs were not updated but once, causing very high-Z models
00433          * to develop numerical oscialltions.  removed if, now evaluate level 2 always */
00434         /* only update their dest/esc probs one time per zone */
00435         /* lgLevel2_OTS_Imp set true in dimacool if ots rates were significant */
00436         /* level 2 heavy element lines in cooling with only g-bar,*/
00437         for( i=0; i < nWindLine; i++ )
00438         {
00439                 /* >>chng 02 sug 11, do not include ions done in iso-sequences */
00440                 /*if( TauLine2[i].IonStg <= TauLine2[i].nelem-1 )*/
00441                 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
00442                 {
00443                         RT_line_one(&TauLine2[i] , lgDoEsc , lgUpdateFineOpac,true);
00444                 }
00445         }
00446 
00447         /* the large H2 molecule */
00448         H2_RTMake( lgDoEsc , lgUpdateFineOpac);
00449 
00450         /* The large model FeII atom */
00451         FeII_RT_Make( lgDoEsc , lgUpdateFineOpac);
00452         return;
00453 }
/home66/gary/public_html/cloudy/c08_branch/source/rt_line_all.cpp
