source: trunk/external/atnf/pks/pks_maths.cc@ 2571

Last change on this file since 2571 was 1720, checked in by Malte Marquarding, 15 years ago

Update from livedata CVS repository

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