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

Last change on this file since 19293 was 19293, checked in by tbretz, 14 months ago
kSolarStandardHorizon for convenience (no need to include libnova.h and no need to copy all these variables. A refernce is enough.
File size: 6.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    const static double kSolarStandardHorizon = LN_SOLAR_STANDART_HORIZON;
30
31    const static ln_lnlat_posn kORM  OBS(-(17.+53./60+26.525/3600), 28.+45./60+42.462/3600); // 2800m
32    const static ln_lnlat_posn kHAWC OBS( -97.3084,                 18.9947               ); // 4100m;
33    const static ln_lnlat_posn kSPM  OBS(-115.4637,                 31.0439               ); // 2800m;
34    const static ln_lnlat_posn kRWTH OBS(   6.0657,                 50.7801               );
35
36    struct LnLatPosn : public ln_lnlat_posn
37    {
38        LnLatPosn() { }
39        LnLatPosn(const ln_lnlat_posn &pos) : ln_lnlat_posn(pos) { }
40        LnLatPosn(std::string obs)
41        {
42#ifndef __CINT__
43            static const std::map<std::string, ln_lnlat_posn> lut =
44            {
45                { "ORM",  kORM  },
46                { "HAWC", kHAWC },
47                { "SPM",  kSPM  },
48                { "RWTH", kRWTH },
49            };
50
51            for (auto it=obs.begin(); it!=obs.end(); it++)
52                *it = toupper(*it);
53
54            const auto it = lut.find(obs);
55            *this = it==lut.end() ? ln_lnlat_posn({std::numeric_limits<double>::quiet_NaN(),std::numeric_limits<double>::quiet_NaN()}) : it->second;
56#endif
57        }
58        bool isValid() const
59        {
60            return std::isfinite(lat) && std::isfinite(lng);
61        }
62    };
63
64    // Warning: 0deg=South, 90deg=W
65    struct HrzPosn : public ln_hrz_posn
66    {
67        HrzPosn() { }
68        HrzPosn(const ZdAzPosn &);
69    };
70
71    struct ZdAzPosn
72    {
73        double zd; // [deg]
74        double az; // [deg]
75
76        ZdAzPosn(double z=0, double a=0) : zd(z), az(a) { }
77        ZdAzPosn(const HrzPosn &hrz) : zd(90-hrz.alt), az(hrz.az-180) { }
78
79        // This crahsed cint, but it could save up the dedicate structure HrzPosn :(
80        //operator HrzPosn() const { const HrzPosn p = { az-180, 90-zd}; return p; }
81    };
82
83    HrzPosn::HrzPosn(const ZdAzPosn &za) { alt = 90-za.zd; az = za.az-180; }
84
85
86    // Note that right ascension is stored in degree
87    struct EquPosn : public ln_equ_posn
88    {
89        EquPosn() { }
90        EquPosn(const RaDecPosn &);
91    };
92
93    struct RaDecPosn
94    {
95        double ra;  // [h]
96        double dec; // [deg]
97
98        RaDecPosn(double r=0, double d=0) : ra(r), dec(d) { }
99        RaDecPosn(const EquPosn &equ) : ra(equ.ra/15), dec(equ.dec) { }
100
101        // This crahsed cint, but it could save up the dedicate structure HrzPosn :(
102        //operator HrzPosn() const { const HrzPosn p = { az-180, 90-zd}; return p; }
103    };
104
105    EquPosn::EquPosn(const RaDecPosn &rd) { ra = rd.ra*15; dec = rd.dec; }
106
107    HrzPosn GetHrzFromEqu(const EquPosn &equ, const LnLatPosn &obs, const double &jd)
108    {
109        HrzPosn hrz;
110        ln_get_hrz_from_equ(const_cast<EquPosn*>(&equ), const_cast<LnLatPosn*>(&obs), jd, &hrz);
111        return hrz;
112    }
113    HrzPosn GetHrzFromEqu(const EquPosn &equ, const double &jd)
114    {
115        return GetHrzFromEqu(equ, kORM, jd);
116    }
117
118    EquPosn GetEquFromHrz(const HrzPosn &hrz, const LnLatPosn &obs, const double &jd)
119    {
120        EquPosn equ;
121        ln_get_equ_from_hrz(const_cast<HrzPosn*>(&hrz), const_cast<LnLatPosn*>(&obs), jd, &equ);
122        return equ;
123    }
124    EquPosn GetEquFromHrz(const HrzPosn &hrz, const double &jd)
125    {
126        return GetEquFromHrz(hrz, kORM, jd);
127    }
128
129    RstTime GetSolarRst(const double &jd, const LnLatPosn &obs, const double &hrz=kSolarStandardHorizon)
130    {
131        RstTime rst;
132        ln_get_solar_rst_horizon(jd, const_cast<LnLatPosn*>(&obs), hrz, &rst);
133        return rst;
134    }
135    RstTime GetSolarRst(const double &jd, const double &hrz=kSolarStandardHorizon)
136    {
137        return GetSolarRst(jd, kORM, hrz);
138    }
139
140    RstTime GetLunarRst(const double &jd, const LnLatPosn &obs=kORM)
141    {
142        RstTime rst;
143        ln_get_lunar_rst(jd, const_cast<LnLatPosn*>(&obs), &rst);
144        return rst;
145    }
146    EquPosn GetSolarEquCoords(const double &jd)
147    {
148        EquPosn equ;
149        ln_get_solar_equ_coords(jd, &equ);
150        return equ;
151    }
152
153    double GetLunarDisk(const double &jd)
154    {
155        return ln_get_lunar_disk(jd);
156    }
157
158    double GetLunarSdiam(const double &jd)
159    {
160        return ln_get_lunar_sdiam(jd);
161    }
162
163    double GetLunarPhase(const double &jd)
164    {
165        return ln_get_lunar_phase(jd);
166    }
167
168    EquPosn GetLunarEquCoords(const double &jd, double precision=0)
169    {
170        EquPosn equ;
171        ln_get_lunar_equ_coords_prec(jd, &equ, precision);
172        return equ;
173    }
174
175    double GetLunarEarthDist(const double &jd)
176    {
177        return ln_get_lunar_earth_dist(jd);
178    }
179
180    double GetAngularSeparation(const EquPosn &p1, const EquPosn &p2)
181    {
182        return ln_get_angular_separation(const_cast<EquPosn*>(&p1), const_cast<EquPosn*>(&p2));
183    }
184
185    double GetAngularSeparation(const HrzPosn &h1, const HrzPosn &h2)
186    {
187        EquPosn p1; p1.ra=h1.az; p1.dec=h1.alt;
188        EquPosn p2; p2.ra=h2.az; p2.dec=h2.alt;
189        return ln_get_angular_separation(&p1, &p2);
190    }
191
192    struct SolarObjects
193    {
194        double  fJD;
195
196        EquPosn fSunEqu;
197        HrzPosn fSunHrz;
198
199        EquPosn fMoonEqu;
200        HrzPosn fMoonHrz;
201        double  fMoonDisk;
202
203        double fEarthDist;
204
205        SolarObjects(const double &jd, const LnLatPosn &obs=Nova::kORM)
206        {
207            fJD = jd;
208
209            // Sun properties
210            fSunEqu    = GetSolarEquCoords(jd);
211            fSunHrz    = GetHrzFromEqu(fSunEqu, obs, jd);
212
213            // Moon properties
214            fMoonEqu   = GetLunarEquCoords(jd, 0.01);
215            fMoonHrz   = GetHrzFromEqu(fMoonEqu, obs, jd);
216
217            fMoonDisk  = GetLunarDisk(jd);
218            fEarthDist = GetLunarEarthDist(jd);
219        }
220    };
221
222}
223
224#endif
Note: See TracBrowser for help on using the repository browser.