1 | //#---------------------------------------------------------------------------
|
---|
2 | //# pks_maths.cc: Mathematical functions for Parkes single-dish data reduction
|
---|
3 | //#---------------------------------------------------------------------------
|
---|
4 | //# livedata - processing pipeline for single-dish, multibeam spectral data.
|
---|
5 | //# Copyright (C) 2004-2009, Australia Telescope National Facility, CSIRO
|
---|
6 | //#
|
---|
7 | //# This file is part of livedata.
|
---|
8 | //#
|
---|
9 | //# livedata is free software: you can redistribute it and/or modify it under
|
---|
10 | //# the terms of the GNU General Public License as published by the Free
|
---|
11 | //# Software Foundation, either version 3 of the License, or (at your option)
|
---|
12 | //# any later version.
|
---|
13 | //#
|
---|
14 | //# livedata is distributed in the hope that it will be useful, but WITHOUT
|
---|
15 | //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
---|
16 | //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
|
---|
17 | //# more details.
|
---|
18 | //#
|
---|
19 | //# You should have received a copy of the GNU General Public License along
|
---|
20 | //# with livedata. If not, see <http://www.gnu.org/licenses/>.
|
---|
21 | //#
|
---|
22 | //# Correspondence concerning livedata may be directed to:
|
---|
23 | //# Internet email: mcalabre@atnf.csiro.au
|
---|
24 | //# Postal address: Dr. Mark Calabretta
|
---|
25 | //# Australia Telescope National Facility, CSIRO
|
---|
26 | //# PO Box 76
|
---|
27 | //# Epping NSW 1710
|
---|
28 | //# AUSTRALIA
|
---|
29 | //#
|
---|
30 | //# http://www.atnf.csiro.au/computing/software/livedata.html
|
---|
31 | //# $Id: pks_maths.cc,v 1.7 2009-09-29 07:45:02 cal103 Exp $
|
---|
32 | //#---------------------------------------------------------------------------
|
---|
33 | //# Original: 2004/07/16 Mark Calabretta
|
---|
34 | //#---------------------------------------------------------------------------
|
---|
35 |
|
---|
36 | // AIPS++ includes.
|
---|
37 | #include <complex>
|
---|
38 | #include <casa/aips.h>
|
---|
39 | #include <casa/math.h>
|
---|
40 | #include <casa/Arrays/ArrayMath.h>
|
---|
41 | #include <casa/Arrays/Vector.h>
|
---|
42 | #include <casa/BasicSL/Constants.h>
|
---|
43 | #include <casa/Utilities/GenSort.h>
|
---|
44 |
|
---|
45 | // Parkes includes.
|
---|
46 | #include <atnf/pks/pks_maths.h>
|
---|
47 |
|
---|
48 |
|
---|
49 | //----------------------------------------------------------------------- nint
|
---|
50 |
|
---|
51 | // Nearest integral value; halfway cases are rounded to the integral value
|
---|
52 | // larger in value. No check is made for integer overflow.
|
---|
53 |
|
---|
54 | Int nint(Double v)
|
---|
55 | {
|
---|
56 | return Int(floor(v + 0.5));
|
---|
57 | }
|
---|
58 |
|
---|
59 | //---------------------------------------------------------------------- anint
|
---|
60 |
|
---|
61 | // Nearest integral value; halfway cases are rounded to the integral value
|
---|
62 | // larger in value.
|
---|
63 |
|
---|
64 | Double anint(Double v)
|
---|
65 | {
|
---|
66 | return floor(v + 0.5);
|
---|
67 | }
|
---|
68 |
|
---|
69 | //---------------------------------------------------------------------- round
|
---|
70 |
|
---|
71 | // Round value v to the nearest integral multiple of precision p.
|
---|
72 |
|
---|
73 | Double round(Double v, Double p)
|
---|
74 | {
|
---|
75 | return p * floor(v/p + 0.5);
|
---|
76 | }
|
---|
77 |
|
---|
78 | //--------------------------------------------------------------------- median
|
---|
79 |
|
---|
80 | // Compute the weighted median value of an array.
|
---|
81 |
|
---|
82 | Float median(const Vector<Float> &v, const Vector<Float> &wgt)
|
---|
83 | {
|
---|
84 | uInt nElem = v.nelements();
|
---|
85 | if (nElem == 0) return 0.0f;
|
---|
86 |
|
---|
87 | // Generate the sort index.
|
---|
88 | Vector<uInt> sortindex(nElem);
|
---|
89 | GenSortIndirect<Float>::sort(sortindex, v);
|
---|
90 |
|
---|
91 | // Find the middle weight.
|
---|
92 | Float wgt_2 = sum(wgt)/2.0f;
|
---|
93 |
|
---|
94 | // Find the corresponding vector element.
|
---|
95 | Float weight = 0.0f;
|
---|
96 | Float accwgt = 0.0f;
|
---|
97 | uInt j1 = 0;
|
---|
98 | uInt j2;
|
---|
99 | for (j2 = 0; j2 < nElem; j2++) {
|
---|
100 | weight = wgt(sortindex(j2));
|
---|
101 | if (weight == 0.0f) {
|
---|
102 | // Ignore zero-weight data;
|
---|
103 | continue;
|
---|
104 | }
|
---|
105 |
|
---|
106 | // The accumulated weight.
|
---|
107 | accwgt += weight;
|
---|
108 |
|
---|
109 | if (accwgt <= wgt_2) {
|
---|
110 | // Keep looping.
|
---|
111 | j1 = j2;
|
---|
112 | } else {
|
---|
113 | break;
|
---|
114 | }
|
---|
115 | }
|
---|
116 |
|
---|
117 | // Compute weighted median.
|
---|
118 | Float v1 = v(sortindex(j1));
|
---|
119 | Float v2 = v(sortindex(j2));
|
---|
120 |
|
---|
121 | // Compute pro-rata value from below.
|
---|
122 | Float dw = wgt_2 - (accwgt - weight);
|
---|
123 | v1 += (v2 - v1) * dw / weight;
|
---|
124 |
|
---|
125 | // Find next non-zero-weight value.
|
---|
126 | for (j2++ ; j2 < nElem; j2++) {
|
---|
127 | weight = wgt(sortindex(j2));
|
---|
128 | if (weight != 0.0f) {
|
---|
129 | break;
|
---|
130 | }
|
---|
131 | }
|
---|
132 |
|
---|
133 | if (j2 < nElem) {
|
---|
134 | // Compute pro-rata value from above.
|
---|
135 | Float v3 = v(sortindex(j2));
|
---|
136 |
|
---|
137 | v2 += (v3 - v2) * dw / weight;
|
---|
138 | }
|
---|
139 |
|
---|
140 | return (v1 + v2)/2.0f;
|
---|
141 | }
|
---|
142 |
|
---|
143 | //---------------------------------------------------------------- angularDist
|
---|
144 |
|
---|
145 | // Determine the angular distance between two directions (angles in radians).
|
---|
146 |
|
---|
147 | Double angularDist(Double lng0, Double lat0, Double lng, Double lat)
|
---|
148 | {
|
---|
149 | Double costheta = sin(lat0)*sin(lat) + cos(lat0)*cos(lat)*cos(lng0-lng);
|
---|
150 | return acos(costheta);
|
---|
151 | }
|
---|
152 |
|
---|
153 | //--------------------------------------------------------------------- distPA
|
---|
154 |
|
---|
155 | void distPA(Double lng0, Double lat0, Double lng, Double lat, Double &dist,
|
---|
156 | Double &pa)
|
---|
157 |
|
---|
158 | // Determine the generalized position angle of the field point (lng,lat) from
|
---|
159 | // the reference point (lng0,lat0) and the angular distance between them
|
---|
160 | // (angles in radians).
|
---|
161 |
|
---|
162 | {
|
---|
163 | // Euler angles which rotate the coordinate frame so that (lng0,lat0) is
|
---|
164 | // at the pole of the new system, with the pole of the old system at zero
|
---|
165 | // longitude in the new.
|
---|
166 | Double phi0 = C::pi_2 + lng0;
|
---|
167 | Double theta = C::pi_2 - lat0;
|
---|
168 | Double phi = -C::pi_2;
|
---|
169 |
|
---|
170 | // Rotate the field point to the new system.
|
---|
171 | Double alpha, beta;
|
---|
172 | eulerx(lng, lat, phi0, theta, phi, alpha, beta);
|
---|
173 |
|
---|
174 | dist = C::pi_2 - beta;
|
---|
175 | pa = -alpha;
|
---|
176 | if (pa < -C::pi) pa = pa + C::_2pi;
|
---|
177 | }
|
---|
178 |
|
---|
179 | //--------------------------------------------------------------------- eulerx
|
---|
180 |
|
---|
181 | void eulerx(Double lng0, Double lat0, Double phi0, Double theta, Double phi,
|
---|
182 | Double &lng1, Double &lat1)
|
---|
183 |
|
---|
184 | // Applies the Euler angle based transformation of spherical coordinates.
|
---|
185 | //
|
---|
186 | // phi0 Longitude of the ascending node in the old system, radians. The
|
---|
187 | // ascending node is the point of intersection of the equators of
|
---|
188 | // the two systems such that the equator of the new system crosses
|
---|
189 | // from south to north as viewed in the old system.
|
---|
190 | //
|
---|
191 | // theta Angle between the poles of the two systems, radians. THETA is
|
---|
192 | // positive for a positive rotation about the ascending node.
|
---|
193 | //
|
---|
194 | // phi Longitude of the ascending node in the new system, radians.
|
---|
195 |
|
---|
196 | {
|
---|
197 | // Compute intermediaries.
|
---|
198 | Double lng0p = lng0 - phi0;
|
---|
199 | Double slng0p = sin(lng0p);
|
---|
200 | Double clng0p = cos(lng0p);
|
---|
201 | Double slat0 = sin(lat0);
|
---|
202 | Double clat0 = cos(lat0);
|
---|
203 | Double ctheta = cos(theta);
|
---|
204 | Double stheta = sin(theta);
|
---|
205 |
|
---|
206 | Double x = clat0*clng0p;
|
---|
207 | Double y = clat0*slng0p*ctheta + slat0*stheta;
|
---|
208 |
|
---|
209 | // Longitude in the new system.
|
---|
210 | if (x != 0.0 || y != 0.0) {
|
---|
211 | lng1 = phi + atan2(y, x);
|
---|
212 | } else {
|
---|
213 | // Longitude at the poles in the new system is consistent with that
|
---|
214 | // specified in the old system.
|
---|
215 | lng1 = phi + lng0p;
|
---|
216 | }
|
---|
217 | lng1 = fmod(lng1, C::_2pi);
|
---|
218 | if (lng1 < 0.0) lng1 += C::_2pi;
|
---|
219 |
|
---|
220 | lat1 = asin(slat0*ctheta - clat0*stheta*slng0p);
|
---|
221 | }
|
---|
222 |
|
---|
223 | //------------------------------------------------------------------------ sol
|
---|
224 |
|
---|
225 | // Low precision coordinates of the Sun (accurate to 1 arcmin between 1800 and
|
---|
226 | // 2200) from http://aa.usno.navy.mil/faq/docs/SunApprox.html matches closely
|
---|
227 | // that in the Astronomical Almanac.
|
---|
228 |
|
---|
229 | void sol(Double mjd, Double &elng, Double &ra, Double &dec)
|
---|
230 | {
|
---|
231 | Double d2r = C::pi/180.0;
|
---|
232 |
|
---|
233 | // Number of days since J2000.0.
|
---|
234 | Double d = mjd - 51544.5;
|
---|
235 |
|
---|
236 | // Mean longitude and mean anomaly of the Sun (deg).
|
---|
237 | Double L = 280.459 + 0.98564736*d;
|
---|
238 | Double g = 357.529 + 0.98560028*d;
|
---|
239 |
|
---|
240 | // Apparent ecliptic longitude corrected for aberration (deg).
|
---|
241 | g *= d2r;
|
---|
242 | elng = L + 1.915*sin(g) + 0.020*sin(g+g);
|
---|
243 | elng = fmod(elng, 360.0);
|
---|
244 | if (elng < 0.0) elng += 360.0;
|
---|
245 |
|
---|
246 | // Obliquity of the ecliptic (deg).
|
---|
247 | Double epsilon = 23.439 - 0.00000036*d;
|
---|
248 |
|
---|
249 | // Transform ecliptic to equatorial coordinates.
|
---|
250 | elng *= d2r;
|
---|
251 | epsilon *= d2r;
|
---|
252 | ra = atan2(cos(epsilon)*sin(elng), cos(elng));
|
---|
253 | dec = asin(sin(epsilon)*sin(elng));
|
---|
254 | if (ra < 0.0) ra += C::_2pi;
|
---|
255 | }
|
---|
256 |
|
---|
257 | //------------------------------------------------------------------------ gst
|
---|
258 |
|
---|
259 | // Greenwich mean sidereal time, and low precision Greenwich apparent sidereal
|
---|
260 | // time, both in radian, from http://aa.usno.navy.mil/faq/docs/GAST.html. UT1
|
---|
261 | // is given in MJD form.
|
---|
262 |
|
---|
263 | void gst(Double ut1, Double &gmst, Double &gast)
|
---|
264 | {
|
---|
265 | Double d2r = C::pi/180.0;
|
---|
266 |
|
---|
267 | Double d = ut1 - 51544.5;
|
---|
268 | Double d0 = int(ut1) - 51544.5;
|
---|
269 | Double h = 24.0*(d - d0);
|
---|
270 | Double t = d / 35625.0;
|
---|
271 |
|
---|
272 | // GMST (hr).
|
---|
273 | gmst = 6.697374558 + 0.06570982441908*d0 + 1.00273790935*h + 0.000026*t*t;
|
---|
274 | gmst = fmod(gmst, 24.0);
|
---|
275 |
|
---|
276 | // Longitude of the ascending node of the Moon (deg).
|
---|
277 | Double Omega = 125.04 - 0.052954*d;
|
---|
278 |
|
---|
279 | // Mean Longitude of the Sun (deg).
|
---|
280 | Double L = 280.47 + 0.98565*d;
|
---|
281 |
|
---|
282 | // Obliquity of the ecliptic (deg).
|
---|
283 | Double epsilon = 23.4393 - 0.0000004*d;
|
---|
284 |
|
---|
285 | // Approximate nutation in longitude (hr).
|
---|
286 | Double dpsi = -0.000319*sin(Omega*d2r) - 0.000024*sin((L+L)*d2r);
|
---|
287 |
|
---|
288 | // Equation of the equinoxes (hr).
|
---|
289 | Double eqeq = dpsi*cos(epsilon*d2r);
|
---|
290 |
|
---|
291 | // GAST (hr).
|
---|
292 | gast = gmst + eqeq;
|
---|
293 | gast = fmod(gast, 24.0);
|
---|
294 |
|
---|
295 | // Convert to radian.
|
---|
296 | gmst *= C::pi/12.0;
|
---|
297 | gast *= C::pi/12.0;
|
---|
298 | }
|
---|
299 |
|
---|
300 | //----------------------------------------------------------------------- azel
|
---|
301 |
|
---|
302 | // Convert (ra,dec) to (az,el). Position as a Cartesian triplet in m, UT1 in
|
---|
303 | // MJD form, and all angles in radian.
|
---|
304 |
|
---|
305 | void azel(const Vector<Double> position, Double ut1, Double ra, Double dec,
|
---|
306 | Double &az, Double &el)
|
---|
307 | {
|
---|
308 | // Get geocentric longitude and latitude (rad).
|
---|
309 | Double x = position(0);
|
---|
310 | Double y = position(1);
|
---|
311 | Double z = position(2);
|
---|
312 | Double r = sqrt(x*x + y*y + z*z);
|
---|
313 | Double lng = atan2(y, x);
|
---|
314 | Double lat = asin(z/r);
|
---|
315 |
|
---|
316 | // Get GAST (rad).
|
---|
317 | Double gast, gmst;
|
---|
318 | gst(ut1, gmst, gast);
|
---|
319 |
|
---|
320 | // Local hour angle (rad).
|
---|
321 | Double ha = (gast + lng) - ra;
|
---|
322 |
|
---|
323 | // Azimuth and elevation (rad).
|
---|
324 | az = atan2(-cos(dec)*sin(ha),
|
---|
325 | sin(dec)*cos(lat) - cos(dec)*sin(lat)*cos(ha));
|
---|
326 | el = asin(sin(dec)*sin(lat) + cos(dec)*cos(lat)*cos(ha));
|
---|
327 |
|
---|
328 | if (az < 0.0) az += C::_2pi;
|
---|
329 | }
|
---|
330 |
|
---|
331 | //---------------------------------------------------------------------- solel
|
---|
332 |
|
---|
333 | // Compute the Solar elevation using the above functions.
|
---|
334 |
|
---|
335 | Double solel(const Vector<Double> position, Double ut1)
|
---|
336 | {
|
---|
337 | Double az, dec, el, elng, gast, gmst, ra;
|
---|
338 | sol(ut1, elng, ra, dec);
|
---|
339 | gst(ut1, gmst, gast);
|
---|
340 | azel(position, ut1, ra, dec, az, el);
|
---|
341 | return el;
|
---|
342 | }
|
---|