RefractionExtinction.cpp   RefractionExtinction.cpp 
skipping to change at line 23 skipping to change at line 23
* GNU General Public License for more details. * GNU General Public License for more details.
* *
* You should have received a copy of the GNU General Public License * You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software * along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA. * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
* *
* Refraction and extinction computations. * Refraction and extinction computations.
* Principal implementation: 2010-03-23 GZ=Georg Zotti, Georg.Zotti@univie. ac.at * Principal implementation: 2010-03-23 GZ=Georg Zotti, Georg.Zotti@univie. ac.at
*/ */
#include <qsettings.h>
#include "StelApp.hpp" #include "StelApp.hpp"
#include "RefractionExtinction.hpp" #include "RefractionExtinction.hpp"
// To be decided: The following should be either 0 or 40 (or 42? ;-) Extinction::Extinction() : ext_coeff(50), undergroundExtinctionMode(Undergr
float Extinction::SUBHORIZONTAL_AIRMASS=0.0f; oundExtinctionMirror)
Extinction::Extinction()
{
QSettings* conf = StelApp::getInstance().getSettings();
SUBHORIZONTAL_AIRMASS = (conf->value("astro/flag_extinction_below_horiz
on", true).toBool()? 42.0f : 0.0f);
ext_coeff=conf->value("landscape/atmospheric_extinction_coefficient", 0
.2f).toFloat();
}
// altAzPos is the NORMALIZED (!!!) star position vector AFTER REFRACTION,
and its z component sin(altitude).
void Extinction::forward(const Vec3d *altAzPos, float *mag, const int num)
const
{
for (int i=0; i<num; ++i) mag[i] += airmass(altAzPos[i][2], true) *
ext_coeff;
}
void Extinction::forward(const Vec3f *altAzPos, float *mag, const int num)
const
{
for (int i=0; i<num; ++i) mag[i] += airmass(altAzPos[i][2], true) *
ext_coeff;
}
// If only sin(altitude) is available:
void Extinction::forward(const double *sinAlt, float *mag, const int num) c
onst
{
for (int i=0; i<num; ++i) mag[i] += airmass(sinAlt[i], true) * ext_c
oeff;
}
void Extinction::forward(const float *sinAlt, float *mag, const int num) co
nst
{
for (int i=0; i<num; ++i) mag[i] += airmass(sinAlt[i], true) * ext_c
oeff;
}
void Extinction::forward(const double *sinAlt, float *mag) const
{
*mag += airmass(*sinAlt, true) * ext_coeff;
}
void Extinction::forward(const float *sinAlt, float *mag) const
{
*mag += airmass(*sinAlt, true) * ext_coeff;
}
// from observed magnitude in apparent (observed) altitude to atmosphere-fr
ee mag, still in apparent, refracted altitude.
void Extinction::backward(const Vec3d *altAzPos, float *mag, const int num)
const
{
for (int i=0; i<num; ++i) mag[i] -= airmass(altAzPos[i][2], true) *
ext_coeff;
}
void Extinction::backward(const Vec3f *altAzPos, float *mag, const int num)
const
{
for (int i=0; i<num; ++i) mag[i] -= airmass(altAzPos[i][2], true) *
ext_coeff;
}
// If only sin(altitude) is available:
void Extinction::backward(const double *sinAlt, float *mag, const int num)
const
{ {
for (int i=0; i<num; ++i) mag[i] -= airmass(sinAlt[i], true) * ext_c
oeff;
}
void Extinction::backward(const float *sinAlt, float *mag, const int num) c
onst
{
for (int i=0; i<num; ++i) mag[i] -= airmass(sinAlt[i], true) * ext_c
oeff;
} }
// airmass computation for cosine of zenith angle z // airmass computation for cosine of zenith angle z
float Extinction::airmass(const float cosZ, const bool apparent_z) const float Extinction::airmass(float cosZ, const bool apparent_z) const
{ {
if (cosZ<-0.035f) // about -2 degrees. Here, RozenbergZ>574 and clim bs fast! if (cosZ<-0.035f) // about -2 degrees. Here, RozenbergZ>574 and clim bs fast!
return Extinction::SUBHORIZONTAL_AIRMASS; // Safety: 0 or 40 for {
below -2 degrees. switch (undergroundExtinctionMode)
{
case UndergroundExtinctionZero:
return 0.f;
case UndergroundExtinctionMax:
return 42.f;
case UndergroundExtinctionMirror:
cosZ = std::min(1.f, -0.035f - (cosZ+0.035f)
);
}
}
float X;
if (apparent_z) if (apparent_z)
{ {
// Rozenberg 1966, reported by Schaefer (1993-2000). // Rozenberg 1966, reported by Schaefer (1993-2000).
X=1.0f/(cosZ+0.025f*std::exp(-11.f*cosZ)); return 1.0f/(cosZ+0.025f*std::exp(-11.f*cosZ));
} }
else else
{ {
//Young 1994 //Young 1994
float nom=(1.002432f*cosZ+0.148386f)*cosZ+0.0096467f; const float nom=(1.002432f*cosZ+0.148386f)*cosZ+0.0096467f;
float denum=((cosZ+0.149864f)*cosZ+0.0102963f)*cosZ+0.000303 const float denum=((cosZ+0.149864f)*cosZ+0.0102963f)*cosZ+0.
978f; 000303978f;
X=nom/denum; return nom/denum;
} }
return X;
} }
/* ************************************************************************ ***************************** */ /* ************************************************************************ ***************************** */
// The following 4 are to be configured, the rest is derived. // The following 4 are to be configured, the rest is derived.
// Recommendations: -4.9/-4.3/0.1/0.1: sharp but continuous transition, no effects below -5. // Recommendations: -4.9/-4.3/0.1/0.1: sharp but continuous transition, no effects below -5.
// -4.3/-4.3/0.7/0.7: sharp but continuous transition, no effects below -5. Maybe better for picking? // -4.3/-4.3/0.7/0.7: sharp but continuous transition, no effects below -5. Maybe better for picking?
// -3/-3/2/2: "strange" zone 2 degrees wide. Both formulae are close near -3. Actually, refraction formulae dubious below 0 // -3/-3/2/2: "strange" zone 2 degrees wide. Both formulae are close near -3. Actually, refraction formulae dubious below 0
// 0/0/1/1: "strange zone 1 degree wide just below horizo n, no effect below -1. Actually, refraction formulae dubious below 0! But i t looks stupid, see the sun! // 0/0/1/1: "strange zone 1 degree wide just below horizo n, no effect below -1. Actually, refraction formulae dubious below 0! But i t looks stupid, see the sun!
// Optimum:-3.54/-3.21783/1.46/1.78217. Here forward/backw ard are almost perfect inverses (-->good picking!), and nothing happens bel ow -5 degrees. // Optimum:-3.54/-3.21783/1.46/1.78217. Here forward/backw ard are almost perfect inverses (-->good picking!), and nothing happens bel ow -5 degrees.
// This must be -5 or higher. // This must be -5 or higher.
const double Refraction::MIN_GEO_ALTITUDE_DEG=-3.54; static const float MIN_GEO_ALTITUDE_DEG=-3.54f;
// this must be -4.3 or higher. -3.21783 is an optimal value to fit against forward refraction! // this must be -4.3 or higher. -3.21783 is an optimal value to fit against forward refraction!
const double Refraction::MIN_APP_ALTITUDE_DEG=-3.21783; static const float MIN_APP_ALTITUDE_DEG=-3.21783f;
// this must be positive. Transition zone goes that far below the values ju st specified. // this must be positive. Transition zone goes that far below the values ju st specified.
const double Refraction::TRANSITION_WIDTH_GEO_DEG=1.46; static const float TRANSITION_WIDTH_GEO_DEG=1.46f;
const double Refraction::TRANSITION_WIDTH_APP_DEG=1.78217; static const float TRANSITION_WIDTH_APP_DEG=1.78217f;
const double Refraction::MIN_GEO_ALTITUDE_RAD=Refraction::MIN_GEO_ALTITUDE_ static const float MIN_GEO_ALTITUDE_SIN=std::sin(MIN_GEO_ALTITUDE_DEG*M_PI/
DEG*M_PI/180.0; 180.f);
const double Refraction::MIN_GEO_ALTITUDE_SIN=std::sin(Refraction::MIN_GEO_ static const float MIN_APP_ALTITUDE_SIN=std::sin(MIN_APP_ALTITUDE_DEG*M_PI/
ALTITUDE_RAD); 180.f);
const double Refraction::MIN_APP_ALTITUDE_RAD=Refraction::MIN_APP_ALTITUDE_
DEG*M_PI/180.0;
const double Refraction::MIN_APP_ALTITUDE_SIN=std::sin(Refraction::MIN_APP_
ALTITUDE_RAD);
const float Refraction::MIN_GEO_ALTITUDE_DEG_F=(float)Refraction::MIN_GEO_A
LTITUDE_DEG;
const float Refraction::MIN_GEO_ALTITUDE_RAD_F=(float)Refraction::MIN_GEO_A
LTITUDE_RAD;
const float Refraction::MIN_GEO_ALTITUDE_SIN_F=(float)Refraction::MIN_GEO_A
LTITUDE_SIN;
const float Refraction::MIN_APP_ALTITUDE_DEG_F=(float)Refraction::MIN_APP_A
LTITUDE_DEG;
const float Refraction::MIN_APP_ALTITUDE_RAD_F=(float)Refraction::MIN_APP_A
LTITUDE_RAD;
const float Refraction::MIN_APP_ALTITUDE_SIN_F=(float)Refraction::MIN_APP_A
LTITUDE_SIN;
const double Refraction::TRANSITION_WIDTH_GEO_DEG_F=(float)Refraction::TRAN
SITION_WIDTH_GEO_DEG;
const double Refraction::TRANSITION_WIDTH_APP_DEG_F=(float)Refraction::TRAN
SITION_WIDTH_APP_DEG;
Refraction::Refraction() : //pressure(1013.f), temperature(10.f), Refraction::Refraction() : pressure(1013.f), temperature(10.f),
preTransfoMat(Mat4d::identity()), invertPreTransfoMat(Mat4d::identit y()), preTransfoMatf(Mat4f::identity()), invertPreTransfoMatf(Mat4f::identi ty()), preTransfoMat(Mat4d::identity()), invertPreTransfoMat(Mat4d::identit y()), preTransfoMatf(Mat4f::identity()), invertPreTransfoMatf(Mat4f::identi ty()),
postTransfoMat(Mat4d::identity()), invertPostTransfoMat(Mat4d::ident ity()), postTransfoMatf(Mat4f::identity()), invertPostTransfoMatf(Mat4f::id entity()) postTransfoMat(Mat4d::identity()), invertPostTransfoMat(Mat4d::ident ity()), postTransfoMatf(Mat4f::identity()), invertPostTransfoMatf(Mat4f::id entity())
{ {
QSettings* conf = StelApp::getInstance().getSettings();
pressure=conf->value("landscape/pressure_mbar", 1013.0f).toFloat();
temperature=conf->value("landscape/temperature_C", 15.0f).toFloat();
updatePrecomputed(); updatePrecomputed();
} }
void Refraction::setPreTransfoMat(const Mat4d& m) void Refraction::setPreTransfoMat(const Mat4d& m)
{ {
preTransfoMat=m; preTransfoMat=m;
invertPreTransfoMat=m.inverse(); invertPreTransfoMat=m.inverse();
preTransfoMatf.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]); preTransfoMatf.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]);
invertPreTransfoMatf.set(invertPreTransfoMat[0], invertPreTransfoMat [1], invertPreTransfoMat[2], invertPreTransfoMat[3], invertPreTransfoMatf.set(invertPreTransfoMat[0], invertPreTransfoMat [1], invertPreTransfoMat[2], invertPreTransfoMat[3],
invertPreTransfoMat [4], invertPreTransfoMat[5], invertPreTransfoMat[6], invertPreTransfoMat[7] , invertPreTransfoMat [4], invertPreTransfoMat[5], invertPreTransfoMat[6], invertPreTransfoMat[7] ,
skipping to change at line 170 skipping to change at line 113
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) / 60.f; press_temp_corr_Bennett=pressure/1010.f * 283.f/(273.f+temperature) / 60.f;
press_temp_corr_Saemundson=1.02f*press_temp_corr_Bennett; press_temp_corr_Saemundson=1.02f*press_temp_corr_Bennett;
} }
void Refraction::forward(Vec3d& altAzPos) const void Refraction::innerRefractionForward(Vec3d& altAzPos) const
{ {
altAzPos.transfo4d(preTransfoMat);
const double length = altAzPos.length(); const double length = altAzPos.length();
double geom_alt_deg=180./M_PI*std::asin(altAzPos[2]/length); double geom_alt_deg=180./M_PI*std::asin(altAzPos[2]/length);
if (geom_alt_deg > Refraction::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; float r=press_temp_corr_Saemundson / std::tan((geom_alt_deg+ 10.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.) geom_alt_deg=90.; // SAFETY if (geom_alt_deg > 90.)
geom_alt_deg=90.;
altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.)*length; altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.)*length;
} }
else if(geom_alt_deg>Refraction::MIN_GEO_ALTITUDE_DEG-Refraction::TR ANSITION_WIDTH_GEO_DEG) else if(geom_alt_deg>MIN_GEO_ALTITUDE_DEG-TRANSITION_WIDTH_GEO_DEG)
{ {
// Avoids the jump near -5 by interpolating linearly between // Avoids the jump below -5 by interpolating linearly betwee
MIN_GEO_ALTITUDE_DEG and bottom of transition zone n MIN_GEO_ALTITUDE_DEG and bottom of transition zone
float r_min=press_temp_corr_Saemundson / std::tan((Refractio float r_m5=press_temp_corr_Saemundson / std::tan((MIN_GEO_AL
n::MIN_GEO_ALTITUDE_DEG+10.3f/(Refraction::MIN_GEO_ALTITUDE_DEG+5.11f))*M_P TITUDE_DEG+10.3f/(MIN_GEO_ALTITUDE_DEG+5.11f))*M_PI/180.f) + 0.0019279f;
I/180.f) + 0.0019279f; geom_alt_deg += r_m5*(geom_alt_deg-(MIN_GEO_ALTITUDE_DEG-TRA
geom_alt_deg += r_min*(geom_alt_deg-(Refraction::MIN_GEO_ALT NSITION_WIDTH_GEO_DEG))/TRANSITION_WIDTH_GEO_DEG;
ITUDE_DEG-Refraction::TRANSITION_WIDTH_GEO_DEG))/Refraction::TRANSITION_WID
TH_GEO_DEG;
altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.)*length; altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.)*length;
} }
altAzPos.transfo4d(postTransfoMat);
} }
//Bennett's formula is not a strict inverse of Saemundsson's. There is a no void Refraction::innerRefractionBackward(Vec3d& altAzPos) const
table discrepancy (alt!=backward(forward(alt))) for
// geometric altitudes <-.3deg. This is not a problem in real life, but if
a user switches off landscape, click-identify may pose a problem.
// Below this altitude, we therefore use a polynomial fit that should repre
sent a close inverse of Saemundsson.
void Refraction::backward(Vec3d& altAzPos) const
{ {
altAzPos.transfo4d(invertPostTransfoMat); // going from observed position/magnitude to geometrical position an
// going from apparent (observed) position to geometrical position. d atmosphere-free mag.
const double length = altAzPos.length(); const double length = altAzPos.length();
double obs_alt_deg=180./M_PI*std::asin(altAzPos[2]/length); float obs_alt_deg=180./M_PI*std::asin(altAzPos[2]/length);
if (obs_alt_deg > 0.22879) if (obs_alt_deg > 0.22879)
{ {
// refraction directly 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 float r=press_temp_corr_Bennett / std::tan((obs_alt_deg+7.31
/(obs_alt_deg+4.4f))*M_PI/180.f) + 0.0013515f; f/(obs_alt_deg+4.4f))*M_PI/180.) + 0.0013515f;
obs_alt_deg -= r; obs_alt_deg -= r;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length; altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.f)*length;
} }
else if (obs_alt_deg > Refraction::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.0444*obs_alt_deg+.7662)*obs_alt_deg+4.9746)*o bs_alt_deg+13.599)*obs_alt_deg+8.052)*obs_alt_deg-11.308)*obs_alt_deg+34.34 1; 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_Bennett*r;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length; altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length;
} }
else if (obs_alt_deg > Refraction::MIN_APP_ALTITUDE_DEG-Refraction:: TRANSITION_WIDTH_APP_DEG) else if (obs_alt_deg > MIN_APP_ALTITUDE_DEG-TRANSITION_WIDTH_APP_DEG )
{ {
// use polynomial fit for topmost value, linear interpolatio // Compute top value from polynome, apply linear interpolati
n inside transition zone on
float r_min=(((((0.0444*Refraction::MIN_APP_ALTITUDE_DEG+.76 static const float r_min=(((((0.0444f*MIN_APP_ALTITUDE_DEG+.
62)*Refraction::MIN_APP_ALTITUDE_DEG 7662f)*MIN_APP_ALTITUDE_DEG
+4.9746)*Refraction::MIN_APP_ALTITUDE_DEG+13 +4.9746f)*MIN_APP_ALTITUDE_DEG+13.599f)*MIN_
.599)*Refraction::MIN_APP_ALTITUDE_DEG APP_ALTITUDE_DEG
+8.052)*Refraction::MIN_APP_ALTITUDE_DEG-11.30 +8.052f)*MIN_APP_ALTITUDE_DEG-11.308f)*MIN_APP
8)*Refraction::MIN_APP_ALTITUDE_DEG+34.341; _ALTITUDE_DEG+34.341f;
r_min*=press_temp_corr_Bennett;
obs_alt_deg -= r_min*(obs_alt_deg-(Refraction::MIN_APP_ALTIT obs_alt_deg -= r_min*press_temp_corr_Bennett*(obs_alt_deg-(M
UDE_DEG-Refraction::TRANSITION_WIDTH_APP_DEG))/Refraction::TRANSITION_WIDTH IN_APP_ALTITUDE_DEG-TRANSITION_WIDTH_APP_DEG))/TRANSITION_WIDTH_APP_DEG;
_APP_DEG;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length; altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length;
} }
}
void Refraction::forward(Vec3d& altAzPos) const
{
altAzPos.transfo4d(preTransfoMat);
innerRefractionForward(altAzPos);
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
// geometric altitudes <-.3deg. This is not a problem in real life, but if
a user switches off landscape, click-identify may pose a problem.
// Below this altitude, we therefore use a polynomial fit that should repre
sent a close inverse of Saemundsson.
void Refraction::backward(Vec3d& altAzPos) const
{
altAzPos.transfo4d(invertPostTransfoMat);
innerRefractionBackward(altAzPos);
altAzPos.transfo4d(invertPreTransfoMat); altAzPos.transfo4d(invertPreTransfoMat);
} }
void Refraction::forward(Vec3f& altAzPos) const void Refraction::forward(Vec3f& altAzPos) const
{ {
// Using doubles internally to avoid jitter. Vec3d vf(altAzPos[0], altAzPos[1], altAzPos[2]);
// (This affects planet drawing - which is done using floats, vf.transfo4d(preTransfoMat);
// as doubles are unsupported/slow on GPUs) innerRefractionForward(vf);
Vec3d altAzPosD(altAzPos[0], altAzPos[1], altAzPos[2]); vf.transfo4d(postTransfoMat);
altAzPosD.transfo4d(preTransfoMat); altAzPos.set(vf[0], vf[1], vf[2]);
const double length = altAzPosD.length();
double geom_alt_deg=180./M_PI*std::asin(altAzPosD[2]/length);
if (geom_alt_deg > Refraction::MIN_GEO_ALTITUDE_DEG)
{
// refraction from Saemundsson, S&T1986 p70 / in Meeus, Astr
.Alg.
double r=press_temp_corr_Saemundson /
std::tan((geom_alt_deg+10.3/(geom_alt_deg+5.11))*M_P
I/180.) + 0.0019279;
geom_alt_deg += r;
if (geom_alt_deg > 90.) geom_alt_deg=90.; // SAFETY
altAzPosD[2]=std::sin(geom_alt_deg*M_PI/180.)*length;
}
else if(geom_alt_deg>Refraction::MIN_GEO_ALTITUDE_DEG-Refraction::TR
ANSITION_WIDTH_GEO_DEG)
{
// Avoids the jump below -5 by interpolating linearly betwee
n MIN_GEO_ALTITUDE_DEG and bottom of transition zone
double r_m5=press_temp_corr_Saemundson /
std::tan((Refraction::MIN_GEO_ALTITUDE_DEG+10.3/
(Refraction::MIN_GEO_ALTITUDE_DEG+5.11))
*M_PI/180.) + 0.0019279;
geom_alt_deg +=
r_m5*(geom_alt_deg-(Refraction::MIN_GEO_ALTITUDE_DEG
-
Refraction::TRANSITION_WIDTH_GEO
_DEG))
/Refraction::TRANSITION_WIDTH_GEO_DEG;
altAzPosD[2]=std::sin(geom_alt_deg*M_PI/180.)*length;
}
altAzPosD.transfo4d(postTransfoMat);
altAzPos = Vec3f(altAzPosD[0], altAzPosD[1], altAzPosD[2]);
} }
void Refraction::backward(Vec3f& altAzPos) const void Refraction::backward(Vec3f& altAzPos) const
{ {
altAzPos.transfo4d(invertPostTransfoMatf); altAzPos.transfo4d(invertPostTransfoMatf);
// going from observed position/magnitude to geometrical position an Vec3d vf(altAzPos[0], altAzPos[1], altAzPos[2]);
d atmosphere-free mag. innerRefractionBackward(vf);
const float length = altAzPos.length(); altAzPos.set(vf[0], vf[1], vf[2]);
float obs_alt_deg=180.f/M_PI*std::asin(altAzPos[2]/length);
if (obs_alt_deg > 0.22879f)
{
// 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.f) + 0.0013515f;
obs_alt_deg -= r;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.f)*length;
}
else if (obs_alt_deg > Refraction::MIN_APP_ALTITUDE_DEG_F)
{
// 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;
obs_alt_deg -= press_temp_corr_Bennett*r;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.f)*length;
}
else if (obs_alt_deg > Refraction::MIN_APP_ALTITUDE_DEG_F-Refraction
::TRANSITION_WIDTH_APP_DEG_F)
{
// Compute top value from polynome, apply linear interpolati
on
float r_min=(((((0.0444f*Refraction::MIN_APP_ALTITUDE_DEG_F+
.7662f)*Refraction::MIN_APP_ALTITUDE_DEG_F
+4.9746f)*Refraction::MIN_APP_ALTITUDE_DEG_F
+13.599f)*Refraction::MIN_APP_ALTITUDE_DEG_F
+8.052f)*Refraction::MIN_APP_ALTITUDE_DEG_F-11
.308f)*Refraction::MIN_APP_ALTITUDE_DEG_F+34.341f;
r_min*=press_temp_corr_Bennett;
obs_alt_deg -= r_min*(obs_alt_deg-(Refraction::MIN_APP_ALTIT
UDE_DEG_F-Refraction::TRANSITION_WIDTH_APP_DEG_F))/Refraction::TRANSITION_W
IDTH_APP_DEG_F;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.f)*length;
}
altAzPos.transfo4d(invertPreTransfoMatf); altAzPos.transfo4d(invertPreTransfoMatf);
} }
void Refraction::setPressure(float p) void Refraction::setPressure(float p)
{ {
pressure=p; pressure=p;
updatePrecomputed(); updatePrecomputed();
} }
void Refraction::setTemperature(float t) void Refraction::setTemperature(float t)
 End of changes. 34 change blocks. 
237 lines changed or deleted 89 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/