RefractionExtinction.cpp   RefractionExtinction.cpp 
skipping to change at line 20 skipping to change at line 20
* This program is distributed in the hope that it will be useful, * This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of * but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, U SA. * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, U SA.
* *
* Refraction and extinction computations. * Refraction and extinction computations.
* 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 "RefractionExtinction.hpp" #include "RefractionExtinction.hpp"
RefractionExtinction::RefractionExtinction() : pressure(1013.f), temperatur // To be decided: The following should be either 0 or 40 (or 42? ;-)
e(10.f), ext_coeff(0.20f) const float Extinction::SUBHORIZONTAL_AIRMASS=0.0f;
Extinction::Extinction() : ext_coeff(0.20f)
{ {
updatePrecomputed();
} }
void RefractionExtinction::updatePrecomputed() // altAzPos is the normalized star position vector AFTER REFRACTION, and i
ts z component sin(altitude).
void Extinction::forward(const Vec3d *altAzPos, float *mag, const int num)
const
{ {
press_temp_corr_Bennett=pressure/1010.f * 283.f/(273.f+temperature) for (int i=0; i<num; ++i) mag[i] += airmass(altAzPos[i][2], true) *
/ 60.f; ext_coeff;
press_temp_corr_Saemundson=1.02f*press_temp_corr_Bennett;
} }
void Extinction::forward(const Vec3f *altAzPos, float *mag, const int num)
void RefractionExtinction::forward(Vec3d* altAzPos, float* mag, int size) const
{ {
// Assuming altAzPos is the normalized star position vector, and its for (int i=0; i<num; ++i) mag[i] += airmass(altAzPos[i][2], true) *
z component sin(altitude). ext_coeff;
for (int i=0; i<size; ++i)
{
float geom_alt_deg=(180.f/M_PI)*std::asin(altAzPos[i][2]);
if (geom_alt_deg > -2.f)
{
// refraction from Saemundsson, S&T1986 p70 / in Mee
us, 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;
geom_alt_deg += r;
if (geom_alt_deg > 90.f)
geom_alt_deg=90.f; // SAFETY, SHOULD NOT BE
NECESSARY
altAzPos[i][2]=std::sin(geom_alt_deg*M_PI/180.f);
// now altitude is corrected, but object still too b
right.
// note that sin(h)=cos(z)...
mag[i] += airmass(altAzPos[i][2], true) * ext_coeff;
}
}
} }
// If only sin(altitude) is available:
void RefractionExtinction::backward(Vec3d* altAzPos, float* mag, int size) void Extinction::forward(const double *sinAlt, float *mag, const int num) c
onst
{ {
// going from observed position/magnitude to geometrical position an for (int i=0; i<num; ++i) mag[i] += airmass(sinAlt[i], true) * ext_c
d atmosphere-free mag. oeff;
for (int i=0; i<size; ++i) }
{ void Extinction::forward(const float *sinAlt, float *mag, const int num) co
float obs_alt_deg=(180.f/M_PI)*std::asin(altAzPos[i][2]); nst
if (obs_alt_deg > -2.f) {
{ for (int i=0; i<num; ++i) mag[i] += airmass(sinAlt[i], true) * ext_c
mag[i] -= airmass(altAzPos[i][2], true) * ext_coeff; oeff;
// refraction from Bennett, in Meeus, Astr.Alg. }
float r=press_temp_corr_Bennett / std::tan((obs_alt_ // from observed magnitude in apparent (observed) altitude to atmosphere-fr
deg+7.31/(obs_alt_deg+4.4f))*M_PI/180.f) + 0.0013515f; ee mag, still in apparent, refracted altitude.
obs_alt_deg -= r; void Extinction::backward(const Vec3d *altAzPos, float *mag, const int num)
altAzPos[i][2]=std::sin(obs_alt_deg*M_PI/180.f); 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 RefractionExtinction::airmass(float cosZ, bool apparent_z) float Extinction::airmass(const float cosZ, const bool apparent_z) const
{ {
if (cosZ<-0.035f) if (cosZ<-0.035f)
return 0.0f; // Safety: Do nothing for below -2 degrees. return Extinction::SUBHORIZONTAL_AIRMASS; // Safety: 0 or 40 for below -2 degrees.
float X; 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)); X=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; float nom=(1.002432f*cosZ+0.148386f)*cosZ+0.0096467f;
float denum=((cosZ+0.149864f)*cosZ+0.0102963f)*cosZ+0.000303 978f; float denum=((cosZ+0.149864f)*cosZ+0.0102963f)*cosZ+0.000303 978f;
X=nom/denum; X=nom/denum;
} }
return X; return X;
} }
void RefractionExtinction::setPressure(float p) // 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.
// -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
// 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.
// This must be -5 or higher.
const double Refraction::MIN_GEO_ALTITUDE_DEG=-3.54;
// 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;
// this must be positive. Transition zone goes that far below the values ju
st specified.
const double Refraction::TRANSITION_WIDTH_GEO_DEG=1.46;
const double Refraction::TRANSITION_WIDTH_APP_DEG=1.78217;
const double Refraction::MIN_GEO_ALTITUDE_RAD=Refraction::MIN_GEO_ALTITUDE_
DEG*M_PI/180.0;
const double Refraction::MIN_GEO_ALTITUDE_SIN=std::sin(Refraction::MIN_GEO_
ALTITUDE_RAD);
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),
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())
{ {
pressure=p;
updatePrecomputed(); updatePrecomputed();
} }
void RefractionExtinction::setTemperature(float t) void Refraction::setPreTransfoMat(const Mat4d& m)
{ {
temperature=t; preTransfoMat=m;
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]);
invertPreTransfoMatf.set(invertPreTransfoMat[0], invertPreTransfoMat
[1], invertPreTransfoMat[2], invertPreTransfoMat[3],
invertPreTransfoMat
[4], invertPreTransfoMat[5], invertPreTransfoMat[6], invertPreTransfoMat[7]
,
invertPreTransfoMat
[8], invertPreTransfoMat[9], invertPreTransfoMat[10], invertPreTransfoMat[1
1],
invertPreTransfoMat
[12], invertPreTransfoMat[13], invertPreTransfoMat[14], invertPreTransfoMat
[15]);
}
void Refraction::setPostTransfoMat(const Mat4d& m)
{
postTransfoMat=m;
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]);
invertPostTransfoMatf.set(invertPostTransfoMat[0], invertPostTransfo
Mat[1], invertPostTransfoMat[2], invertPostTransfoMat[3],
invertPostTransfoMa
t[4], invertPostTransfoMat[5], invertPostTransfoMat[6], invertPostTransfoMa
t[7],
invertPostTransfoMa
t[8], invertPostTransfoMat[9], invertPostTransfoMat[10], invertPostTransfoM
at[11],
invertPostTransfoMa
t[12], invertPostTransfoMat[13], invertPostTransfoMat[14], invertPostTransf
oMat[15]);
}
void Refraction::updatePrecomputed()
{
press_temp_corr_Bennett=pressure/1010.f * 283.f/(273.f+temperature)
/ 60.f;
press_temp_corr_Saemundson=1.02f*press_temp_corr_Bennett;
}
void Refraction::forward(Vec3d& altAzPos) const
{
altAzPos.transfo4d(preTransfoMat);
const double length = altAzPos.length();
double geom_alt_deg=180./M_PI*std::asin(altAzPos[2]/length);
if (geom_alt_deg > Refraction::MIN_GEO_ALTITUDE_DEG)
{
// 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;
geom_alt_deg += r;
if (geom_alt_deg > 90.) geom_alt_deg=90.; // SAFETY
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)
{
// Avoids the jump near -5 by interpolating linearly between
MIN_GEO_ALTITUDE_DEG and bottom of transition zone
float r_min=press_temp_corr_Saemundson / std::tan((Refractio
n::MIN_GEO_ALTITUDE_DEG+10.3f/(Refraction::MIN_GEO_ALTITUDE_DEG+5.11f))*M_P
I/180.f) + 0.0019279f;
geom_alt_deg += r_min*(geom_alt_deg-(Refraction::MIN_GEO_ALT
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.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);
// going from apparent (observed) position to geometrical position.
const double length = altAzPos.length();
double obs_alt_deg=180./M_PI*std::asin(altAzPos[2]/length);
if (obs_alt_deg > 0.22879)
{
// refraction directly from Bennett, in Meeus, Astr.Alg.
float r=press_temp_corr_Bennett / std::tan((obs_alt_deg+7.31
/(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.)*length;
}
else if (obs_alt_deg > Refraction::MIN_APP_ALTITUDE_DEG)
{
// 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;
obs_alt_deg -= press_temp_corr_Bennett*r;
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)
{
// use polynomial fit for topmost value, linear interpolatio
n inside transition zone
float r_min=(((((0.0444*Refraction::MIN_APP_ALTITUDE_DEG+.76
62)*Refraction::MIN_APP_ALTITUDE_DEG
+4.9746)*Refraction::MIN_APP_ALTITUDE_DEG+13
.599)*Refraction::MIN_APP_ALTITUDE_DEG
+8.052)*Refraction::MIN_APP_ALTITUDE_DEG-11.30
8)*Refraction::MIN_APP_ALTITUDE_DEG+34.341;
r_min*=press_temp_corr_Bennett;
obs_alt_deg -= r_min*(obs_alt_deg-(Refraction::MIN_APP_ALTIT
UDE_DEG-Refraction::TRANSITION_WIDTH_APP_DEG))/Refraction::TRANSITION_WIDTH
_APP_DEG;
altAzPos[2]=std::sin(obs_alt_deg*M_PI/180.)*length;
}
altAzPos.transfo4d(invertPreTransfoMat);
}
void Refraction::forward(Vec3f& altAzPos) const
{
altAzPos.transfo4d(preTransfoMatf);
const float length = altAzPos.length();
float geom_alt_deg=180.f/M_PI*std::asin(altAzPos[2]/length);
if (geom_alt_deg > Refraction::MIN_GEO_ALTITUDE_DEG_F)
{
// 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;
geom_alt_deg += r;
if (geom_alt_deg > 90.f) geom_alt_deg=90.f; // SAFETY
altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.f)*length;
}
else if(geom_alt_deg>Refraction::MIN_GEO_ALTITUDE_DEG_F-Refraction::
TRANSITION_WIDTH_GEO_DEG_F)
{
// Avoids the jump below -5 by interpolating linearly betwee
n MIN_GEO_ALTITUDE_DEG_F and bottom of transition zone
float r_m5=press_temp_corr_Saemundson / std::tan((Refraction
::MIN_GEO_ALTITUDE_DEG_F+10.3f/(Refraction::MIN_GEO_ALTITUDE_DEG_F+5.11f))*
M_PI/180.f) + 0.0019279f;
geom_alt_deg += r_m5*(geom_alt_deg-(Refraction::MIN_GEO_ALTI
TUDE_DEG_F-Refraction::TRANSITION_WIDTH_GEO_DEG_F))/Refraction::TRANSITION_
WIDTH_GEO_DEG_F;
altAzPos[2]=std::sin(geom_alt_deg*M_PI/180.)*length;
}
altAzPos.transfo4d(postTransfoMatf);
}
void Refraction::backward(Vec3f& altAzPos) const
{
altAzPos.transfo4d(invertPostTransfoMatf);
// going from observed position/magnitude to geometrical position an
d atmosphere-free mag.
const float length = altAzPos.length();
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);
}
void Refraction::setPressure(float p)
{
pressure=p;
updatePrecomputed(); updatePrecomputed();
} }
void RefractionExtinction::setExtinctionCoefficient(float k) void Refraction::setTemperature(float t)
{ {
ext_coeff=k; temperature=t;
updatePrecomputed();
} }
 End of changes. 17 change blocks. 
58 lines changed or deleted 318 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/