RefractionExtinction.cpp   RefractionExtinction.cpp 
skipping to change at line 115 skipping to change at line 115
} }
void Refraction::updatePrecomputed() void Refraction::updatePrecomputed()
{ {
press_temp_corr=pressure/1010.f * 283.f/(273.f+temperature) / 60.f; press_temp_corr=pressure/1010.f * 283.f/(273.f+temperature) / 60.f;
} }
void Refraction::innerRefractionForward(Vec3d& altAzPos) const void Refraction::innerRefractionForward(Vec3d& altAzPos) const
{ {
const double length = altAzPos.length(); const double length = altAzPos.length();
if (length==0.0)
{
// Under some circumstances there are zero coordinates. Just
leave them alone.
//qDebug() << "Refraction::innerRefractionForward(): Zero ve
ctor detected - Continue with zero vector.";
return;
}
Q_ASSERT(length>0.0);
const double sinGeo = altAzPos[2]/length; const double sinGeo = altAzPos[2]/length;
Q_ASSERT(fabs(sinGeo)<=1.0);
double geom_alt_rad = std::asin(sinGeo); double geom_alt_rad = std::asin(sinGeo);
double geom_alt_deg = 180./M_PI*geom_alt_rad; float 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.
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); float r=press_temp_corr * ( 1.02f / std::tan((geom_alt_deg+1 0.3f/(geom_alt_deg+5.11f))*M_PI/180.f) + 0.0019279f);
geom_alt_deg += r; geom_alt_deg += r;
if (geom_alt_deg > 90.) if (geom_alt_deg > 90.f)
geom_alt_deg=90.; geom_alt_deg=90.f;
} }
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 * ( 1.02f / std::tan((MIN_GEO_ALT ITUDE_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;
} }
else return; 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. // 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)) // 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 refr_alt_rad=geom_alt_deg*M_PI/180.;
const double sinRef=std::sin(refr_alt_rad); 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! const double shortenxy=((fabs(sinGeo)>=1.0) ? 1.0 :
std::sqrt((1.-sinRef*sinRef)/(1.-sinGeo*sinGeo))); /
/ we need double's mantissa length here, sorry!
altAzPos[0]*=shortenxy; altAzPos[0]*=shortenxy;
altAzPos[1]*=shortenxy; altAzPos[1]*=shortenxy;
altAzPos[2]=sinRef*length; altAzPos[2]=sinRef*length;
} }
// going from observed position to geometrical position.
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.
const double length = altAzPos.length(); const double length = altAzPos.length();
if (length==0.0)
{
// Under some circumstances there are zero coordinates. Just
leave them alone.
//qDebug() << "Refraction::innerRefractionBackward(): Zero v
ector detected - Continue with zero vector.";
return;
}
Q_ASSERT(length>0.0);
const double sinObs = altAzPos[2]/length; const double sinObs = altAzPos[2]/length;
double obs_alt_deg=180./M_PI*std::asin(sinObs); Q_ASSERT(fabs(sinObs)<=1.0);
float obs_alt_deg=180./M_PI*std::asin(sinObs);
if (obs_alt_deg > 0.22879f) if (obs_alt_deg > 0.22879f)
{ {
// refraction from Bennett, in Meeus, Astr.Alg. // refraction from Bennett, in Meeus, Astr.Alg.
double r=press_temp_corr * (1. / std::tan((obs_alt_deg+7.31/ (obs_alt_deg+4.4))*M_PI/180.) + 0.0013515); float r=press_temp_corr * (1.f / std::tan((obs_alt_deg+7.31f /(obs_alt_deg+4.4f))*M_PI/180.f) + 0.0013515f);
obs_alt_deg -= r; obs_alt_deg -= r;
} }
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*r; obs_alt_deg -= press_temp_corr*r;
} }
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 )
{ {
skipping to change at line 179 skipping to change at line 197
+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*(obs_alt_deg-(MIN_APP_A LTITUDE_DEG-TRANSITION_WIDTH_APP_DEG))/TRANSITION_WIDTH_APP_DEG; obs_alt_deg -= r_min*press_temp_corr*(obs_alt_deg-(MIN_APP_A LTITUDE_DEG-TRANSITION_WIDTH_APP_DEG))/TRANSITION_WIDTH_APP_DEG;
} }
else return; 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. // 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)) // 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 geo_alt_rad=obs_alt_deg*M_PI/180.;
const double sinGeo=std::sin(geo_alt_rad); const double sinGeo=std::sin(geo_alt_rad);
const double longerxy=std::sqrt((1.-sinGeo*sinGeo)/(1.-sinObs*sinObs const double longerxy=((fabs(sinObs)>=1.0) ? 1.0 :
)); std::sqrt((1.-sinGeo*sinGeo)/(1.-sinObs*sinObs)));
altAzPos[0]*=longerxy; altAzPos[0]*=longerxy;
altAzPos[1]*=longerxy; altAzPos[1]*=longerxy;
altAzPos[2]=sinGeo*length; 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);
 End of changes. 13 change blocks. 
12 lines changed or deleted 34 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/