cloudy trunk
Loading...
Searching...
No Matches
cont_pump.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/*DrvContPump local continuum pumping rate radiative transfer for all lines */
4/*con_pump_op routine used to get continuum pumping of lines
5 * used in DrvContPump in call to qg32 */
6#include "cddefines.h"
7#include "rfield.h"
8#include "doppvel.h"
9#include "radius.h"
10#include "continuum.h"
11#include "transition.h"
12#include "rt.h"
13
14/*con_pump_op routine used to get continuum pumping of lines
15 * used in DrvContPump in call to qg32 */
17{
18public:
19 /* damping constant used for pumping */
21 /* variable used for inward optical depth for pumping */
23
24 double operator() (double x)
25 {
26 realnum v, rx = realnum(x);
27 VoigtH(damp,&rx,&v,1);
28 double opfun_v = sexp(PumpTau*v)*v;
29 return opfun_v;
30 }
31};
32
33/* fit to results for tau less than 10 */
34inline double fitted(double t)
35{
36 return (0.98925439 + 0.084594094*t)/(1. + t*(0.64794212 + t*0.44743976));
37}
38
40double DrvContPump(const TransitionProxy& t, realnum DopplerWidth)
41{
42 double a0,
43 ContPump_v,
44 tau,
45 yinc1,
46 yinc2;
47
48 DEBUG_ENTRY( "DrvContPump()" );
49
50 if( !rfield.lgInducProcess )
51 {
52 /* option to turn off continuum pumping with no fluorescence */
53 ContPump_v = 0.;
54 }
55 else
56 {
57 /* tau used will be optical depth in center of next zone */
58 tau = t.Emis().TauIn() + t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth * radius.dRNeff;
59 /* compute pumping probability */
60 if( tau <= 10. )
61 {
62 /* for tau<10 a does not matter, and one fit for all */
63 ContPump_v = fitted(tau);
64 }
65 else if( tau > 1e6 )
66 {
67 /* this far in winds line opacity well below electron scattering
68 * so ignore for this problem */
69 ContPump_v = 0.;
70 }
71 else
72 {
74 if( t.Emis().iRedisFun() > 0 )
75 {
76 func.damp = t.Emis().damp();
77 }
78 else
79 {
80 func.damp = 0.;
81 }
82 func.PumpTau = (realnum)tau;
83
85 static const double BREAK = 3.;
86 yinc1 = con_pump_op.sum(0.,BREAK,func);
87 yinc2 = con_pump_op.sum(BREAK,100.,func);
88
89 a0 = 0.886227*(1. + func.damp);
90 ContPump_v = (yinc1 + yinc2)/a0;
91 }
92 }
93
94 /* EscProb is escape probability, will not allow ContPump to be greater than it
95 * on second iteration with thick lines, pump prob=1 and esc=0.5
96 * ContPump = MIN( ContPump , t.t(ipLnEscP) )
97 * */
98 return ContPump_v;
99}
static double a0[83]
sys_float sexp(sys_float x)
Definition service.cpp:914
float realnum
Definition cddefines.h:103
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
realnum & damp() const
Definition emission.h:563
double & PopOpc() const
Definition emission.h:603
int & iRedisFun() const
Definition emission.h:403
realnum & TauIn() const
Definition emission.h:423
realnum & opacity() const
Definition emission.h:593
double sum(double min, double max, Integrand func)
Definition cddefines.h:1531
EmissionList::reference Emis() const
Definition transition.h:408
double operator()(double x)
Definition cont_pump.cpp:24
double fitted(double t)
Definition cont_pump.cpp:34
double DrvContPump(const TransitionProxy &t, realnum DopplerWidth)
Definition cont_pump.cpp:40
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
static const double BREAK
void VoigtH(realnum a, const realnum v[], realnum y[], int n)
Definition thirdparty.h:350