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

Last change on this file since 19394 was 19394, checked in by tbretz, 6 years ago
Minor fix.
File size: 7.1 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 GetSolarRst(const double &jd, const LnLatPosn &obs, const double &hrz=kSolarStandardHorizon)
170 {
171 RstTime rst;
172 ln_get_solar_rst_horizon(jd, const_cast<LnLatPosn*>(&obs), hrz, &rst);
173 return rst;
174 }
175 inline RstTime GetSolarRst(const double &jd, const double &hrz=kSolarStandardHorizon)
176 {
177 return GetSolarRst(jd, PRESET_OBSERVATORY, hrz);
178 }
179
180 inline RstTime GetLunarRst(const double &jd, const LnLatPosn &obs=PRESET_OBSERVATORY)
181 {
182 RstTime rst;
183 ln_get_lunar_rst(jd, const_cast<LnLatPosn*>(&obs), &rst);
184 return rst;
185 }
186 inline EquPosn GetSolarEquCoords(const double &jd)
187 {
188 EquPosn equ;
189 ln_get_solar_equ_coords(jd, &equ);
190 return equ;
191 }
192
193 inline double GetLunarDisk(const double &jd)
194 {
195 return ln_get_lunar_disk(jd);
196 }
197
198 inline double GetLunarSdiam(const double &jd)
199 {
200 return ln_get_lunar_sdiam(jd);
201 }
202
203 inline double GetLunarPhase(const double &jd)
204 {
205 return ln_get_lunar_phase(jd);
206 }
207
208 inline EquPosn GetLunarEquCoords(const double &jd, double precision=0)
209 {
210 EquPosn equ;
211 ln_get_lunar_equ_coords_prec(jd, &equ, precision);
212 return equ;
213 }
214
215 inline double GetLunarEarthDist(const double &jd)
216 {
217 return ln_get_lunar_earth_dist(jd);
218 }
219
220 inline double GetAngularSeparation(const EquPosn &p1, const EquPosn &p2)
221 {
222 return ln_get_angular_separation(const_cast<EquPosn*>(&p1), const_cast<EquPosn*>(&p2));
223 }
224
225 inline double GetAngularSeparation(const HrzPosn &h1, const HrzPosn &h2)
226 {
227 EquPosn p1; p1.ra=h1.az; p1.dec=h1.alt;
228 EquPosn p2; p2.ra=h2.az; p2.dec=h2.alt;
229 return ln_get_angular_separation(&p1, &p2);
230 }
231
232 struct SolarObjects
233 {
234 double fJD;
235
236 EquPosn fSunEqu;
237 HrzPosn fSunHrz;
238
239 EquPosn fMoonEqu;
240 HrzPosn fMoonHrz;
241 double fMoonDisk;
242
243 double fEarthDist;
244
245 SolarObjects(const double &jd, const LnLatPosn &obs=Nova::PRESET_OBSERVATORY)
246 {
247 fJD = jd;
248
249 // Sun properties
250 fSunEqu = GetSolarEquCoords(jd);
251 fSunHrz = GetHrzFromEqu(fSunEqu, obs, jd);
252
253 // Moon properties
254 fMoonEqu = GetLunarEquCoords(jd, 0.01);
255 fMoonHrz = GetHrzFromEqu(fMoonEqu, obs, jd);
256
257 fMoonDisk = GetLunarDisk(jd);
258 fEarthDist = GetLunarEarthDist(jd);
259 }
260 };
261}
262
263#endif
Note: See TracBrowser for help on using the repository browser.