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_recom_effic generate escape probability function for continua, */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "rfield.h" 00007 #include "phycon.h" 00008 #include "opacity.h" 00009 #include "rt.h" 00010 00011 double RT_recom_effic(long int ip) 00012 { 00013 long int i; 00014 double dEner, 00015 denom, 00016 escin, 00017 escout, 00018 hnukt, 00019 receff_v, 00020 sum, 00021 tin, 00022 tout; 00023 00024 DEBUG_ENTRY( "RT_recom_effic()" ); 00025 00026 /* escape probability function for continua, 00027 * formally correct for photoelectric absorption only */ 00028 00029 ASSERT( ip > 0 && ip <= rfield.nupper ); 00030 00031 if( ip > rfield.nflux ) 00032 { 00033 /* >>chng 01 dec 18, return had been zero, but this did not 00034 * work for case where gas much hotter than continuum, as in a 00035 * coronal plasma. change to return of unity */ 00036 receff_v = 1.; 00037 return( receff_v ); 00038 } 00039 00040 /* bug in following statement unocvered June 93 S. Schaefer */ 00041 hnukt = TE1RYD*rfield.anu[ip-1]/phycon.te; 00042 00043 /* rfield.chDffTrns = "OU2" by default */ 00044 /* inward optical depth and escape prob */ 00045 if( strcmp(rfield.chDffTrns,"OSS") == 0 ) 00046 { 00047 /* which side of Lyman limit is this? */ 00048 if( rfield.anu[ip] > 0.99 ) 00049 { 00050 /* this is a simple OTS, turned on with DIFFUSE OTS SIMPLE */ 00051 receff_v = SMALLFLOAT; 00052 } 00053 else 00054 { 00055 receff_v = 1.; 00056 } 00057 } 00058 else if( strcmp(rfield.chDffTrns,"OTS") == 0 ) 00059 { 00060 tin = opac.TauAbsGeo[0][ip-1]; 00061 if( tin < 5. ) 00062 { 00063 escin = esccon(tin,hnukt); 00064 } 00065 else 00066 { 00067 escin = 1e-4; 00068 } 00069 00070 /* outward optical depth */ 00071 tout = opac.TauAbsGeo[1][ip-1] - tin; 00072 00073 if( iteration > 1 ) 00074 { 00075 /* check whether we have overrun the optical depth scale */ 00076 if( tout > 0. ) 00077 { 00078 /* good optical depths in both directions, take mean */ 00079 if( tout < 5. ) 00080 { 00081 escout = esccon(tout,hnukt); 00082 } 00083 else 00084 { 00085 escout = 1e-4; 00086 } 00087 receff_v = 0.5*(escin + escout); 00088 } 00089 else 00090 { 00091 /* >>chng 91 apr add logic to prevent big change in 00092 * esc prob, resulting in terminal oscillations, when optical 00093 * depth scale overrun 00094 * tau was negative, use 5% of inward optical depth */ 00095 escout = esccon(tin*0.05,hnukt); 00096 receff_v = 0.5*(escin + escout); 00097 } 00098 } 00099 else 00100 { 00101 receff_v = escin; 00102 } 00103 } 00104 else if( strcmp(rfield.chDffTrns,"OU1") == 0 ) 00105 { 00106 receff_v = opac.ExpZone[ip+1]; 00107 } 00108 else if( strcmp(rfield.chDffTrns,"OU2") == 0 ) 00109 { 00110 /* this is the default rt method, as set in zero 00111 * E2TauAbsFace is optical depth to illuminated face */ 00112 receff_v = opac.E2TauAbsFace[ip+1]; 00113 } 00114 else if( strcmp(rfield.chDffTrns,"OU3") == 0 ) 00115 { 00116 receff_v = 1.; 00117 } 00118 else if( strcmp(rfield.chDffTrns,"OU4") == 0 ) 00119 { 00120 /* this cannot happen, was the former outward treat 00121 * optical depth for this zone */ 00122 if( rfield.ContBoltz[ip-1] > 0. ) 00123 { 00124 i = ip; 00125 dEner = phycon.te/TE1RYD*0.5; 00126 sum = 0.; 00127 denom = 0.; 00128 while( rfield.ContBoltz[i-1] > 0. && 00129 rfield.anu[i-1]-rfield.anu[ip-1] < (realnum)dEner && 00130 i <= rfield.nflux ) 00131 { 00132 sum += rfield.ContBoltz[i-1]*opac.tmn[i-1]; 00133 denom += rfield.ContBoltz[i-1]; 00134 i += 1; 00135 } 00136 receff_v = sum/denom; 00137 } 00138 else 00139 { 00140 receff_v = opac.tmn[ip-1]; 00141 } 00142 } 00143 else 00144 { 00145 fprintf( ioQQQ, " RECEFF does not understand the transfer method=%3.3s\n", 00146 rfield.chDffTrns ); 00147 cdEXIT(EXIT_FAILURE); 00148 } 00149 00150 receff_v = MAX2((double)opac.otsmin,receff_v); 00151 /* can get epsilon above unity on cray */ 00152 receff_v = MIN2(1.,receff_v); 00153 return( receff_v ); 00154 }