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_stark compute stark broadening escape probabilities using Puetter formalism */ 00004 #include "cddefines.h" 00005 #include "taulines.h" 00006 #include "iso.h" 00007 #include "dense.h" 00008 #include "hydrogenic.h" 00009 #include "phycon.h" 00010 #include "rt.h" 00011 00012 void RT_stark(void) 00013 { 00014 long int ipLo, 00015 ipHi, 00016 nelem, 00017 ipISO; 00018 00019 double aa , ah, 00020 stark, 00021 strkla; 00022 00023 DEBUG_ENTRY( "RT_stark()" ); 00024 00025 /* only evaluate one time per zone */ 00026 static long int nZoneEval=-1; 00027 if( nzone==nZoneEval ) 00028 return; 00029 nZoneEval = nzone; 00030 00031 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00032 { 00033 /* loop over all iso-electronic sequences */ 00034 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00035 { 00036 if( nelem >= 2 && !dense.lgElmtOn[nelem] ) 00037 continue; 00038 00039 if( !rt.lgStarkON || dense.eden < 1e8 ) 00040 { 00041 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00042 { 00043 for( ipLo=0; ipLo < iso.numLevels_max[ipISO][nelem]; ipLo++ ) 00044 { 00045 iso.pestrk[ipISO][nelem][ipHi][ipLo] = 0.; 00046 iso.pestrk[ipISO][nelem][ipLo][ipHi] = 0.; 00047 } 00048 } 00049 continue; 00050 } 00051 00052 /* evaluate Stark escape probability from 00053 * >>ref Puetter Ap.J. 251, 446. */ 00054 00055 /* coefficients for Stark broadening escape probability 00056 * to be Puetters AH, equation 9b, needs factor of (Z^-4.5 * (nu*nl)^3 * xl) */ 00057 ah = 6.9e-6*1000./1e12/(phycon.sqrte*phycon.te10*phycon.te10* 00058 phycon.te03*phycon.te01*phycon.te01)*dense.eden; 00059 00060 /* include Z factor */ 00061 ah *= pow( (double)(nelem+1), -4.5 ); 00062 00063 /* coefficient for all lines except Ly alpha */ 00064 /* equation 10b, except missing tau^-0.6 */ 00065 stark = 0.264*pow(ah,0.4); 00066 00067 /* coefficient for Ly alpha */ 00068 /* first few factors resemble equation 13c...what about the rest? */ 00069 strkla = 0.538*ah*4.*9.875*(phycon.sqrte/phycon.te10/phycon.te03); 00070 00071 /* Lyman lines always have outer optical depths */ 00072 /*ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn> 0. );*/ 00073 /* >>chng 02 mar 31, put in max, crashed on some first iteration 00074 * with negative optical depths, 00075 * NB did not understand why neg optical depths started */ 00076 aa = (realnum)SDIV(Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn); 00077 aa = pow( aa ,-0.75); 00078 iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]] = strkla/2.*MAX2(1.,aa); 00079 00084 iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]] = MIN2(.01,iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]]); 00085 iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]] = 0.; 00086 iso.pestrk[ipISO][nelem][iso.nLyaLevel[ipISO]][0] = 00087 iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]]*Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Aul; 00088 00089 00090 /* >>chng 06 aug 28, from numLevels_max to _local. */ 00091 for( ipHi=3; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00092 { 00093 if( Transitions[ipISO][nelem][ipHi][ipH1s].ipCont <= 0 ) 00094 continue; 00095 00096 iso.pestrk[ipISO][nelem][0][ipHi] = stark*iso.strkar[ipISO][nelem][0][ipHi]/2.*pow(MAX2(1., 00097 Transitions[ipISO][nelem][ipHi][ipH1s].Emis->TauIn),-0.75); 00098 00099 iso.pestrk[ipISO][nelem][0][ipHi] = MIN2(.01,iso.pestrk[ipISO][nelem][0][ipHi]); 00100 iso.pestrk[ipISO][nelem][ipHi][0] = Transitions[ipISO][nelem][ipHi][ipH1s].Emis->Aul* 00101 iso.pestrk[ipISO][nelem][0][ipHi]; 00102 } 00103 00104 /* zero out rates above iso.numLevels_local */ 00105 for( ipHi=iso.numLevels_local[ipISO][nelem]; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00106 { 00107 iso.pestrk[ipISO][nelem][0][ipHi] = 0.; 00108 iso.pestrk[ipISO][nelem][ipHi][0] = 0.; 00109 } 00110 00111 /* all other lines */ 00112 for( ipLo=ipH2s; ipLo < (iso.numLevels_local[ipISO][nelem] - 1); ipLo++ ) 00113 { 00114 for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00115 { 00116 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 ) 00117 continue; 00118 00119 aa = stark*iso.strkar[ipISO][nelem][ipLo][ipHi]* 00120 pow(MAX2(1.,Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn),-0.75); 00121 iso.pestrk[ipISO][nelem][ipLo][ipHi] = 00122 (realnum)MIN2(.01,aa); 00123 00124 iso.pestrk[ipISO][nelem][ipHi][ipLo] = Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul* 00125 iso.pestrk[ipISO][nelem][ipLo][ipHi]; 00126 } 00127 } 00128 00129 /* zero out rates above iso.numLevels_local */ 00130 for( ipLo=(iso.numLevels_local[ipISO][nelem] - 1); ipLo<(iso.numLevels_max[ipISO][nelem] - 1); ipLo++ ) 00131 { 00132 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00133 { 00134 iso.pestrk[ipISO][nelem][ipLo][ipHi] = 0.; 00135 iso.pestrk[ipISO][nelem][ipHi][ipLo] = 0.; 00136 } 00137 } 00138 } 00139 } 00140 00141 return; 00142 }