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

Last change on this file since 1449 was 1427, checked in by Malte Marquarding, 16 years ago

sync with livedata/implement/atnf

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