cloudy trunk
Loading...
Searching...
No Matches
rt_recom_effic.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*RT_recom_effic generate escape probability function for continua, */
4#include "cddefines.h"
5#include "physconst.h"
6#include "rfield.h"
7#include "phycon.h"
8#include "opacity.h"
9#include "rt.h"
10
11double RT_recom_effic(long int ip)
12{
13 long int i;
14 double dEner,
15 denom,
16 escin,
17 escout,
18 hnukt,
19 receff_v,
20 sum,
21 tin,
22 tout;
23
24 DEBUG_ENTRY( "RT_recom_effic()" );
25
26 /* escape probability function for continua,
27 * formally correct for photoelectric absorption only */
28
29 ASSERT( ip > 0 && ip <= rfield.nupper );
30
31 if( ip > rfield.nflux )
32 {
33 /* >>chng 01 dec 18, return had been zero, but this did not
34 * work for case where gas much hotter than continuum, as in a
35 * coronal plasma. change to return of unity */
36 receff_v = 1.;
37 return( receff_v );
38 }
39
40 /* bug in following statement unocvered June 93 S. Schaefer */
41 hnukt = TE1RYD*rfield.anu[ip-1]/phycon.te;
42
43 /* rfield.chDffTrns = "OU2" by default */
44 /* inward optical depth and escape prob */
45 if( strcmp(rfield.chDffTrns,"OSS") == 0 )
46 {
47 /* which side of Lyman limit is this? */
48 if( rfield.anu[ip] > 0.99 )
49 {
50 /* this is a simple OTS, turned on with DIFFUSE OTS SIMPLE */
51 receff_v = SMALLFLOAT;
52 }
53 else
54 {
55 receff_v = 1.;
56 }
57 }
58 else if( strcmp(rfield.chDffTrns,"OTS") == 0 )
59 {
60 tin = opac.TauAbsGeo[0][ip-1];
61 if( tin < 5. )
62 {
63 escin = esccon(tin,hnukt);
64 }
65 else
66 {
67 escin = 1e-4;
68 }
69
70 /* outward optical depth */
71 tout = opac.TauAbsGeo[1][ip-1] - tin;
72
73 if( iteration > 1 )
74 {
75 /* check whether we have overrun the optical depth scale */
76 if( tout > 0. )
77 {
78 /* good optical depths in both directions, take mean */
79 if( tout < 5. )
80 {
81 escout = esccon(tout,hnukt);
82 }
83 else
84 {
85 escout = 1e-4;
86 }
87 receff_v = 0.5*(escin + escout);
88 }
89 else
90 {
91 /* >>chng 91 apr add logic to prevent big change in
92 * esc prob, resulting in terminal oscillations, when optical
93 * depth scale overrun
94 * tau was negative, use 5% of inward optical depth */
95 escout = esccon(tin*0.05,hnukt);
96 receff_v = 0.5*(escin + escout);
97 }
98 }
99 else
100 {
101 receff_v = escin;
102 }
103 }
104 else if( strcmp(rfield.chDffTrns,"OU1") == 0 )
105 {
106 receff_v = opac.ExpZone[ip+1];
107 }
108 else if( strcmp(rfield.chDffTrns,"OU2") == 0 )
109 {
110 /* this is the default rt method, as set in zero
111 * E2TauAbsFace is optical depth to illuminated face */
112 receff_v = opac.E2TauAbsFace[ip+1];
113 }
114 else if( strcmp(rfield.chDffTrns,"OU3") == 0 )
115 {
116 receff_v = 1.;
117 }
118 else if( strcmp(rfield.chDffTrns,"OU4") == 0 )
119 {
120 /* this cannot happen, was the former outward treat
121 * optical depth for this zone */
122 if( rfield.ContBoltz[ip-1] > 0. )
123 {
124 i = ip;
125 dEner = phycon.te/TE1RYD*0.5;
126 sum = 0.;
127 denom = 0.;
128 while( rfield.ContBoltz[i-1] > 0. &&
129 rfield.anu[i-1]-rfield.anu[ip-1] < (realnum)dEner &&
130 i <= rfield.nflux )
131 {
132 sum += rfield.ContBoltz[i-1]*opac.tmn[i-1];
133 denom += rfield.ContBoltz[i-1];
134 i += 1;
135 }
136 receff_v = sum/denom;
137 }
138 else
139 {
140 receff_v = opac.tmn[ip-1];
141 }
142 }
143#if 0
144 else if( strcmp(rfield.chDffTrns,"SOB") == 0 )
145 {
146 long int ipRecombEdgeFine = rfield.ipnt_coarse_2_fine[ip];
147 double OpacityEffective, EffectiveThickness;
148 realnum tau;
149
150 /* find line center opacity - use fine opacity if array indices are OK */
151 if( ipRecombEdgeFine>=0 && ipRecombEdgeFine<rfield.nfine && rfield.lgOpacityFine )
152 {
153 /* use fine opacities fine grid fine mesh to get optical depth
154 * to continuum source */
155 /* total line center optical depth, all overlapping lines included */
156 OpacityEffective = rfield.fine_opac_zone[ipRecombEdgeFine];
157 }
158 else
159 {
160 OpacityEffective = opac.opacity_abs[ip];
161 }
162
163 if( cosmology.lgDo )
164 {
165 /* dv/dr (s-1), equal to dv/dt / v */
166 /* in this case, dv/dr is just the Hubble factor */
167 realnum dvdr = GetHubbleFactor(cosmology.redshift_current);
168 /* here, we assume that recombining electrons have energy kT
169 * and that photons resulting from recombination must be redshifted
170 * by kT before they are too weak to photoionize the same lower level again */
171 realnum width_to_shift = phycon.te_ryd/(rfield.anu[ip]+phycon.te_ryd) * SPEEDLIGHT;
172 EffectiveThickness = width_to_shift / dvdr;
173 tau = (realnum)(OpacityEffective * EffectiveThickness);
174 }
175 else
176 tau = opac.taumin;
177
178 tau = MAX2((double)opac.taumin,tau);
179
180 ASSERT( tau >= 0. );
181
182 if( tau < 1e-5 )
183 receff_v = (1. - tau/2.);
184 else
185 receff_v = (1. - sexp( tau ) )/ tau;
186 ASSERT( receff_v >= 0.f );
187 ASSERT( receff_v <= 1.f );
188 }
189#endif
190 else
191 {
192 fprintf( ioQQQ, " RECEFF does not understand the transfer method=%3.3s\n",
193 rfield.chDffTrns );
195 }
196
197 receff_v = MAX2((double)opac.otsmin,receff_v);
198 /* can get epsilon above unity on cray */
199 receff_v = MIN2(1.,receff_v);
200 return( receff_v );
201}
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define MIN2
Definition cddefines.h:761
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_cosmology cosmology
Definition cosmology.cpp:11
realnum GetHubbleFactor(realnum z)
Definition cosmology.cpp:13
const realnum SMALLFLOAT
Definition cpu.h:191
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double TE1RYD
Definition physconst.h:183
t_rfield rfield
Definition rfield.cpp:8
double esccon(double tau, double hnukt)
double RT_recom_effic(long int ip)