cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_recom_effic.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 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 "rt.h"
6 #include "rfield.h"
7 #include "phycon.h"
8 #include "opacity.h"
9 #include "rt_escprob.h"
10 #include "cosmology.h"
11 
12 double RT_recom_effic(long int ip)
13 {
14  long int i;
15  double dEner,
16  denom,
17  escin,
18  escout,
19  hnukt,
20  receff_v,
21  sum,
22  tin,
23  tout;
24 
25  DEBUG_ENTRY( "RT_recom_effic()" );
26 
27  /* escape probability function for continua,
28  * formally correct for photoelectric absorption only */
29 
30  ASSERT( ip > 0 && ip <= rfield.nflux_with_check );
31 
32  if( ip > rfield.nflux )
33  {
34  /* >>chng 01 dec 18, return had been zero, but this did not
35  * work for case where gas much hotter than continuum, as in a
36  * coronal plasma. change to return of unity */
37  receff_v = 1.;
38  return( receff_v );
39  }
40 
41  /* bug in following statement unocvered June 93 S. Schaefer */
42  hnukt = TE1RYD*rfield.anu(ip-1)/phycon.te;
43 
44  /* rfield.chDffTrns = "OU2" by default */
45  /* inward optical depth and escape prob */
46  if( strcmp(rfield.chDffTrns,"OSS") == 0 )
47  {
48  /* which side of Lyman limit is this? */
49  if( rfield.anu(ip) > 0.99 )
50  {
51  /* this is a simple OTS, turned on with DIFFUSE OTS SIMPLE */
52  receff_v = SMALLFLOAT;
53  }
54  else
55  {
56  receff_v = 1.;
57  }
58  }
59  else if( strcmp(rfield.chDffTrns,"OTS") == 0 )
60  {
61  tin = opac.TauAbsGeo[0][ip-1];
62  if( tin < 5. )
63  {
64  escin = esccon(tin,hnukt);
65  }
66  else
67  {
68  escin = 1e-4;
69  }
70 
71  /* outward optical depth */
72  tout = opac.TauAbsGeo[1][ip-1] - tin;
73 
74  if( iteration > 1 )
75  {
76  /* check whether we have overrun the optical depth scale */
77  if( tout > 0. )
78  {
79  /* good optical depths in both directions, take mean */
80  if( tout < 5. )
81  {
82  escout = esccon(tout,hnukt);
83  }
84  else
85  {
86  escout = 1e-4;
87  }
88  receff_v = 0.5*(escin + escout);
89  }
90  else
91  {
92  /* >>chng 91 apr add logic to prevent big change in
93  * esc prob, resulting in terminal oscillations, when optical
94  * depth scale overrun
95  * tau was negative, use 5% of inward optical depth */
96  escout = esccon(tin*0.05,hnukt);
97  receff_v = 0.5*(escin + escout);
98  }
99  }
100  else
101  {
102  receff_v = escin;
103  }
104  }
105  else if( strcmp(rfield.chDffTrns,"OU1") == 0 )
106  {
107  receff_v = opac.ExpZone[ip+1];
108  }
109  else if( strcmp(rfield.chDffTrns,"OU2") == 0 )
110  {
111  /* this is the default rt method, as set in zero
112  * E2TauAbsFace is optical depth to illuminated face */
113  receff_v = opac.E2TauAbsFace[ip+1];
114  }
115  else if( strcmp(rfield.chDffTrns,"OU3") == 0 )
116  {
117  receff_v = 1.;
118  }
119  else if( strcmp(rfield.chDffTrns,"OU4") == 0 )
120  {
121  /* this cannot happen, was the former outward treat
122  * optical depth for this zone */
123  if( rfield.ContBoltz[ip-1] > 0. )
124  {
125  i = ip;
126  dEner = phycon.te/TE1RYD*0.5;
127  sum = 0.;
128  denom = 0.;
129  while( rfield.ContBoltz[i-1] > 0. &&
130  rfield.anu(i-1)-rfield.anu(ip-1) < (realnum)dEner &&
131  i <= rfield.nflux )
132  {
133  sum += rfield.ContBoltz[i-1]*opac.tmn[i-1];
134  denom += rfield.ContBoltz[i-1];
135  i += 1;
136  }
137  receff_v = sum/denom;
138  }
139  else
140  {
141  receff_v = opac.tmn[ip-1];
142  }
143  }
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 */
168  EffectiveThickness = SPEEDLIGHT / dvdr;
169  tau = (realnum)(OpacityEffective * EffectiveThickness);
170  }
171  else
172  tau = opac.taumin;
173 
174  tau = MAX2((double)opac.taumin,tau);
175 
176  ASSERT( tau >= 0. );
177 
178  if( tau < 1e-5 )
179  receff_v = (1. - tau/2.);
180  else
181  receff_v = (1. - dsexp( tau ) )/ tau;
182  ASSERT( receff_v >= 0.f );
183  ASSERT( receff_v <= 1.f );
184  }
185  else
186  {
187  fprintf( ioQQQ, " RECEFF does not understand the transfer method=%3.3s\n",
188  rfield.chDffTrns );
190  }
191 
192  receff_v = MAX2((double)opac.otsmin,receff_v);
193  /* can get epsilon above unity on cray */
194  receff_v = MIN2(1.,receff_v);
195  return( receff_v );
196 }
realnum otsmin
Definition: opacity.h:314
double RT_recom_effic(long int ip)
vector< double, allocator_avx< double > > ContBoltz
Definition: rfield.h:124
double * opacity_abs
Definition: opacity.h:104
t_opac opac
Definition: opacity.cpp:5
const realnum SMALLFLOAT
Definition: cpu.h:246
char chDffTrns[4]
Definition: rfield.h:217
bool lgOpacityFine
Definition: rfield.h:404
t_phycon phycon
Definition: phycon.cpp:6
FILE * ioQQQ
Definition: cddefines.cpp:7
#define MIN2(a, b)
Definition: cddefines.h:803
bool lgDo
Definition: cosmology.h:44
double anu(size_t i) const
Definition: mesh.h:120
double dsexp(double x)
Definition: service.cpp:1038
long int nflux_with_check
Definition: rfield.h:49
long int iteration
Definition: cddefines.cpp:16
double esccon(double tau, double hnukt)
Definition: rt_escprob.cpp:552
realnum redshift_current
Definition: cosmology.h:26
t_rfield rfield
Definition: rfield.cpp:9
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:482
realnum * fine_opac_zone
Definition: rfield.h:391
#define ASSERT(exp)
Definition: cddefines.h:613
long nfine
Definition: rfield.h:387
vector< double, allocator_avx< double > > ExpZone
Definition: opacity.h:133
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
realnum * E2TauAbsFace
Definition: opacity.h:137
t_cosmology cosmology
Definition: cosmology.cpp:8
long int * ipnt_coarse_2_fine
Definition: rfield.h:382
realnum ** TauAbsGeo
Definition: opacity.h:91
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
realnum GetHubbleFactor(realnum z)
Definition: cosmology.cpp:10
realnum * tmn
Definition: opacity.h:149
double te
Definition: phycon.h:21
long int nflux
Definition: rfield.h:46
realnum taumin
Definition: opacity.h:167