cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
two_photon.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 
4 #include "cddefines.h"
5 #include "ipoint.h"
6 #include "rfield.h"
7 #include "transition.h"
8 #include "two_photon.h"
9 
10 void TwoPhotonSetup( vector<two_photon> &tnu_vec, const long &ipHi, const long &ipLo, const double &Aul, const TransitionProxy &tr, const long ipISO, const long nelem )
11 {
12  DEBUG_ENTRY( "TwoPhotonSetup()" );
13 
14  tnu_vec.resize( tnu_vec.size() + 1 );
15  two_photon &tnu = tnu_vec.back();
16 
17  tnu.ipHi = ipHi;
18  tnu.ipLo = ipLo;
19  tnu.AulTotal = Aul;
20  tnu.Pop = &(*tr.Hi()).Pop();
21  tnu.E2nu = tr.EnergyRyd();
22  tnu.induc_dn_max = 0.;
23 
24  /* pointer to the Lya energy */
25  tnu.ipTwoPhoE = ipoint(tnu.E2nu);
26  while( rfield.anu(tnu.ipTwoPhoE) > tnu.E2nu )
27  {
28  --tnu.ipTwoPhoE;
29  }
30  tnu.ipSym2nu.resize( tnu.ipTwoPhoE );
31  tnu.As2nu.resize( tnu.ipTwoPhoE );
32  tnu.local_emis.resize( tnu.ipTwoPhoE );
33 
34  /* >>chng 02 aug 14, change upper limit to full Lya energy */
35  for( long i=0; i < tnu.ipTwoPhoE; i++ )
36  {
37  /* energy is symmetric energy, the other side of half E2nu */
38  double energy = tnu.E2nu - rfield.anu(i);
39  /* this is needed since mirror image of cell next to two-nu energy
40  * may be slightly negative */
41  energy = MAX2( energy, rfield.anumax(0) );
42  /* find index for this symmetric energy */
43  tnu.ipSym2nu[i] = ipoint(energy);
44  while( rfield.anu(tnu.ipSym2nu[i]) > tnu.E2nu ||
45  tnu.ipSym2nu[i] >= tnu.ipTwoPhoE)
46  {
47  --tnu.ipSym2nu[i];
48  }
49  ASSERT( tnu.ipSym2nu[i] >= 0 );
50  }
51 
52  double SumShapeFunction = 0., Renorm= 0.;
53 
54  /* ipTwoPhoE is the cell holding the 2nu energy itself, and we do not want
55  * to include that in the following sum */
56  ASSERT( rfield.anu(tnu.ipTwoPhoE-1)<=tnu.E2nu );
57  for( long i=0; i < tnu.ipTwoPhoE; i++ )
58  {
59  double ShapeFunction;
60 
61  ShapeFunction = atmdat_2phot_shapefunction( rfield.anu(i)/tnu.E2nu, ipISO, nelem )*rfield.widflx(i)/tnu.E2nu;
62  SumShapeFunction += ShapeFunction;
63 
64  /* >>refer HI 2nu Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
65  /* As2nu will add up to the A, so its units are s-1 */
66  tnu.As2nu[i] = (realnum)( tnu.AulTotal * ShapeFunction );
67  }
68 
69  /* The spline function in twophoton.c causes a bit of an error in the integral of the
70  * shape function. So we renormalize the integral to 1. */
71  Renorm = 1./SumShapeFunction;
72 
73  for( long i=0; i < tnu.ipTwoPhoE; i++ )
74  {
75  tnu.As2nu[i] *= (realnum)Renorm;
76  }
77 
78  /* The result should be VERY close to 1. */
79  ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );
80 
81  return;
82 }
83 
84 void CalcTwoPhotonRates( two_photon& tnu, bool lgDoInduced )
85 {
86  DEBUG_ENTRY( "CalcTwoPhotonRates()" );
87 
88  /* this could fail when pops very low and Pop2Ion is zero */
89  ASSERT( tnu.ipTwoPhoE>0 && tnu.E2nu>0. );
90 
91  double sum = 0.;
92  tnu.induc_up = 0.;
93  tnu.induc_dn = 0.;
94  /* two photon emission, ipTwoPhoE is
95  * continuum array index for Lya energy */
96  ASSERT( rfield.anu(tnu.ipTwoPhoE-1) < 1.01*tnu.E2nu || rfield.anu(tnu.ipTwoPhoE-2)<tnu.E2nu );
97  for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
98  {
99  // We do not assert rfield.anu(nu)<=tnu.E2nu because the maximum
100  // index could be set to point to the next highest bin.
101 
102  // As2nu[nu] is transition probability A per bin
103  // So sum is the total transition probability
104  sum += tnu.As2nu[nu];
105 
106  // only include this if induced processes turned on,
107  // otherwise inconsistent with rate solver treatment.
108  if( lgDoInduced )
109  {
110  double rate_up = tnu.As2nu[nu] *
111  rfield.SummedOcc[nu] * rfield.SummedOcc[tnu.ipSym2nu[nu]-1];
112  tnu.induc_up += rate_up;
113  tnu.induc_dn += rate_up + tnu.As2nu[nu] *
114  (rfield.SummedOcc[nu] + rfield.SummedOcc[tnu.ipSym2nu[nu]-1]);
115  }
116  }
117 
118  /* a sanity check on the code, see Spitzer and Greenstein (1951), eqn 4. */
119  /* >>refer HI 2nu Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
120  ASSERT( fabs( 1.f - (realnum)sum/tnu.AulTotal ) < 0.01f );
121 
122  return;
123 }
124 
125 void CalcTwoPhotonEmission( two_photon& tnu, bool lgDoInduced )
126 {
127  DEBUG_ENTRY( "CalcTwoPhotonEmission()" );
128 
129  /* this could fail when pops very low and Pop2Ion is zero */
130  ASSERT( tnu.ipTwoPhoE>0 );
131 
132  /* two photon emission, ipTwoPhoE is
133  * continuum array index for Lya energy */
134  for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
135  {
136  // Pop has dimension cm^-3. The factor of 2 is for two photons per
137  // transition. Thus two_photon_emiss has dimension photons cm-3 s-1.
138  tnu.local_emis[nu] = 2.f * (realnum)(*tnu.Pop) * tnu.As2nu[nu];
139  }
140 
141  // only include this if induced processes turned on,
142  // otherwise inconsistent with rate solver treatment.
143  if( lgDoInduced )
144  {
145  for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
146  {
147  // this is the total rate (in this energy bin)
148  tnu.local_emis[nu] *= (1.f + rfield.SummedOcc[nu]) *
149  (1.f+rfield.SummedOcc[tnu.ipSym2nu[nu]-1]);
150  }
151  }
152 
153  return;
154 }
155 
156 /* option to print hydrogen and helium two-photon emission coefficients? */
157 void PrtTwoPhotonEmissCoef( const two_photon& tnu, const double& densityProduct )
158 {
159  DEBUG_ENTRY( "PrtTwoPhotonEmissCoef()" );
160 
161  fprintf( ioQQQ, "\ny\tGammaNot(2q)\n");
162 
163  for( long yTimes20=1; yTimes20<=10; yTimes20++ )
164  {
165  double y = yTimes20/20.;
166 
167  fprintf( ioQQQ, "%.3e\t", (double)y );
168 
169  long i = ipoint(y*tnu.E2nu);
170  fprintf( ioQQQ, "%.3e\n",
171  8./3.*HPLANCK*(*tnu.Pop)/densityProduct*y*tnu.As2nu[i]*tnu.E2nu/rfield.widflx(i) );
172  }
173 
174  return;
175 }
176 
double induc_dn
Definition: two_photon.h:53
double induc_up
Definition: two_photon.h:51
realnum AulTotal
Definition: two_photon.h:38
double widflx(size_t i) const
Definition: mesh.h:156
double induc_dn_max
Definition: two_photon.h:55
long ipTwoPhoE
Definition: two_photon.h:41
realnum * SummedOcc
Definition: rfield.h:163
FILE * ioQQQ
Definition: cddefines.cpp:7
double anu(size_t i) const
Definition: mesh.h:120
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:15
double energy(const genericState &gs)
void PrtTwoPhotonEmissCoef(const two_photon &tnu, const double &densityProduct)
Definition: two_photon.cpp:157
t_rfield rfield
Definition: rfield.cpp:9
long ipHi
Definition: two_photon.h:35
float realnum
Definition: cddefines.h:124
double * Pop
Definition: two_photon.h:36
qList::iterator Hi() const
Definition: transition.h:435
vector< long > ipSym2nu
Definition: two_photon.h:44
void CalcTwoPhotonRates(two_photon &tnu, bool lgDoInduced)
Definition: two_photon.cpp:84
long ipLo
Definition: two_photon.h:35
#define ASSERT(exp)
Definition: cddefines.h:613
void TwoPhotonSetup(vector< two_photon > &tnu_vec, const long &ipHi, const long &ipLo, const double &Aul, const TransitionProxy &tr, const long ipISO, const long nelem)
Definition: two_photon.cpp:10
vector< realnum > local_emis
Definition: two_photon.h:48
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
double EnergyRyd() const
Definition: transition.h:95
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
vector< realnum > As2nu
Definition: two_photon.h:46
double anumax(size_t i) const
Definition: mesh.h:152
double E2nu
Definition: two_photon.h:37
void CalcTwoPhotonEmission(two_photon &tnu, bool lgDoInduced)
Definition: two_photon.cpp:125
double atmdat_2phot_shapefunction(double EbyE2nu, long ipISO, long nelem)