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