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

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