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

Last change on this file since 19388 was 19388, checked in by tbretz, 13 months ago
Declare all functions inline to avoid multiple definitions when including the header twice. Use the precompiler directive PRESET_OBSERVATORY as default observatory instead of kORM (keep kORM as default if PRESET_OBSERVATORY was not set. Implement a function which allows to return the name corresponding to some observatory location. Implemented a function which return the corresponding name of the PRESET_OBSERVATORY. If the string in the constructor is empty, return the PRESET_OBSERVATORY. This allows to overwrite the default during compile time.
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                return PRESET_OBSERVATORY;
72
73            for (auto it=obs.begin(); it!=obs.end(); it++)
74                *it = toupper(*it);
75
76            const auto it = lut().find(obs);
77            *this = it==lut().end() ? ln_lnlat_posn({std::numeric_limits<double>::quiet_NaN(),std::numeric_limits<double>::quiet_NaN()}) : it->second;
78#endif
79        }
80        bool isValid() const
81        {
82            return std::isfinite(lat) && std::isfinite(lng);
83        }
84
85        const std::string &name() const
86        {
87            for (auto it=lut().begin(); it!=lut().end(); it++)
88                if (LnLatPosn(it->second)==*this)
89                    return it->first;
90
91            static const std::string dummy;
92            return dummy;
93        }
94
95        static const std::string preset()
96        {
97            return LnLatPosn(PRESET_OBSERVATORY).name();
98        }
99    };
100
101    // Warning: 0deg=South, 90deg=W
102    struct HrzPosn : public ln_hrz_posn
103    {
104        HrzPosn() { }
105        HrzPosn(const ZdAzPosn &);
106    };
107
108    struct ZdAzPosn
109    {
110        double zd; // [deg]
111        double az; // [deg]
112
113        ZdAzPosn(double z=0, double a=0) : zd(z), az(a) { }
114        ZdAzPosn(const HrzPosn &hrz) : zd(90-hrz.alt), az(hrz.az-180) { }
115
116        // This crahsed cint, but it could save up the dedicate structure HrzPosn :(
117        //operator HrzPosn() const { const HrzPosn p = { az-180, 90-zd}; return p; }
118    };
119
120    inline HrzPosn::HrzPosn(const ZdAzPosn &za) { alt = 90-za.zd; az = za.az-180; }
121
122
123    // Note that right ascension is stored in degree
124    struct EquPosn : public ln_equ_posn
125    {
126        EquPosn() { }
127        EquPosn(const RaDecPosn &);
128    };
129
130    struct RaDecPosn
131    {
132        double ra;  // [h]
133        double dec; // [deg]
134
135        RaDecPosn(double r=0, double d=0) : ra(r), dec(d) { }
136        RaDecPosn(const EquPosn &equ) : ra(equ.ra/15), dec(equ.dec) { }
137
138        // This crahsed cint, but it could save up the dedicate structure HrzPosn :(
139        //operator HrzPosn() const { const HrzPosn p = { az-180, 90-zd}; return p; }
140    };
141
142    inline EquPosn::EquPosn(const RaDecPosn &rd) { ra = rd.ra*15; dec = rd.dec; }
143
144    inline HrzPosn GetHrzFromEqu(const EquPosn &equ, const LnLatPosn &obs, const double &jd)
145    {
146        HrzPosn hrz;
147        ln_get_hrz_from_equ(const_cast<EquPosn*>(&equ), const_cast<LnLatPosn*>(&obs), jd, &hrz);
148        return hrz;
149    }
150    inline HrzPosn GetHrzFromEqu(const EquPosn &equ, const double &jd)
151    {
152        return GetHrzFromEqu(equ, PRESET_OBSERVATORY, jd);
153    }
154
155    inline EquPosn GetEquFromHrz(const HrzPosn &hrz, const LnLatPosn &obs, const double &jd)
156    {
157        EquPosn equ;
158        ln_get_equ_from_hrz(const_cast<HrzPosn*>(&hrz), const_cast<LnLatPosn*>(&obs), jd, &equ);
159        return equ;
160    }
161    inline EquPosn GetEquFromHrz(const HrzPosn &hrz, const double &jd)
162    {
163        return GetEquFromHrz(hrz, PRESET_OBSERVATORY, jd);
164    }
165
166    inline RstTime GetSolarRst(const double &jd, const LnLatPosn &obs, const double &hrz=kSolarStandardHorizon)
167    {
168        RstTime rst;
169        ln_get_solar_rst_horizon(jd, const_cast<LnLatPosn*>(&obs), hrz, &rst);
170        return rst;
171    }
172    inline RstTime GetSolarRst(const double &jd, const double &hrz=kSolarStandardHorizon)
173    {
174        return GetSolarRst(jd, PRESET_OBSERVATORY, hrz);
175    }
176
177    inline RstTime GetLunarRst(const double &jd, const LnLatPosn &obs=PRESET_OBSERVATORY)
178    {
179        RstTime rst;
180        ln_get_lunar_rst(jd, const_cast<LnLatPosn*>(&obs), &rst);
181        return rst;
182    }
183    inline EquPosn GetSolarEquCoords(const double &jd)
184    {
185        EquPosn equ;
186        ln_get_solar_equ_coords(jd, &equ);
187        return equ;
188    }
189
190    inline double GetLunarDisk(const double &jd)
191    {
192        return ln_get_lunar_disk(jd);
193    }
194
195    inline double GetLunarSdiam(const double &jd)
196    {
197        return ln_get_lunar_sdiam(jd);
198    }
199
200    inline double GetLunarPhase(const double &jd)
201    {
202        return ln_get_lunar_phase(jd);
203    }
204
205    inline EquPosn GetLunarEquCoords(const double &jd, double precision=0)
206    {
207        EquPosn equ;
208        ln_get_lunar_equ_coords_prec(jd, &equ, precision);
209        return equ;
210    }
211
212    inline double GetLunarEarthDist(const double &jd)
213    {
214        return ln_get_lunar_earth_dist(jd);
215    }
216
217    inline double GetAngularSeparation(const EquPosn &p1, const EquPosn &p2)
218    {
219        return ln_get_angular_separation(const_cast<EquPosn*>(&p1), const_cast<EquPosn*>(&p2));
220    }
221
222    inline double GetAngularSeparation(const HrzPosn &h1, const HrzPosn &h2)
223    {
224        EquPosn p1; p1.ra=h1.az; p1.dec=h1.alt;
225        EquPosn p2; p2.ra=h2.az; p2.dec=h2.alt;
226        return ln_get_angular_separation(&p1, &p2);
227    }
228
229    struct SolarObjects
230    {
231        double  fJD;
232
233        EquPosn fSunEqu;
234        HrzPosn fSunHrz;
235
236        EquPosn fMoonEqu;
237        HrzPosn fMoonHrz;
238        double  fMoonDisk;
239
240        double fEarthDist;
241
242        SolarObjects(const double &jd, const LnLatPosn &obs=Nova::PRESET_OBSERVATORY)
243        {
244            fJD = jd;
245
246            // Sun properties
247            fSunEqu    = GetSolarEquCoords(jd);
248            fSunHrz    = GetHrzFromEqu(fSunEqu, obs, jd);
249
250            // Moon properties
251            fMoonEqu   = GetLunarEquCoords(jd, 0.01);
252            fMoonHrz   = GetHrzFromEqu(fMoonEqu, obs, jd);
253
254            fMoonDisk  = GetLunarDisk(jd);
255            fEarthDist = GetLunarEarthDist(jd);
256        }
257    };
258}
259
260#endif
Note: See TracBrowser for help on using the repository browser.