source: trunk/FACT++/src/moon.cc@ 15658

Last change on this file since 15658 was 14308, checked in by tbretz, 12 years ago
Added the distance to the sun... the sun's zd needs to be checked... it seems it is not correctly calculated below the horizon.
File size: 5.1 KB
Line 
1#include <libnova/solar.h>
2#include <libnova/lunar.h>
3#include <libnova/rise_set.h>
4#include <libnova/transform.h>
5
6#include "Time.h"
7#include "Configuration.h"
8
9using namespace std;
10
11// ------------------------------------------------------------------------
12
13class Moon
14{
15public:
16 Time time;
17
18 double ra;
19 double dec;
20
21 double zd;
22 double az;
23
24 double disk;
25
26 bool visible;
27
28 Time fRise;
29 Time fTransit;
30 Time fSet;
31
32 int state;
33
34 Moon() : time(Time::none)
35 {
36 }
37
38 // Could be done more efficient: Only recalcuate if
39 // the current time exceeds at least on of the stored times
40 Moon(double lon, double lat, const Time &t=Time()) : time(t)
41 {
42 const double JD = time.JD();
43
44 ln_lnlat_posn observer;
45 observer.lng = lon;
46 observer.lat = lat;
47
48 ln_rst_time moon;
49 ln_get_lunar_rst(JD-0.5, &observer, &moon);
50
51 fRise = Time(moon.rise);
52 fTransit = Time(moon.transit);
53 fSet = Time(moon.set);
54
55 //visible =
56 // ((JD>moon.rise && JD<moon.set ) && moon.rise<moon.set) ||
57 // ((JD<moon.set || JD>moon.rise) && moon.rise>moon.set);
58
59 const bool is_up = JD>moon.rise;
60 const bool is_sinking = JD>moon.transit;
61 const bool is_dn = JD>moon.set;
62
63 ln_get_lunar_rst(JD+0.5, &observer, &moon);
64 if (is_up)
65 fRise = Time(moon.rise);
66 if (is_sinking)
67 fTransit = Time(moon.transit);
68 if (is_dn)
69 fSet = Time(moon.set);
70
71 ln_equ_posn pos;
72 ln_get_lunar_equ_coords(JD, &pos);
73
74 ln_hrz_posn hrz;
75 ln_get_hrz_from_equ (&pos, &observer, JD, &hrz);
76 az = hrz.az;
77 zd = 90-hrz.alt;
78
79 ra = pos.ra/15;
80 dec = pos.dec;
81
82 disk = ln_get_lunar_disk(JD)*100;
83 state = 0;
84 if (fRise <fTransit && fRise <fSet) state = 0; // not visible
85 if (fTransit<fSet && fTransit<fRise) state = 1; // before culm
86 if (fSet <fRise && fSet <fTransit) state = 2; // after culm
87
88 visible = state!=0;
89
90 // 0: not visible
91 // 1: visible before cul
92 // 2: visible after cul
93 }
94
95 double Angle(double r, double d) const
96 {
97 const double theta0 = M_PI/2-d*M_PI/180;
98 const double phi0 = r*M_PI/12;
99
100 const double theta1 = M_PI/2-dec*M_PI/180;
101 const double phi1 = ra*M_PI/12;
102
103 const double x0 = sin(theta0) * cos(phi0);
104 const double y0 = sin(theta0) * sin(phi0);
105 const double z0 = cos(theta0);
106
107 const double x1 = sin(theta1) * cos(phi1);
108 const double y1 = sin(theta1) * sin(phi1);
109 const double z1 = cos(theta1);
110
111 double arg = x0*x1 + y0*y1 + z0*z1;
112 if(arg > 1.0) arg = 1.0;
113 if(arg < -1.0) arg = -1.0;
114
115 return acos(arg) * 180/M_PI;
116 }
117};
118
119// ========================================================================
120// ========================================================================
121// ========================================================================
122
123void SetupConfiguration(Configuration &conf)
124{
125 po::options_description control("Smart FACT");
126 control.add_options()
127 ("ra", var<double>(), "Source right ascension")
128 ("dec", var<double>(), "Source declination")
129 ("date-time", var<string>()
130#if BOOST_VERSION >= 104200
131 ->required()
132#endif
133 , "SQL time (UTC)")
134 ;
135
136 po::positional_options_description p;
137 p.add("date-time", 1); // The first positional options
138
139 conf.AddOptions(control);
140 conf.SetArgumentPositions(p);
141}
142
143void PrintUsage()
144{
145 cout <<
146 "moon - The moon calculator\n"
147 "\n"
148 "Usage: moon sql-datetime [--ra={ra} --dec={dec}]\n";
149 cout << endl;
150}
151
152int main(int argc, const char* argv[])
153{
154 Configuration conf(argv[0]);
155 conf.SetPrintUsage(PrintUsage);
156 SetupConfiguration(conf);
157
158 if (!conf.DoParse(argc, argv))
159 return 127;
160
161 if (conf.Has("ra")^conf.Has("dec"))
162 {
163 cout << "ERROR - Either --ra or --dec missing." << endl;
164 return 1;
165 }
166
167 Time time;
168 time.SetFromStr(conf.Get<string>("date-time"));
169
170 ln_lnlat_posn observer;
171 observer.lng = -(17.+53./60+26.525/3600);
172 observer.lat = 28.+45./60+42.462/3600;
173
174 Moon moon(observer.lng, observer.lat, time);
175
176 cout << setprecision(15);
177 cout << time.GetAsStr() << '\n';
178
179 ln_equ_posn pos;
180 ln_hrz_posn hrz;
181 ln_get_solar_equ_coords(time.JD(), &pos);
182 ln_get_hrz_from_equ(&pos, &observer, time.JD(), &hrz);
183 cout << 90-hrz.alt << '\n';
184
185 cout << moon.visible << '\n';
186 cout << moon.disk << '\n';
187 cout << moon.zd << '\n';
188
189 if (conf.Has("ra") && conf.Has("dec"))
190 {
191 pos.ra = conf.Get<double>("ra")*15;
192 pos.dec = conf.Get<double>("dec");
193
194 cout << moon.Angle(pos.ra/15, pos.dec) << '\n';
195
196 // Trick 17
197 moon.ra = pos.ra;
198 moon.dec = pos.dec;
199
200 // Sun distance
201 cout << moon.Angle(pos.ra/15, pos.dec) << '\n';
202 }
203
204 cout << endl;
205
206 return 0;
207}
Note: See TracBrowser for help on using the repository browser.