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