00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "rt.h"
00008 #include "transition.h"
00009 #include "thirdparty.h"
00010
00011
00012 STATIC double conpmp(const TransitionProxy& t);
00013 STATIC inline double FITTED( double t );
00014
00015 namespace {
00016 class my_Integrand
00017 {
00018 public:
00019 double PumpDamp, PumpTau;
00020
00021 double operator() (double x)
00022 {
00023 realnum v, rx = realnum(x);
00024 VoigtH(realnum(PumpDamp),&rx,&v,1);
00025 double opfun_v = sexp(PumpTau*v)*v;
00026 return opfun_v;
00027 }
00028 };
00029 }
00030
00031
00032
00033 double RT_continuum_shield_fcn( const TransitionProxy& t )
00034 {
00035 double value;
00036
00037 DEBUG_ENTRY( "rt_continuum_shield_fcn()" );
00038
00039 value = -1.;
00040
00041 ASSERT( t.Emis().damp() > 0. );
00042
00043 if( rt.nLineContShield == LINE_CONT_SHIELD_PESC )
00044 {
00045
00046 if( t.Emis().iRedisFun() == ipPRD )
00047 {
00048 value = esc_PRD_1side(t.Emis().TauCon(),t.Emis().damp());
00049 }
00050 else if( t.Emis().iRedisFun() == ipCRD )
00051 {
00052 value = esca0k2(t.Emis().TauCon());
00053 }
00054 else if( t.Emis().iRedisFun() == ipCRDW )
00055 {
00056 value = esc_CRDwing_1side(t.Emis().TauCon(),t.Emis().damp());
00057 }
00058 else if( t.Emis().iRedisFun() == ipLY_A )
00059 {
00060 value = esc_PRD_1side(t.Emis().TauCon(),t.Emis().damp());
00061 }
00062 else
00063 TotalInsanity();
00064 }
00065 else if( rt.nLineContShield == LINE_CONT_SHIELD_FEDERMAN )
00066 {
00067
00068 double core, wings;
00069
00070
00071
00072
00073
00074 if( t.Emis().TauCon() < 2. )
00075 {
00076 core = sexp( t.Emis().TauCon() * 0.66666 );
00077 }
00078 else if( t.Emis().TauCon() < 10. )
00079 {
00080 core = 0.638 * pow(t.Emis().TauCon(),(realnum)-1.25f );
00081 }
00082 else if( t.Emis().TauCon() < 100. )
00083 {
00084 core = 0.505 * pow(t.Emis().TauCon(),(realnum)-1.15f );
00085 }
00086 else
00087 {
00088 core = 0.344 * pow(t.Emis().TauCon(),(realnum)-1.0667f );
00089 }
00090
00091
00092 wings = 0.;
00093 if( t.Emis().TauCon() > 1.f && t.Emis().damp()>0. )
00094 {
00095
00096 double t1 = 3.02*pow(t.Emis().damp()*1e3,-0.064 );
00097 double u1 = sqrt(t.Emis().TauCon()*t.Emis().damp() )/SDIV(t1);
00098 wings = t.Emis().damp()/SDIV(t1)/sqrt( 0.78540 + POW2(u1) );
00099
00100
00101
00102
00103 if( t.Emis().TauCon()>1e7 )
00104 wings *= pow( t.Emis().TauCon()/1e7,-1.1 );
00105 }
00106 value = core + wings;
00107
00108
00109
00110
00111 if( t.Emis().TauCon()>0. )
00112 value = MIN2(1., value );
00113 }
00114 else if( rt.nLineContShield == LINE_CONT_SHIELD_FERLAND )
00115 {
00116
00117 value = conpmp( t );
00118 }
00119 else if( rt.nLineContShield == 0 )
00120 {
00121
00122 value = 1.;
00123 }
00124 else
00125 {
00126 TotalInsanity();
00127 }
00128
00129
00130
00131 ASSERT( value>=0 && (value<=1.||t.Emis().TauCon()<0.) );
00132
00133 return value;
00134 }
00135
00138 static const double BREAK = 3.;
00139
00140 STATIC inline double FITTED( double t )
00141 {
00142 return (0.98925439 + 0.084594094*t)/(1. + t*(0.64794212 + t*0.44743976));
00143 }
00144
00145
00146 STATIC double conpmp(const TransitionProxy& t)
00147 {
00148 double a0,
00149 conpmp_v,
00150 tau,
00151 yinc1,
00152 yinc2;
00153
00154 DEBUG_ENTRY( "conpmp()" );
00155
00156
00157
00158 tau = t.Emis().TauCon();
00159
00160 if( tau <= 10. )
00161 {
00162
00163 conpmp_v = FITTED(tau);
00164 }
00165 else if( tau > 1e6 )
00166 {
00167
00168
00169 conpmp_v = 0.;
00170 }
00171 else
00172 {
00173 my_Integrand func;
00174 func.PumpDamp = t.Emis().damp();
00175 func.PumpTau = tau;
00176 Integrator<my_Integrand,Gaussian32> opfun;
00177
00178 yinc1 = opfun.sum( 0., BREAK, func );
00179 yinc2 = opfun.sum( BREAK, 100., func );
00180
00181 a0 = 0.886227*(1. + func.PumpDamp);
00182 conpmp_v = (yinc1 + yinc2)/a0;
00183 }
00184
00185
00186
00187
00188
00189 return conpmp_v;
00190 }