cloudy trunk
Loading...
Searching...
No Matches
two_photon.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
4#include "cddefines.h"
5#include "ipoint.h"
6#include "rfield.h"
7#include "two_photon.h"
8
9void TwoPhotonSetup( vector<two_photon> &tnu_vec, const long &ipHi, const long &ipLo, const double &Aul, const TransitionProxy &tr, const long ipISO, const long nelem )
10{
11 DEBUG_ENTRY( "TwoPhotonSetup()" );
12
13 tnu_vec.resize( tnu_vec.size() + 1 );
14 two_photon &tnu = tnu_vec.back();
15
16 tnu.ipHi = ipHi;
17 tnu.ipLo = ipLo;
18 tnu.AulTotal = Aul;
19 tnu.Pop = &(*tr.Hi()).Pop();
20 tnu.E2nu = tr.EnergyRyd();
21 tnu.induc_dn_max = 0.;
22
23 /* pointer to the Lya energy */
24 tnu.ipTwoPhoE = ipoint(tnu.E2nu);
25 while( rfield.anu[tnu.ipTwoPhoE] > tnu.E2nu )
26 {
27 --tnu.ipTwoPhoE;
28 }
29 tnu.ipSym2nu.resize( tnu.ipTwoPhoE );
30 tnu.As2nu.resize( tnu.ipTwoPhoE );
31 tnu.local_emis.resize( tnu.ipTwoPhoE );
32
33 /* >>chng 02 aug 14, change upper limit to full Lya energy */
34 for( long i=0; i < tnu.ipTwoPhoE; i++ )
35 {
36 /* energy is symmetric energy, the other side of half E2nu */
37 double energy = tnu.E2nu - rfield.anu[i];
38 /* this is needed since mirror image of cell next to two-nu energy
39 * may be slightly negative */
40 energy = MAX2( energy, rfield.anu[0] + rfield.widflx[0]/2. );
41 /* find index for this symmetric energy */
42 tnu.ipSym2nu[i] = ipoint(energy);
43 while( rfield.anu[tnu.ipSym2nu[i]] > tnu.E2nu ||
44 tnu.ipSym2nu[i] >= tnu.ipTwoPhoE)
45 {
46 --tnu.ipSym2nu[i];
47 }
48 ASSERT( tnu.ipSym2nu[i] >= 0 );
49 }
50
51 double SumShapeFunction = 0., Renorm= 0.;
52
53 /* ipTwoPhoE is the cell holding the 2nu energy itself, and we do not want
54 * to include that in the following sum */
55 for( long i=0; i < tnu.ipTwoPhoE; i++ )
56 {
57 double ShapeFunction;
58
59 ASSERT( rfield.anu[i]<=tnu.E2nu );
60 ShapeFunction = atmdat_2phot_shapefunction( rfield.AnuOrg[i]/tnu.E2nu, ipISO, nelem )*rfield.widflx[i]/tnu.E2nu;
61 SumShapeFunction += ShapeFunction;
62
63 /* >>refer HI 2nu Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
64 /* As2nu will add up to the A, so its units are s-1 */
65 tnu.As2nu[i] = (realnum)( tnu.AulTotal * ShapeFunction );
66 }
67
68 /* The spline function in twophoton.c causes a bit of an error in the integral of the
69 * shape function. So we renormalize the integral to 1. */
70 Renorm = 1./SumShapeFunction;
71
72 for( long i=0; i < tnu.ipTwoPhoE; i++ )
73 {
74 tnu.As2nu[i] *= (realnum)Renorm;
75 }
76
77 /* The result should be VERY close to 1. */
78 ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );
79
80 return;
81}
82
83void CalcTwoPhotonRates( two_photon& tnu, bool lgDoInduced )
84{
85 DEBUG_ENTRY( "CalcTwoPhotonRates()" );
86
87 /* this could fail when pops very low and Pop2Ion is zero */
88 ASSERT( tnu.ipTwoPhoE>0 && tnu.E2nu>0. );
89
90 double sum = 0.;
91 tnu.induc_up = 0.;
92 tnu.induc_dn = 0.;
93 /* two photon emission, ipTwoPhoE is
94 * continuum array index for Lya energy */
95 for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
96 {
97 // We do not assert rfield.anu[nu]<=tnu.E2nu because the maximum
98 // index could be set to point to the next highest bin.
99 ASSERT( rfield.anu[nu] < 1.01*tnu.E2nu || rfield.anu[nu-1]<tnu.E2nu );
100
101 // As2nu[nu] is transition probability A per bin
102 // So sum is the total transition probability
103 sum += tnu.As2nu[nu];
104
105 // only include this if induced processes turned on,
106 // otherwise inconsistent with rate solver treatment.
107 if( lgDoInduced )
108 {
109 // factor 0.5 accounts for fact that photons must be anti-aligned.
110 double rate_up = 0.5 * 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
125void 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? */
157void 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 atmdat_2phot_shapefunction(double EbyE2nu, long ipISO, long nelem)
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
double EnergyRyd() const
Definition transition.h:83
qList::iterator Hi() const
Definition transition.h:396
vector< realnum > As2nu
Definition two_photon.h:46
vector< long > ipSym2nu
Definition two_photon.h:44
double induc_dn
Definition two_photon.h:53
double induc_up
Definition two_photon.h:51
long ipTwoPhoE
Definition two_photon.h:41
realnum AulTotal
Definition two_photon.h:38
vector< realnum > local_emis
Definition two_photon.h:48
double * Pop
Definition two_photon.h:36
double E2nu
Definition two_photon.h:37
double induc_dn_max
Definition two_photon.h:55
long ipoint(double energy_ryd)
UNUSED const double HPLANCK
Definition physconst.h:103
t_rfield rfield
Definition rfield.cpp:8
void CalcTwoPhotonEmission(two_photon &tnu, bool lgDoInduced)
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:9
void CalcTwoPhotonRates(two_photon &tnu, bool lgDoInduced)
void PrtTwoPhotonEmissCoef(const two_photon &tnu, const double &densityProduct)