source: trunk/Mars/mcore/nova.h@ 20022

Last change on this file since 20022 was 19655, checked in by tbretz, 6 years ago
Added to more functions to be able to access the returncode for the object.
File size: 8.3 KB
Line 
1#ifndef MARS_Nova
2#define MARS_Nova
3
4#include <map>
5#include <cmath>
6#include <string>
7#include <limits>
8#include <algorithm>
9
10#include <libnova/solar.h>
11#include <libnova/lunar.h>
12#include <libnova/rise_set.h>
13#include <libnova/transform.h>
14#include <libnova/angular_separation.h>
15
16namespace Nova
17{
18 typedef ln_rst_time RstTime; // Note that times are in JD not Mjd
19
20 struct ZdAzPosn;
21 struct RaDecPosn;
22
23#ifndef __CINT__
24#define OBS(lat,lng) ={lat,lng}
25#else
26#define OBS(lat,lng)
27#endif
28
29#ifndef PRESET_OBSERVATORY
30#define PRESET_OBSERVATORY kORM
31#endif
32
33 const static double kSolarStandardHorizon = LN_SOLAR_STANDART_HORIZON;
34
35 const static ln_lnlat_posn kORM OBS(-(17.+53./60+26.525/3600), 28.+45./60+42.462/3600); // 2800m
36 const static ln_lnlat_posn kHAWC OBS( -97.3084, 18.9947 ); // 4100m;
37 const static ln_lnlat_posn kSPM OBS(-115.4637, 31.0439 ); // 2800m;
38 const static ln_lnlat_posn kRWTH OBS( 6.0657, 50.7801 );
39
40 struct LnLatPosn : public ln_lnlat_posn
41 {
42 typedef std::map<std::string, ln_lnlat_posn> LUT;
43
44 bool operator==(const LnLatPosn &eq) const
45 {
46 return lng==eq.lng && lat==eq.lat;
47 }
48
49 static const LUT &lut()
50 {
51#ifndef __CINT__
52 static const LUT lut =
53 {
54 { "ORM", kORM },
55 { "HAWC", kHAWC },
56 { "SPM", kSPM },
57 { "RWTH", kRWTH },
58 };
59 return lut;
60#else
61 return LUT();
62#endif
63 }
64
65 LnLatPosn() { }
66 LnLatPosn(const ln_lnlat_posn &pos) : ln_lnlat_posn(pos) { }
67 LnLatPosn(std::string obs)
68 {
69#ifndef __CINT__
70 if (obs.empty())
71 {
72 *this = PRESET_OBSERVATORY;
73 return;
74 }
75
76 for (auto it=obs.begin(); it!=obs.end(); it++)
77 *it = toupper(*it);
78
79 const auto it = lut().find(obs);
80 *this = it==lut().end() ? ln_lnlat_posn({std::numeric_limits<double>::quiet_NaN(),std::numeric_limits<double>::quiet_NaN()}) : it->second;
81#endif
82 }
83 bool isValid() const
84 {
85 return std::isfinite(lat) && std::isfinite(lng);
86 }
87
88 const std::string &name() const
89 {
90 for (auto it=lut().begin(); it!=lut().end(); it++)
91 if (LnLatPosn(it->second)==*this)
92 return it->first;
93
94 static const std::string dummy;
95 return dummy;
96 }
97
98 static const std::string preset()
99 {
100 return LnLatPosn(PRESET_OBSERVATORY).name();
101 }
102 };
103
104 // Warning: 0deg=South, 90deg=W
105 struct HrzPosn : public ln_hrz_posn
106 {
107 HrzPosn() { }
108 HrzPosn(const ZdAzPosn &);
109 };
110
111 struct ZdAzPosn
112 {
113 double zd; // [deg]
114 double az; // [deg]
115
116 ZdAzPosn(double z=0, double a=0) : zd(z), az(a) { }
117 ZdAzPosn(const HrzPosn &hrz) : zd(90-hrz.alt), az(hrz.az-180) { }
118
119 // This crahsed cint, but it could save up the dedicate structure HrzPosn :(
120 //operator HrzPosn() const { const HrzPosn p = { az-180, 90-zd}; return p; }
121 };
122
123 inline HrzPosn::HrzPosn(const ZdAzPosn &za) { alt = 90-za.zd; az = za.az-180; }
124
125
126 // Note that right ascension is stored in degree
127 struct EquPosn : public ln_equ_posn
128 {
129 EquPosn() { }
130 EquPosn(const RaDecPosn &);
131 };
132
133 struct RaDecPosn
134 {
135 double ra; // [h]
136 double dec; // [deg]
137
138 RaDecPosn(double r=0, double d=0) : ra(r), dec(d) { }
139 RaDecPosn(const EquPosn &equ) : ra(equ.ra/15), dec(equ.dec) { }
140
141 // This crahsed cint, but it could save up the dedicate structure HrzPosn :(
142 //operator HrzPosn() const { const HrzPosn p = { az-180, 90-zd}; return p; }
143 };
144
145 inline EquPosn::EquPosn(const RaDecPosn &rd) { ra = rd.ra*15; dec = rd.dec; }
146
147 inline HrzPosn GetHrzFromEqu(const EquPosn &equ, const LnLatPosn &obs, const double &jd)
148 {
149 HrzPosn hrz;
150 ln_get_hrz_from_equ(const_cast<EquPosn*>(&equ), const_cast<LnLatPosn*>(&obs), jd, &hrz);
151 return hrz;
152 }
153 inline HrzPosn GetHrzFromEqu(const EquPosn &equ, const double &jd)
154 {
155 return GetHrzFromEqu(equ, PRESET_OBSERVATORY, jd);
156 }
157
158 inline EquPosn GetEquFromHrz(const HrzPosn &hrz, const LnLatPosn &obs, const double &jd)
159 {
160 EquPosn equ;
161 ln_get_equ_from_hrz(const_cast<HrzPosn*>(&hrz), const_cast<LnLatPosn*>(&obs), jd, &equ);
162 return equ;
163 }
164 inline EquPosn GetEquFromHrz(const HrzPosn &hrz, const double &jd)
165 {
166 return GetEquFromHrz(hrz, PRESET_OBSERVATORY, jd);
167 }
168
169 inline RstTime GetObjectRst(const EquPosn &equ, const LnLatPosn &obs, const double &jd, const double &hrz=0)
170 {
171 RstTime rst;
172 // 0 for success, 1 for circumpolar (above the horizon), -1 for circumpolar (bellow the horizon)
173 ln_get_object_next_rst_horizon(jd, const_cast<LnLatPosn*>(&obs), const_cast<EquPosn*>(&equ), hrz, &rst);
174 return rst;
175 }
176 inline RstTime GetObjectRst(const EquPosn &equ, const double &jd, const double &hrz=0)
177 {
178 return GetObjectRst(equ, PRESET_OBSERVATORY, jd, hrz);
179 }
180
181 inline int GetObjectRst(RstTime &rst, const EquPosn &equ, const LnLatPosn &obs, const double &jd, const double &hrz=0)
182 {
183 // 0 for success, 1 for circumpolar (above the horizon), -1 for circumpolar (bellow the horizon)
184 return ln_get_object_next_rst_horizon(jd, const_cast<LnLatPosn*>(&obs), const_cast<EquPosn*>(&equ), hrz, &rst);
185 }
186 inline int GetObjectRst(RstTime &rst, const EquPosn &equ, const double &jd, const double &hrz=0)
187 {
188 // 0 for success, 1 for circumpolar (above the horizon), -1 for circumpolar (bellow the horizon)
189 return GetObjectRst(rst, equ, PRESET_OBSERVATORY, jd, hrz);
190 }
191
192 inline RstTime GetSolarRst(const double &jd, const LnLatPosn &obs, const double &hrz=kSolarStandardHorizon)
193 {
194 RstTime rst;
195 ln_get_solar_rst_horizon(jd, const_cast<LnLatPosn*>(&obs), hrz, &rst);
196 return rst;
197 }
198
199 inline RstTime GetSolarRst(const double &jd, const double &hrz=kSolarStandardHorizon)
200 {
201 return GetSolarRst(jd, PRESET_OBSERVATORY, hrz);
202 }
203
204 inline RstTime GetLunarRst(const double &jd, const LnLatPosn &obs=PRESET_OBSERVATORY)
205 {
206 RstTime rst;
207 ln_get_lunar_rst(jd, const_cast<LnLatPosn*>(&obs), &rst);
208 return rst;
209 }
210 inline EquPosn GetSolarEquCoords(const double &jd)
211 {
212 EquPosn equ;
213 ln_get_solar_equ_coords(jd, &equ);
214 return equ;
215 }
216
217 inline double GetLunarDisk(const double &jd)
218 {
219 return ln_get_lunar_disk(jd);
220 }
221
222 inline double GetLunarSdiam(const double &jd)
223 {
224 return ln_get_lunar_sdiam(jd);
225 }
226
227 inline double GetLunarPhase(const double &jd)
228 {
229 return ln_get_lunar_phase(jd);
230 }
231
232 inline EquPosn GetLunarEquCoords(const double &jd, double precision=0)
233 {
234 EquPosn equ;
235 ln_get_lunar_equ_coords_prec(jd, &equ, precision);
236 return equ;
237 }
238
239 inline double GetLunarEarthDist(const double &jd)
240 {
241 return ln_get_lunar_earth_dist(jd);
242 }
243
244 inline double GetAngularSeparation(const EquPosn &p1, const EquPosn &p2)
245 {
246 return ln_get_angular_separation(const_cast<EquPosn*>(&p1), const_cast<EquPosn*>(&p2));
247 }
248
249 inline double GetAngularSeparation(const HrzPosn &h1, const HrzPosn &h2)
250 {
251 EquPosn p1; p1.ra=h1.az; p1.dec=h1.alt;
252 EquPosn p2; p2.ra=h2.az; p2.dec=h2.alt;
253 return ln_get_angular_separation(&p1, &p2);
254 }
255
256 struct SolarObjects
257 {
258 double fJD;
259
260 EquPosn fSunEqu;
261 HrzPosn fSunHrz;
262
263 EquPosn fMoonEqu;
264 HrzPosn fMoonHrz;
265 double fMoonDisk;
266
267 double fEarthDist;
268
269 SolarObjects() { }
270
271 SolarObjects(const double &jd, const LnLatPosn &obs=Nova::PRESET_OBSERVATORY)
272 {
273 fJD = jd;
274
275 // Sun properties
276 fSunEqu = GetSolarEquCoords(jd);
277 fSunHrz = GetHrzFromEqu(fSunEqu, obs, jd);
278
279 // Moon properties
280 fMoonEqu = GetLunarEquCoords(jd, 0.01);
281 fMoonHrz = GetHrzFromEqu(fMoonEqu, obs, jd);
282
283 fMoonDisk = GetLunarDisk(jd);
284 fEarthDist = GetLunarEarthDist(jd);
285 }
286 };
287}
288
289#endif
Note: See TracBrowser for help on using the repository browser.