1 /*
2 * Stellarium
3 * Copyright (C) 2010 Fabien Chereau
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License
7 * as published by the Free Software Foundation; either version 2
8 * of the License, or (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
18 *
19 * Refraction and extinction computations.
20 * Principal implementation: 2010-03-23 GZ=Georg Zotti, Georg.Zotti@univie.ac.at
21 */
22
23 #include "StelApp.hpp"
24 #include "StelUtils.hpp"
25 #include "RefractionExtinction.hpp"
26 #include "StelUtils.hpp"
27
Extinction()28 Extinction::Extinction() : ext_coeff(50), undergroundExtinctionMode(UndergroundExtinctionMirror)
29 {
30 }
31
32 // airmass computation for cosine of zenith angle z
airmass(float cosZ,const bool apparent_z) const33 float Extinction::airmass(float cosZ, const bool apparent_z) const
34 {
35 if (cosZ<-0.035f) // about -2 degrees. Here, RozenbergZ>574 and climbs fast!
36 {
37 switch (undergroundExtinctionMode)
38 {
39 case UndergroundExtinctionZero:
40 return 0.f;
41 case UndergroundExtinctionMax:
42 return 42.f;
43 case UndergroundExtinctionMirror:
44 cosZ = std::min(1.f, -0.035f - (cosZ+0.035f));
45 }
46 }
47
48 if (apparent_z)
49 {
50 // Rozenberg 1966, reported by Schaefer (1993-2000).
51 return 1.0f/(cosZ+0.025f*std::exp(-11.f*cosZ));
52 }
53 else
54 {
55 //Young 1994
56 const float nom=(1.002432f*cosZ+0.148386f)*cosZ+0.0096467f;
57 const float denum=((cosZ+0.149864f)*cosZ+0.0102963f)*cosZ+0.000303978f;
58 return nom/denum;
59 }
60 }
61
62 /* ***************************************************************************************************** */
63
64 // The following 4 are to be configured, the rest is derived.
65 // Recommendations: -4.9/-4.3/0.1/0.1: sharp but continuous transition, no effects below -5.
66 // -4.3/-4.3/0.7/0.7: sharp but continuous transition, no effects below -5. Maybe better for picking?
67 // -3/-3/2/2: "strange" zone 2 degrees wide. Both formulae are close near -3. Actually, refraction formulae dubious below 0
68 // 0/0/1/1: "strange zone 1 degree wide just below horizon, no effect below -1. Actually, refraction formulae dubious below 0! But it looks stupid, see the sun!
69 // Optimum:-3.54/-3.21783/1.46/1.78217. Here forward/backward are almost perfect inverses (-->good picking!), and nothing happens below -5 degrees.
70 // This must be -5 or higher.
71 static const float MIN_GEO_ALTITUDE_DEG=-3.54f;
72 // this must be -4.3 or higher. -3.21783 is an optimal value to fit against forward refraction!
73 static const float MIN_APP_ALTITUDE_DEG=-3.21783f;
74 // this must be positive. Transition zone goes that far below the values just specified.
75 static const float TRANSITION_WIDTH_GEO_DEG=1.46f;
76 static const float TRANSITION_WIDTH_APP_DEG=1.78217f;
77
Refraction()78 Refraction::Refraction() : pressure(1013.f), temperature(10.f),
79 preTransfoMat(Mat4d::identity()), invertPreTransfoMat(Mat4d::identity()), preTransfoMatf(Mat4f::identity()), invertPreTransfoMatf(Mat4f::identity()),
80 postTransfoMat(Mat4d::identity()), invertPostTransfoMat(Mat4d::identity()), postTransfoMatf(Mat4f::identity()), invertPostTransfoMatf(Mat4f::identity())
81 {
82 updatePrecomputed();
83 }
84
setPreTransfoMat(const Mat4d & m)85 void Refraction::setPreTransfoMat(const Mat4d& m)
86 {
87 preTransfoMat=m;
88 invertPreTransfoMat=m.inverse();
89 preTransfoMatf.set(static_cast<float>(m[0]), static_cast<float>(m[1]), static_cast<float>(m[2]), static_cast<float>(m[3]),
90 static_cast<float>(m[4]), static_cast<float>(m[5]), static_cast<float>(m[6]), static_cast<float>(m[7]),
91 static_cast<float>(m[8]), static_cast<float>(m[9]), static_cast<float>(m[10]), static_cast<float>(m[11]),
92 static_cast<float>(m[12]), static_cast<float>(m[13]), static_cast<float>(m[14]), static_cast<float>(m[15]));
93 invertPreTransfoMatf.set(static_cast<float>(invertPreTransfoMat[0]), static_cast<float>(invertPreTransfoMat[1]), static_cast<float>(invertPreTransfoMat[2]), static_cast<float>(invertPreTransfoMat[3]),
94 static_cast<float>(invertPreTransfoMat[4]), static_cast<float>(invertPreTransfoMat[5]), static_cast<float>(invertPreTransfoMat[6]), static_cast<float>(invertPreTransfoMat[7]),
95 static_cast<float>(invertPreTransfoMat[8]), static_cast<float>(invertPreTransfoMat[9]), static_cast<float>(invertPreTransfoMat[10]), static_cast<float>(invertPreTransfoMat[11]),
96 static_cast<float>(invertPreTransfoMat[12]), static_cast<float>(invertPreTransfoMat[13]), static_cast<float>(invertPreTransfoMat[14]), static_cast<float>(invertPreTransfoMat[15]));
97 }
98
setPostTransfoMat(const Mat4d & m)99 void Refraction::setPostTransfoMat(const Mat4d& m)
100 {
101 postTransfoMat=m;
102 invertPostTransfoMat=m.inverse();
103 postTransfoMatf.set(static_cast<float>(m[0]), static_cast<float>(m[1]), static_cast<float>(m[2]), static_cast<float>(m[3]),
104 static_cast<float>(m[4]), static_cast<float>(m[5]), static_cast<float>(m[6]), static_cast<float>(m[7]),
105 static_cast<float>(m[8]), static_cast<float>(m[9]), static_cast<float>(m[10]), static_cast<float>(m[11]),
106 static_cast<float>(m[12]), static_cast<float>(m[13]), static_cast<float>(m[14]), static_cast<float>(m[15]));
107 invertPostTransfoMatf.set(static_cast<float>(invertPostTransfoMat[0]), static_cast<float>(invertPostTransfoMat[1]), static_cast<float>(invertPostTransfoMat[2]), static_cast<float>(invertPostTransfoMat[3]),
108 static_cast<float>(invertPostTransfoMat[4]), static_cast<float>(invertPostTransfoMat[5]), static_cast<float>(invertPostTransfoMat[6]), static_cast<float>(invertPostTransfoMat[7]),
109 static_cast<float>(invertPostTransfoMat[8]), static_cast<float>(invertPostTransfoMat[9]), static_cast<float>(invertPostTransfoMat[10]), static_cast<float>(invertPostTransfoMat[11]),
110 static_cast<float>(invertPostTransfoMat[12]), static_cast<float>(invertPostTransfoMat[13]), static_cast<float>(invertPostTransfoMat[14]), static_cast<float>(invertPostTransfoMat[15]));
111 }
112
updatePrecomputed()113 void Refraction::updatePrecomputed()
114 {
115 press_temp_corr=pressure/1010.f * 283.f/(273.f+temperature) / 60.f;
116 }
117
innerRefractionForward(Vec3d & altAzPos) const118 void Refraction::innerRefractionForward(Vec3d& altAzPos) const
119 {
120 const double length = altAzPos.length();
121 if (length==0.0)
122 {
123 // Under some circumstances there are zero coordinates. Just leave them alone.
124 //qDebug() << "Refraction::innerRefractionForward(): Zero vector detected - Continue with zero vector.";
125 return;
126 }
127
128 Q_ASSERT(length>0.0);
129 const double sinGeo = altAzPos[2]/length;
130 Q_ASSERT(fabs(sinGeo)<=1.0);
131 double geom_alt_rad = std::asin(sinGeo);
132 float geom_alt_deg = M_180_PIf*static_cast<float>(geom_alt_rad);
133 if (geom_alt_deg > MIN_GEO_ALTITUDE_DEG)
134 {
135 // refraction from Saemundsson, S&T1986 p70 / in Meeus, Astr.Alg.
136 float r=press_temp_corr * ( 1.02f / std::tan((geom_alt_deg+10.3f/(geom_alt_deg+5.11f))*M_PI_180f) + 0.0019279f);
137 geom_alt_deg += r;
138 if (geom_alt_deg > 90.f)
139 geom_alt_deg=90.f;
140 }
141 else if(geom_alt_deg>MIN_GEO_ALTITUDE_DEG-TRANSITION_WIDTH_GEO_DEG)
142 {
143 // Avoids the jump below -5 by interpolating linearly between MIN_GEO_ALTITUDE_DEG and bottom of transition zone
144 float r_m5=press_temp_corr * ( 1.02f / std::tan((MIN_GEO_ALTITUDE_DEG+10.3f/(MIN_GEO_ALTITUDE_DEG+5.11f))*M_PI_180f) + 0.0019279f);
145 geom_alt_deg += r_m5*(geom_alt_deg-(MIN_GEO_ALTITUDE_DEG-TRANSITION_WIDTH_GEO_DEG))/TRANSITION_WIDTH_GEO_DEG;
146 }
147 else return;
148 // At this point we have corrected geometric altitude. Note that if we just change altAzPos[2], we would change vector length, so this would change our angles.
149 // We have to shorten X,Y components of the vector as well by the change in cosines of altitude, or (sqrt(1-sin(alt))
150
151 const double refr_alt_rad=static_cast<double>(geom_alt_deg)*M_PI_180;
152 const double sinRef=std::sin(refr_alt_rad);
153
154 const double shortenxy=((fabs(sinGeo)>=1.0) ? 1.0 :
155 std::sqrt((1.-sinRef*sinRef)/(1.-sinGeo*sinGeo))); // we need double's mantissa length here, sorry!
156
157 altAzPos[0]*=shortenxy;
158 altAzPos[1]*=shortenxy;
159 altAzPos[2]=sinRef*length;
160 }
161
162 // going from observed position to geometrical position.
innerRefractionBackward(Vec3d & altAzPos) const163 void Refraction::innerRefractionBackward(Vec3d& altAzPos) const
164 {
165 const double length = altAzPos.length();
166 if (length==0.0)
167 {
168 // Under some circumstances there are zero coordinates. Just leave them alone.
169 //qDebug() << "Refraction::innerRefractionBackward(): Zero vector detected - Continue with zero vector.";
170 return;
171 }
172 Q_ASSERT(length>0.0);
173 const double sinObs = altAzPos[2]/length;
174 Q_ASSERT(fabs(sinObs)<=1.0);
175 float obs_alt_deg=static_cast<float>(M_180_PI*std::asin(sinObs));
176 if (obs_alt_deg > 0.22879f)
177 {
178 // refraction from Bennett, in Meeus, Astr.Alg.
179 float r=press_temp_corr * (1.f / std::tan((obs_alt_deg+7.31f/(obs_alt_deg+4.4f))*M_PI_180f) + 0.0013515f);
180 obs_alt_deg -= r;
181 }
182 else if (obs_alt_deg > MIN_APP_ALTITUDE_DEG)
183 {
184 // backward refraction from polynomial fit against Saemundson[-5...-0.3]
185 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;
186 obs_alt_deg -= press_temp_corr*r;
187 }
188 else if (obs_alt_deg > MIN_APP_ALTITUDE_DEG-TRANSITION_WIDTH_APP_DEG)
189 {
190 // Compute top value from polynome, apply linear interpolation
191 static const float r_min=(((((0.0444f*MIN_APP_ALTITUDE_DEG+.7662f)*MIN_APP_ALTITUDE_DEG
192 +4.9746f)*MIN_APP_ALTITUDE_DEG+13.599f)*MIN_APP_ALTITUDE_DEG
193 +8.052f)*MIN_APP_ALTITUDE_DEG-11.308f)*MIN_APP_ALTITUDE_DEG+34.341f;
194
195 obs_alt_deg -= r_min*press_temp_corr*(obs_alt_deg-(MIN_APP_ALTITUDE_DEG-TRANSITION_WIDTH_APP_DEG))/TRANSITION_WIDTH_APP_DEG;
196 }
197 else return;
198 // At this point we have corrected observed altitude. Note that if we just change altAzPos[2], we would change vector length, so this would change our angles.
199 // 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))
200
201 const double geo_alt_rad=static_cast<double>(obs_alt_deg)*M_PI_180;
202 const double sinGeo=std::sin(geo_alt_rad);
203 const double longerxy=((fabs(sinObs)>=1.0) ? 1.0 :
204 std::sqrt((1.-sinGeo*sinGeo)/(1.-sinObs*sinObs)));
205 altAzPos[0]*=longerxy;
206 altAzPos[1]*=longerxy;
207 altAzPos[2]=sinGeo*length;
208 }
209
forward(Vec3d & altAzPos) const210 void Refraction::forward(Vec3d& altAzPos) const
211 {
212 altAzPos.transfo4d(preTransfoMat);
213 innerRefractionForward(altAzPos);
214 altAzPos.transfo4d(postTransfoMat);
215 }
216
217 //Bennett's formula is not a strict inverse of Saemundsson's. There is a notable discrepancy (alt!=backward(forward(alt))) for
218 // geometric altitudes <-.3deg. This is not a problem in real life, but if a user switches off landscape, click-identify may pose a problem.
219 // Below this altitude, we therefore use a polynomial fit that should represent a close inverse of Saemundsson.
backward(Vec3d & altAzPos) const220 void Refraction::backward(Vec3d& altAzPos) const
221 {
222 altAzPos.transfo4d(invertPostTransfoMat);
223 innerRefractionBackward(altAzPos);
224 altAzPos.transfo4d(invertPreTransfoMat);
225 }
226
forward(Vec3f & altAzPos) const227 void Refraction::forward(Vec3f& altAzPos) const
228 {
229 Vec3d vf=altAzPos.toVec3d();
230 vf.transfo4d(preTransfoMat);
231 innerRefractionForward(vf);
232 vf.transfo4d(postTransfoMat);
233 altAzPos.set(static_cast<float>(vf[0]), static_cast<float>(vf[1]), static_cast<float>(vf[2]));
234 }
235
backward(Vec3f & altAzPos) const236 void Refraction::backward(Vec3f& altAzPos) const
237 {
238 altAzPos.transfo4d(invertPostTransfoMatf);
239 Vec3d vf=altAzPos.toVec3d();
240 innerRefractionBackward(vf);
241 altAzPos.set(static_cast<float>(vf[0]), static_cast<float>(vf[1]), static_cast<float>(vf[2]));
242 altAzPos.transfo4d(invertPreTransfoMatf);
243 }
244
setPressure(float p)245 void Refraction::setPressure(float p)
246 {
247 pressure=p;
248 updatePrecomputed();
249 }
250
setTemperature(float t)251 void Refraction::setTemperature(float t)
252 {
253 temperature=t;
254 updatePrecomputed();
255 }
256