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

Last change on this file since 19388 was 19388, checked in by tbretz, 8 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.