00001
00002
00003
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
00027
00028
00029 ASSERT( ip > 0 && ip <= rfield.nupper );
00030
00031 if( ip > rfield.nflux )
00032 {
00033
00034
00035
00036 receff_v = 1.;
00037 return( receff_v );
00038 }
00039
00040
00041 hnukt = TE1RYD*rfield.anu[ip-1]/phycon.te;
00042
00043
00044
00045 if( strcmp(rfield.chDffTrns,"OSS") == 0 )
00046 {
00047
00048 if( rfield.anu[ip] > 0.99 )
00049 {
00050
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
00071 tout = opac.TauAbsGeo[1][ip-1] - tin;
00072
00073 if( iteration > 1 )
00074 {
00075
00076 if( tout > 0. )
00077 {
00078
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
00092
00093
00094
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
00111
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
00121
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 #if 0
00144 else if( strcmp(rfield.chDffTrns,"SOB") == 0 )
00145 {
00146 long int ipRecombEdgeFine = rfield.ipnt_coarse_2_fine[ip];
00147 double OpacityEffective, EffectiveThickness;
00148 realnum tau;
00149
00150
00151 if( ipRecombEdgeFine>=0 && ipRecombEdgeFine<rfield.nfine && rfield.lgOpacityFine )
00152 {
00153
00154
00155
00156 OpacityEffective = rfield.fine_opac_zone[ipRecombEdgeFine];
00157 }
00158 else
00159 {
00160 OpacityEffective = opac.opacity_abs[ip];
00161 }
00162
00163 if( cosmology.lgDo )
00164 {
00165
00166
00167 realnum dvdr = GetHubbleFactor(cosmology.redshift_current);
00168
00169
00170
00171 realnum width_to_shift = phycon.te_ryd/(rfield.anu[ip]+phycon.te_ryd) * SPEEDLIGHT;
00172 EffectiveThickness = width_to_shift / dvdr;
00173 tau = (realnum)(OpacityEffective * EffectiveThickness);
00174 }
00175 else
00176 tau = opac.taumin;
00177
00178 tau = MAX2((double)opac.taumin,tau);
00179
00180 ASSERT( tau >= 0. );
00181
00182 if( tau < 1e-5 )
00183 receff_v = (1. - tau/2.);
00184 else
00185 receff_v = (1. - sexp( tau ) )/ tau;
00186 ASSERT( receff_v >= 0.f );
00187 ASSERT( receff_v <= 1.f );
00188 }
00189 #endif
00190 else
00191 {
00192 fprintf( ioQQQ, " RECEFF does not understand the transfer method=%3.3s\n",
00193 rfield.chDffTrns );
00194 cdEXIT(EXIT_FAILURE);
00195 }
00196
00197 receff_v = MAX2((double)opac.otsmin,receff_v);
00198
00199 receff_v = MIN2(1.,receff_v);
00200 return( receff_v );
00201 }