cloudy  trunk
 All Data Structures 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-2008 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 
11 double 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  else
144  {
145  fprintf( ioQQQ, " RECEFF does not understand the transfer method=%3.3s\n",
146  rfield.chDffTrns );
147  cdEXIT(EXIT_FAILURE);
148  }
149 
150  receff_v = MAX2((double)opac.otsmin,receff_v);
151  /* can get epsilon above unity on cray */
152  receff_v = MIN2(1.,receff_v);
153  return( receff_v );
154 }

Generated for cloudy by doxygen 1.8.4