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_inc increment optical depths once per zone, called after radius_increment */ 00004 #include "cddefines.h" 00005 #include "taulines.h" 00006 #include "iso.h" 00007 #include "rfield.h" 00008 #include "trace.h" 00009 #include "dense.h" 00010 #include "hyperfine.h" 00011 #include "wind.h" 00012 #include "prt.h" 00013 #include "h2.h" 00014 #include "hmi.h" 00015 #include "opacity.h" 00016 #include "radius.h" 00017 #include "atomfeii.h" 00018 #include "rt.h" 00019 00020 /*RT_tau_inc increment optical depths once per zone, called after radius_increment */ 00021 void RT_tau_inc(void) 00022 { 00023 00024 long int i, 00025 ipHi, 00026 nelem, 00027 ipLo, 00028 ipISO; 00029 00030 double factor; 00031 00032 DEBUG_ENTRY( "RT_tau_inc()" ); 00033 00034 if( trace.lgTrace ) 00035 { 00036 fprintf( ioQQQ, " RT_tau_inc called.\n" ); 00037 } 00038 00039 /* call RT_line_all one last time in this zone, to get fine opacities defined */ 00040 RT_line_all( false , true); 00041 00042 opac.telec += (realnum)(radius.drad_x_fillfac*dense.eden*6.65e-25); 00043 opac.thmin += (realnum)(radius.drad_x_fillfac*hmi.Hmolec[ipMHm]*3.9e-17* 00044 (1. - rfield.ContBoltz[hmi.iphmin-1]/ hmi.hmidep)); 00045 00046 /* prevent maser runaway */ 00047 rt.dTauMase = 0; 00048 rt.mas_species = 0; 00049 rt.mas_ion = 0; 00050 rt.mas_hi = 0; 00051 rt.mas_lo = 0; 00052 00053 /* all lines in iso sequences */ 00054 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00055 { 00056 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00057 { 00058 /* this is the parent ion, for HI lines, is 1, 00059 * for element He is 1 for He-like (HeI) and 2 for H-like (HeII) */ 00060 int ion = nelem+1-ipISO; 00061 /* do not evaluate in case where trivial parent ion */ 00062 if( ion <=dense.IonHigh[nelem] && dense.xIonDense[nelem][ion] > dense.density_low_limit ) 00063 { 00064 factor = dense.xIonDense[nelem][ion]; 00065 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00066 { 00067 if( iso.lgDielRecom[ipISO] ) 00068 RT_line_one_tauinc(&SatelliteLines[ipISO][nelem][ipHi], 00069 ipISO , nelem , -1 , ipHi ); 00070 00071 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00072 { 00073 /* must temporarily make ipLnPopOpc physical */ 00074 double save; 00075 00076 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 ) 00077 continue; 00078 00079 save = Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc; 00080 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc *= factor; 00081 /* actually do the work */ 00082 RT_line_one_tauinc(&Transitions[ipISO][nelem][ipHi][ipLo] , 00083 ipISO , nelem , ipHi , ipLo ); 00084 /* go back to original units so that final correction ok */ 00085 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = save; 00086 } 00087 } 00088 ipLo = 0; 00089 /* these are the extra Lyman lines */ 00090 for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ ) 00091 { 00092 /* must make ipLnPopOpc physical */ 00093 ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc = 00094 StatesElem[ipISO][nelem][0].Pop * factor; 00095 00096 /* actually do the work */ 00097 RT_line_one_tauinc(&ExtraLymanLines[ipISO][nelem][ipHi] , 00098 -1 ,ipISO , nelem , ipHi ); 00099 00100 /* reset to population of ground state */ 00101 ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc = 00102 StatesElem[ipISO][nelem][0].Pop; 00103 } 00104 } 00105 } 00106 } 00107 00108 /* increment optical depths for all heavy element lines 00109 * same routine does wind and static, 00110 * does not start from 0 since first line is dummy */ 00111 for( i=1; i <= nLevel1; i++ ) 00112 { 00113 RT_line_one_tauinc(&TauLines[i], 00114 -2 , -2 , -2 , i ); 00115 } 00116 00117 /* all lines in cooling with g-bar */ 00118 for( i=0; i < nWindLine; i++ ) 00119 { 00120 /* do not include H-like or He-like in the level two lines since 00121 * these are already counted in iso sequences */ 00122 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00123 { 00124 RT_line_one_tauinc(&TauLine2[i] , 00125 -3 , -3 , -3 , i ); 00126 } 00127 } 00128 00129 /* the block of inner shell lines */ 00130 for( i=0; i < nUTA; i++ ) 00131 { 00132 if( UTALines[i].Emis->Aul > 0. ) 00133 { 00134 /* populations have not been set */ 00135 UTALines[i].Emis->PopOpc = dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1]; 00136 UTALines[i].Lo->Pop = dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1]; 00137 UTALines[i].Hi->Pop = 0.; 00138 RT_line_one_tauinc(&UTALines[i], -4 , -4 , -4 , i ); 00139 } 00140 } 00141 00142 /* all hyper fine structure lines */ 00143 for( i=0; i < nHFLines; i++ ) 00144 { 00145 /* remember current gas-phase abundances */ 00146 realnum save = dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1]; 00147 00148 /* bail if no abundance */ 00149 if( save<=0. ) continue; 00150 00151 /* set gas-phase abundance to total times isotope ratio */ 00152 dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1] *= hyperfine.HFLabundance[i]; 00153 00154 RT_line_one_tauinc(&HFLines[i] , -5 , -5 , -5 , i ); 00155 00156 /* put the correct gas-phase abundance back in the array */ 00157 dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1] = save; 00158 } 00159 00160 /* carbon monoxide CO lines */ 00161 for( i=0; i < nCORotate; i++ ) 00162 { 00163 RT_line_one_tauinc(&C12O16Rotate[i], 00164 -6 , -6 , -6 , i ); 00165 RT_line_one_tauinc(&C13O16Rotate[i], 00166 -7 , -7 , -7 , i ); 00167 } 00168 00169 /* do large FeII atom if this is enabled */ 00170 FeII_RT_TauInc(); 00171 00172 /* increment optical depth for the H2 molecule */ 00173 H2_RT_tau_inc(); 00174 00175 /* database Lines*/ 00176 for(i=0; i < linesAdded2; i++) 00177 { 00178 RT_line_one_tauinc(atmolEmis[i].tran, 00179 -10 , -10 , -10 , i ); 00180 } 00181 00182 /* following is for static atmosphere */ 00183 if( wind.windv == 0. ) 00184 { 00185 /* iron fe feii fe2 - overlapping feii lines */ 00186 t_fe2ovr_la::Inst().tau_inc(); 00187 } 00188 00189 if( trace.lgTrace && trace.lgOptcBug ) 00190 { 00191 fprintf( ioQQQ, " RT_tau_inc updated optical depths:\n" ); 00192 prtmet(); 00193 } 00194 00195 if( trace.lgTrace ) 00196 fprintf( ioQQQ, " RT_tau_inc returns.\n" ); 00197 00198 return; 00199 }