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_continuum_shield_fcn computing continuum shielding due to single line */ 00004 /*conpmp local continuum pumping rate radiative transfer for all lines */ 00005 /*conpmp local continuum pumping rate radiative transfer for all lines */ 00006 #include "cddefines.h" 00007 #include "rt.h" 00008 /*conpmp local continuum pumping rate radiative transfer for all lines */ 00009 STATIC double conpmp(transition * t); 00010 STATIC inline double FITTED( double t ); 00011 static double PumpDamp , PumpTau; 00012 00013 /*rt_continuum_shield_fcn computing continuum shielding due to single line */ 00014 double RT_continuum_shield_fcn( transition *t ) 00015 { 00016 double value; 00017 00018 DEBUG_ENTRY( "rt_continuum_shield_fcn()" ); 00019 00020 value = -1.; 00021 00022 if( rt.nLineContShield == LINE_CONT_SHIELD_PESC ) 00023 { 00024 /* set continuum shielding pesc - shieding based on escape probability */ 00025 if( t->Emis->iRedisFun == ipPRD ) 00026 { 00027 value = esc_PRD_1side(t->Emis->TauCon,t->Emis->damp); 00028 } 00029 else if( t->Emis->iRedisFun == ipCRD ) 00030 { 00031 value = esca0k2(t->Emis->TauCon); 00032 } 00033 else if( t->Emis->iRedisFun == ipCRDW ) 00034 { 00035 value = esc_CRDwing_1side(t->Emis->TauCon,t->Emis->damp); 00036 } 00037 else if( t->Emis->iRedisFun == ipLY_A ) 00038 { 00039 value = esc_PRD_1side(t->Emis->TauCon,t->Emis->damp); 00040 } 00041 else 00042 TotalInsanity(); 00043 } 00044 else if( rt.nLineContShield == LINE_CONT_SHIELD_FEDERMAN ) 00045 { 00046 /* set continuum shielding Federman - this is the default */ 00047 double core, wings; 00048 00049 /* these expressions implement the appendix of 00050 * >>refer line shielding Federman, S.R., Glassgold, A.E., & 00051 * >>refercon Kwan, J. 1979, ApJ, 227, 466 */ 00052 /* doppler core - equation A8 */ 00053 if( t->Emis->TauCon < 2. ) 00054 { 00055 core = sexp( t->Emis->TauCon * 0.66666 ); 00056 } 00057 else if( t->Emis->TauCon < 10. ) 00058 { 00059 core = 0.638 * pow(t->Emis->TauCon,(realnum)-1.25f ); 00060 } 00061 else if( t->Emis->TauCon < 100. ) 00062 { 00063 core = 0.505 * pow(t->Emis->TauCon,(realnum)-1.15f ); 00064 } 00065 else 00066 { 00067 core = 0.344 * pow(t->Emis->TauCon,(realnum)-1.0667f ); 00068 } 00069 00070 /* do we add damping wings? */ 00071 wings = 0.; 00072 if( t->Emis->TauCon > 1.f && t->Emis->damp>0. ) 00073 { 00074 /* equation A6 */ 00075 double t1 = 3.02*pow(t->Emis->damp*1e3,-0.064 ); 00076 double u1 = sqrt(t->Emis->TauCon*t->Emis->damp )/SDIV(t1); 00077 wings = t->Emis->damp/SDIV(t1)/sqrt( 0.78540 + POW2(u1) ); 00078 /* add very large optical depth tail to converge this with respect 00079 * to escape probabilities - if this function falls off more slowly 00080 * than escape probability then upper level will become overpopulated. 00081 * original paper was not intended for this regime */ 00082 if( t->Emis->TauCon>1e7 ) 00083 wings *= pow( t->Emis->TauCon/1e7,-1.1 ); 00084 } 00085 value = core + wings; 00086 /* some x-ray lines have vastly large damping constants, greater than 1. 00087 * in these cases the previous wings value does not work - approximation 00088 * is for small damping constant - do not let pump efficiency exceed unity 00089 * in this case */ 00090 if( t->Emis->TauCon>0. ) 00091 value = MIN2(1., value ); 00092 } 00093 else if( rt.nLineContShield == LINE_CONT_SHIELD_FERLAND ) 00094 { 00095 /* set continuum shielding ferland */ 00096 value = conpmp( t ); 00097 } 00098 else if( rt.nLineContShield == 0 ) 00099 { 00100 /* set continuum shielding none */ 00101 value = 1.; 00102 } 00103 else 00104 { 00105 TotalInsanity(); 00106 } 00107 00108 /* the returned pump shield function must be greater than zero, 00109 * and less than 1 if a maser did not occur */ 00110 ASSERT( value>=0 && (value<=1.||t->Emis->TauCon<0.) ); 00111 00112 return value; 00113 } 00114 00115 /*opfun routine used to get continuum pumping of lines 00116 * used in conpmp in call to qg32 */ 00117 STATIC double opfun(double x) 00118 { 00119 double opfun_v, 00120 v; 00121 00122 DEBUG_ENTRY( "opfun()" ); 00123 00124 v = vfun(PumpDamp,x); 00125 opfun_v = sexp(PumpTau*v)*v; 00126 return( opfun_v ); 00127 } 00128 00129 static const double BREAK = 3.; 00130 /* fit to results for tau less than 10 */ 00131 // #define FITTED(t) ((0.98925439 + 0.084594094*(t))/(1. + (t)*(0.64794212 + (t)*0.44743976))) 00132 STATIC inline double FITTED( double t ) 00133 { 00134 return (0.98925439 + 0.084594094*t)/(1. + t*(0.64794212 + t*0.44743976)); 00135 } 00136 00137 /*conpmp local continuum pumping rate radiative transfer for all lines */ 00138 STATIC double conpmp(transition * t) 00139 { 00140 double a0, 00141 conpmp_v, 00142 tau, 00143 yinc1, 00144 yinc2; 00145 00146 DEBUG_ENTRY( "conpmp()" ); 00147 00148 /* tau used will be optical depth in center of next zone 00149 * >>chng 96 july 6, had been ipLnTauIn, did not work when sphere static set */ 00150 tau = t->Emis->TauCon; 00151 /* compute pumping probability */ 00152 if( tau <= 10. ) 00153 { 00154 /* for tau<10 a does not matter, and one fit for all */ 00155 conpmp_v = FITTED(tau); 00156 } 00157 else if( tau > 1e6 ) 00158 { 00159 /* this far in winds line opacity well below electron scattering 00160 * so ignore for this problem */ 00161 conpmp_v = 0.; 00162 } 00163 else 00164 { 00165 /* following two are passed on to later subs */ 00166 PumpDamp = t->Emis->damp; 00167 PumpTau = tau; 00168 a0 = 0.886227*(1. + PumpDamp); 00169 yinc1 = qg32(0.,BREAK,opfun); 00170 yinc2 = qg32(BREAK,100.,opfun); 00171 conpmp_v = (yinc1 + yinc2)/a0; 00172 } 00173 00174 /* EscProb is escape probability, will not allow conpmp to be greater than it 00175 * on second iteration with thick lines, pump prob=1 and esc=0.5 00176 * conpmp = MIN( conpmp , t->t(ipLnEscP) ) 00177 * */ 00178 return( conpmp_v ); 00179 }