source: trunk/external-alma/atnf/pks/pks_maths.cc@ 3146

Last change on this file since 3146 was 3029, checked in by Kana Sugimoto, 10 years ago

New Development: Yes

JIRA Issue: Yes (CAS-6929)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): asap as a whole

Description: committing Darrell's changes to make asap work with merged casacore.


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 <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
54Int 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
64Double 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
73Double 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
82Float 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
147Double 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
155void 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
181void 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
229void 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
263void 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
305void 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
335Double 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}
Note: See TracBrowser for help on using the repository browser.