RefractionExtinction.cpp   RefractionExtinction.cpp 
skipping to change at line 109 skipping to change at line 109
invertPostTransfoMat=m.inverse(); invertPostTransfoMat=m.inverse();
postTransfoMatf.set(m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13], m[14], m[15]); postTransfoMatf.set(m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13], m[14], m[15]);
invertPostTransfoMatf.set(invertPostTransfoMat[0], invertPostTransfo Mat[1], invertPostTransfoMat[2], invertPostTransfoMat[3], invertPostTransfoMatf.set(invertPostTransfoMat[0], invertPostTransfo Mat[1], invertPostTransfoMat[2], invertPostTransfoMat[3],
invertPostTransfoMa t[4], invertPostTransfoMat[5], invertPostTransfoMat[6], invertPostTransfoMa t[7], invertPostTransfoMa t[4], invertPostTransfoMat[5], invertPostTransfoMat[6], invertPostTransfoMa t[7],
invertPostTransfoMa t[8], invertPostTransfoMat[9], invertPostTransfoMat[10], invertPostTransfoM at[11], invertPostTransfoMa t[8], invertPostTransfoMat[9], invertPostTransfoMat[10], invertPostTransfoM at[11],
invertPostTransfoMa t[12], invertPostTransfoMat[13], invertPostTransfoMat[14], invertPostTransf oMat[15]); invertPostTransfoMa t[12], invertPostTransfoMat[13], invertPostTransfoMat[14], invertPostTransf oMat[15]);
} }
void Refraction::updatePrecomputed() void Refraction::updatePrecomputed()
{ {
press_temp_corr_Bennett=pressure/1010.f * 283.f/(273.f+temperature) press_temp_corr=pressure/1010.f * 283.f/(273.f+temperature) / 60.f;
/ 60.f;
press_temp_corr_Saemundson=1.02f*press_temp_corr_Bennett;
} }
void Refraction::innerRefractionForward(Vec3d& altAzPos) const void Refraction::innerRefractionForward(Vec3d& altAzPos) const
{ {
const double length = altAzPos.length(); const double length = altAzPos.length();
double geom_alt_deg=180./M_PI*std::asin(altAzPos[2]/length); const double sinGeo = altAzPos[2]/length;
double geom_alt_rad = std::asin(sinGeo);
double geom_alt_deg = 180./M_PI*geom_alt_rad;
if (geom_alt_deg > MIN_GEO_ALTITUDE_DEG) if (geom_alt_deg > MIN_GEO_ALTITUDE_DEG)
{ {
// refraction from Saemundsson, S&T1986 p70 / in Meeus, Astr .Alg. // refraction from Saemundsson, S&T1986 p70 / in Meeus, Astr .Alg.
float r=press_temp_corr_Saemundson / std::tan((geom_alt_deg+ 10.3f/(geom_alt_deg+5.11f))*M_PI/180.f) + 0.0019279f; double r=press_temp_corr * ( 1.02 / std::tan((geom_alt_deg+1 0.3/(geom_alt_deg+5.11))*M_PI/180.) + 0.0019279);
geom_alt_deg += r; geom_alt_deg += r;
if (geom_alt_deg > 90.) if (geom_alt_deg > 90.)
geom_alt_deg=90.; geom_alt_deg=90.;
altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.)*length;
} }
else if(geom_alt_deg>MIN_GEO_ALTITUDE_DEG-TRANSITION_WIDTH_GEO_DEG) else if(geom_alt_deg>MIN_GEO_ALTITUDE_DEG-TRANSITION_WIDTH_GEO_DEG)
{ {
// Avoids the jump below -5 by interpolating linearly betwee n MIN_GEO_ALTITUDE_DEG and bottom of transition zone // Avoids the jump below -5 by interpolating linearly betwee n MIN_GEO_ALTITUDE_DEG and bottom of transition zone
float r_m5=press_temp_corr_Saemundson / std::tan((MIN_GEO_AL TITUDE_DEG+10.3f/(MIN_GEO_ALTITUDE_DEG+5.11f))*M_PI/180.f) + 0.0019279f; float r_m5=press_temp_corr * ( 1.02f / std::tan((MIN_GEO_ALT ITUDE_DEG+10.3f/(MIN_GEO_ALTITUDE_DEG+5.11f))*M_PI/180.f) + 0.0019279f);
geom_alt_deg += r_m5*(geom_alt_deg-(MIN_GEO_ALTITUDE_DEG-TRA NSITION_WIDTH_GEO_DEG))/TRANSITION_WIDTH_GEO_DEG; geom_alt_deg += r_m5*(geom_alt_deg-(MIN_GEO_ALTITUDE_DEG-TRA NSITION_WIDTH_GEO_DEG))/TRANSITION_WIDTH_GEO_DEG;
altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.)*length;
} }
else return;
// At this point we have corrected geometric altitude. Note that if
we just change altAzPos[2], we would change vector length, so this would ch
ange our angles.
// We have to shorten X,Y components of the vector as well by the ch
ange in cosines of altitude, or (sqrt(1-sin(alt))
const double refr_alt_rad=geom_alt_deg*M_PI/180.;
const double sinRef=std::sin(refr_alt_rad);
const double shortenxy=std::sqrt((1.-sinRef*sinRef)/(1.-sinGeo*sinGe
o)); // we need double's mantissa length here, sorry!
altAzPos[0]*=shortenxy;
altAzPos[1]*=shortenxy;
altAzPos[2]=sinRef*length;
} }
void Refraction::innerRefractionBackward(Vec3d& altAzPos) const void Refraction::innerRefractionBackward(Vec3d& altAzPos) const
{ {
// going from observed position/magnitude to geometrical position an d atmosphere-free mag. // going from observed position/magnitude to geometrical position an d atmosphere-free mag.
const double length = altAzPos.length(); const double length = altAzPos.length();
float obs_alt_deg=180./M_PI*std::asin(altAzPos[2]/length); const double sinObs = altAzPos[2]/length;
if (obs_alt_deg > 0.22879) double obs_alt_deg=180./M_PI*std::asin(sinObs);
if (obs_alt_deg > 0.22879f)
{ {
// refraction from Bennett, in Meeus, Astr.Alg. // refraction from Bennett, in Meeus, Astr.Alg.
float r=press_temp_corr_Bennett / std::tan((obs_alt_deg+7.31 f/(obs_alt_deg+4.4f))*M_PI/180.) + 0.0013515f; double r=press_temp_corr * (1. / std::tan((obs_alt_deg+7.31/ (obs_alt_deg+4.4))*M_PI/180.) + 0.0013515);
obs_alt_deg -= r; obs_alt_deg -= r;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.f)*length;
} }
else if (obs_alt_deg > MIN_APP_ALTITUDE_DEG) else if (obs_alt_deg > MIN_APP_ALTITUDE_DEG)
{ {
// backward refraction from polynomial fit against Saemundso n[-5...-0.3] // backward refraction from polynomial fit against Saemundso n[-5...-0.3]
float r=(((((0.0444f*obs_alt_deg+.7662f)*obs_alt_deg+4.9746f )*obs_alt_deg+13.599f)*obs_alt_deg+8.052f)*obs_alt_deg-11.308f)*obs_alt_deg +34.341f; float r=(((((0.0444f*obs_alt_deg+.7662f)*obs_alt_deg+4.9746f )*obs_alt_deg+13.599f)*obs_alt_deg+8.052f)*obs_alt_deg-11.308f)*obs_alt_deg +34.341f;
obs_alt_deg -= press_temp_corr_Bennett*r; obs_alt_deg -= press_temp_corr*r;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length;
} }
else if (obs_alt_deg > MIN_APP_ALTITUDE_DEG-TRANSITION_WIDTH_APP_DEG ) else if (obs_alt_deg > MIN_APP_ALTITUDE_DEG-TRANSITION_WIDTH_APP_DEG )
{ {
// Compute top value from polynome, apply linear interpolati on // Compute top value from polynome, apply linear interpolati on
static const float r_min=(((((0.0444f*MIN_APP_ALTITUDE_DEG+. 7662f)*MIN_APP_ALTITUDE_DEG static const float r_min=(((((0.0444f*MIN_APP_ALTITUDE_DEG+. 7662f)*MIN_APP_ALTITUDE_DEG
+4.9746f)*MIN_APP_ALTITUDE_DEG+13.599f)*MIN_ APP_ALTITUDE_DEG +4.9746f)*MIN_APP_ALTITUDE_DEG+13.599f)*MIN_ APP_ALTITUDE_DEG
+8.052f)*MIN_APP_ALTITUDE_DEG-11.308f)*MIN_APP _ALTITUDE_DEG+34.341f; +8.052f)*MIN_APP_ALTITUDE_DEG-11.308f)*MIN_APP _ALTITUDE_DEG+34.341f;
obs_alt_deg -= r_min*press_temp_corr_Bennett*(obs_alt_deg-(M obs_alt_deg -= r_min*press_temp_corr*(obs_alt_deg-(MIN_APP_A
IN_APP_ALTITUDE_DEG-TRANSITION_WIDTH_APP_DEG))/TRANSITION_WIDTH_APP_DEG; LTITUDE_DEG-TRANSITION_WIDTH_APP_DEG))/TRANSITION_WIDTH_APP_DEG;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length;
} }
else return;
// At this point we have corrected observed altitude. Note that if w
e just change altAzPos[2], we would change vector length, so this would cha
nge our angles.
// We have to make X,Y components of the vector a bit longer as well
by the change in cosines of altitude, or (sqrt(1-sin(alt))
const double geo_alt_rad=obs_alt_deg*M_PI/180.;
const double sinGeo=std::sin(geo_alt_rad);
const double longerxy=std::sqrt((1.-sinGeo*sinGeo)/(1.-sinObs*sinObs
));
altAzPos[0]*=longerxy;
altAzPos[1]*=longerxy;
altAzPos[2]=sinGeo*length;
} }
void Refraction::forward(Vec3d& altAzPos) const void Refraction::forward(Vec3d& altAzPos) const
{ {
altAzPos.transfo4d(preTransfoMat); altAzPos.transfo4d(preTransfoMat);
innerRefractionForward(altAzPos); innerRefractionForward(altAzPos);
altAzPos.transfo4d(postTransfoMat); altAzPos.transfo4d(postTransfoMat);
} }
//Bennett's formula is not a strict inverse of Saemundsson's. There is a no table discrepancy (alt!=backward(forward(alt))) for //Bennett's formula is not a strict inverse of Saemundsson's. There is a no table discrepancy (alt!=backward(forward(alt))) for
 End of changes. 13 change blocks. 
17 lines changed or deleted 43 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/