1 /*============================================================================
2 WCSLIB 7.3 - an implementation of the FITS WCS standard.
3 Copyright (C) 1995-2020, Mark Calabretta
4
5 This file is part of WCSLIB.
6
7 WCSLIB is free software: you can redistribute it and/or modify it under the
8 terms of the GNU Lesser General Public License as published by the Free
9 Software Foundation, either version 3 of the License, or (at your option)
10 any later version.
11
12 WCSLIB is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with WCSLIB. If not, see http://www.gnu.org/licenses.
19
20 Direct correspondence concerning WCSLIB to mark@calabretta.id.au
21
22 Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
23 http://www.atnf.csiro.au/people/Mark.Calabretta
24 $Id: sph.c,v 7.3.1.2 2020/08/17 11:19:09 mcalabre Exp mcalabre $
25 *===========================================================================*/
26
27 #include <math.h>
28 #include "wcstrig.h"
29 #include "sph.h"
30
31 #define copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
32
33 #define tol 1.0e-5
34
35 //----------------------------------------------------------------------------
36
sphx2s(const double eul[5],int nphi,int ntheta,int spt,int sll,const double phi[],const double theta[],double lng[],double lat[])37 int sphx2s(
38 const double eul[5],
39 int nphi,
40 int ntheta,
41 int spt,
42 int sll,
43 const double phi[],
44 const double theta[],
45 double lng[],
46 double lat[])
47
48 {
49 int jphi, mphi, mtheta, rowlen, rowoff;
50 double cosphi, costhe, costhe3, costhe4, dlng, dphi, sinphi, sinthe,
51 sinthe3, sinthe4, x, y, z;
52 register int iphi, itheta;
53 register const double *phip, *thetap;
54 register double *latp, *lngp;
55
56 if (ntheta > 0) {
57 mphi = nphi;
58 mtheta = ntheta;
59 } else {
60 mphi = 1;
61 mtheta = 1;
62 ntheta = nphi;
63 }
64
65
66 // Check for special-case rotations.
67 if (eul[4] == 0.0) {
68 if (eul[1] == 0.0) {
69 // Simple change in origin of longitude.
70 dlng = fmod(eul[0] + 180.0 - eul[2], 360.0);
71
72 jphi = 0;
73 thetap = theta;
74 lngp = lng;
75 latp = lat;
76 for (itheta = 0; itheta < ntheta; itheta++, thetap += spt) {
77 phip = phi + (jphi%nphi)*spt;
78 for (iphi = 0; iphi < mphi; iphi++, phip += spt, jphi++) {
79 *lngp = *phip + dlng;
80 *latp = *thetap;
81
82 // Normalize the celestial longitude.
83 if (eul[0] >= 0.0) {
84 if (*lngp < 0.0) *lngp += 360.0;
85 } else {
86 if (*lngp > 0.0) *lngp -= 360.0;
87 }
88
89 if (*lngp > 360.0) {
90 *lngp -= 360.0;
91 } else if (*lngp < -360.0) {
92 *lngp += 360.0;
93 }
94
95 lngp += sll;
96 latp += sll;
97 }
98 }
99
100 } else {
101 // Pole-flip with change in origin of longitude.
102 dlng = fmod(eul[0] + eul[2], 360.0);
103
104 jphi = 0;
105 thetap = theta;
106 lngp = lng;
107 latp = lat;
108 for (itheta = 0; itheta < ntheta; itheta++, thetap += spt) {
109 phip = phi + (jphi%nphi)*spt;
110 for (iphi = 0; iphi < mphi; iphi++, phip += spt, jphi++) {
111 *lngp = dlng - *phip;
112 *latp = -(*thetap);
113
114 // Normalize the celestial longitude.
115 if (eul[0] >= 0.0) {
116 if (*lngp < 0.0) *lngp += 360.0;
117 } else {
118 if (*lngp > 0.0) *lngp -= 360.0;
119 }
120
121 if (*lngp > 360.0) {
122 *lngp -= 360.0;
123 } else if (*lngp < -360.0) {
124 *lngp += 360.0;
125 }
126
127 lngp += sll;
128 latp += sll;
129 }
130 }
131 }
132
133 return 0;
134 }
135
136
137 // Do phi dependency.
138 phip = phi;
139 rowoff = 0;
140 rowlen = nphi*sll;
141 for (iphi = 0; iphi < nphi; iphi++, rowoff += sll, phip += spt) {
142 dphi = *phip - eul[2];
143
144 lngp = lng + rowoff;
145 for (itheta = 0; itheta < mtheta; itheta++) {
146 *lngp = dphi;
147 lngp += rowlen;
148 }
149 }
150
151
152 // Do theta dependency.
153 thetap = theta;
154 lngp = lng;
155 latp = lat;
156 for (itheta = 0; itheta < ntheta; itheta++, thetap += spt) {
157 sincosd(*thetap, &sinthe, &costhe);
158 costhe3 = costhe*eul[3];
159 costhe4 = costhe*eul[4];
160 sinthe3 = sinthe*eul[3];
161 sinthe4 = sinthe*eul[4];
162
163 for (iphi = 0; iphi < mphi; iphi++, lngp += sll, latp += sll) {
164 dphi = *lngp;
165 sincosd(dphi, &sinphi, &cosphi);
166
167 // Compute the celestial longitude.
168 x = sinthe4 - costhe3*cosphi;
169 if (fabs(x) < tol) {
170 // Rearrange formula to reduce roundoff errors.
171 x = -cosd(*thetap + eul[1]) + costhe3*(1.0 - cosphi);
172 }
173
174 y = -costhe*sinphi;
175 if (x != 0.0 || y != 0.0) {
176 dlng = atan2d(y, x);
177 } else {
178 // Change of origin of longitude.
179 if (eul[1] < 90.0) {
180 dlng = dphi + 180.0;
181 } else {
182 dlng = -dphi;
183 }
184 }
185 *lngp = eul[0] + dlng;
186
187 // Normalize the celestial longitude.
188 if (eul[0] >= 0.0) {
189 if (*lngp < 0.0) *lngp += 360.0;
190 } else {
191 if (*lngp > 0.0) *lngp -= 360.0;
192 }
193
194 if (*lngp > 360.0) {
195 *lngp -= 360.0;
196 } else if (*lngp < -360.0) {
197 *lngp += 360.0;
198 }
199
200 // Compute the celestial latitude.
201 if (fmod(dphi,180.0) == 0.0) {
202 *latp = *thetap + cosphi*eul[1];
203 if (*latp > 90.0) *latp = 180.0 - *latp;
204 if (*latp < -90.0) *latp = -180.0 - *latp;
205 } else {
206 z = sinthe3 + costhe4*cosphi;
207 if (fabs(z) > 0.99) {
208 // Use an alternative formula for greater accuracy.
209 *latp = copysign(acosd(sqrt(x*x+y*y)), z);
210 } else {
211 *latp = asind(z);
212 }
213 }
214 }
215 }
216
217 return 0;
218 }
219
220 //----------------------------------------------------------------------------
221
sphs2x(const double eul[5],int nlng,int nlat,int sll,int spt,const double lng[],const double lat[],double phi[],double theta[])222 int sphs2x(
223 const double eul[5],
224 int nlng,
225 int nlat,
226 int sll,
227 int spt,
228 const double lng[],
229 const double lat[],
230 double phi[],
231 double theta[])
232
233 {
234 int jlng, mlat, mlng, rowlen, rowoff;
235 double coslat, coslat3, coslat4, coslng, dlng, dphi, sinlat, sinlat3,
236 sinlat4, sinlng, x, y, z;
237 register int ilat, ilng;
238 register const double *latp, *lngp;
239 register double *phip, *thetap;
240
241 if (nlat > 0) {
242 mlng = nlng;
243 mlat = nlat;
244 } else {
245 mlng = 1;
246 mlat = 1;
247 nlat = nlng;
248 }
249
250
251 // Check for special-case rotations.
252 if (eul[4] == 0.0) {
253 if (eul[1] == 0.0) {
254 // Simple change in origin of longitude.
255 dphi = fmod(eul[2] - 180.0 - eul[0], 360.0);
256
257 jlng = 0;
258 latp = lat;
259 phip = phi;
260 thetap = theta;
261 for (ilat = 0; ilat < nlat; ilat++, latp += sll) {
262 lngp = lng + (jlng%nlng)*sll;
263 for (ilng = 0; ilng < mlng; ilng++, lngp += sll, jlng++) {
264 *phip = fmod(*lngp + dphi, 360.0);
265 *thetap = *latp;
266
267 // Normalize the native longitude.
268 if (*phip > 180.0) {
269 *phip -= 360.0;
270 } else if (*phip < -180.0) {
271 *phip += 360.0;
272 }
273
274 phip += spt;
275 thetap += spt;
276 }
277 }
278
279 } else {
280 // Pole-flip with change in origin of longitude.
281 dphi = fmod(eul[2] + eul[0], 360.0);
282
283 jlng = 0;
284 latp = lat;
285 phip = phi;
286 thetap = theta;
287 for (ilat = 0; ilat < nlat; ilat++, latp += sll) {
288 lngp = lng + (jlng%nlng)*sll;
289 for (ilng = 0; ilng < mlng; ilng++, lngp += sll, jlng++) {
290 *phip = fmod(dphi - *lngp, 360.0);
291 *thetap = -(*latp);
292
293 // Normalize the native longitude.
294 if (*phip > 180.0) {
295 *phip -= 360.0;
296 } else if (*phip < -180.0) {
297 *phip += 360.0;
298 }
299
300 phip += spt;
301 thetap += spt;
302 }
303 }
304 }
305
306 return 0;
307 }
308
309
310 // Do lng dependency.
311 lngp = lng;
312 rowoff = 0;
313 rowlen = nlng*spt;
314 for (ilng = 0; ilng < nlng; ilng++, rowoff += spt, lngp += sll) {
315 dlng = *lngp - eul[0];
316
317 phip = phi + rowoff;
318 thetap = theta;
319 for (ilat = 0; ilat < mlat; ilat++) {
320 *phip = dlng;
321 phip += rowlen;
322 }
323 }
324
325
326 // Do lat dependency.
327 latp = lat;
328 phip = phi;
329 thetap = theta;
330 for (ilat = 0; ilat < nlat; ilat++, latp += sll) {
331 sincosd(*latp, &sinlat, &coslat);
332 coslat3 = coslat*eul[3];
333 coslat4 = coslat*eul[4];
334 sinlat3 = sinlat*eul[3];
335 sinlat4 = sinlat*eul[4];
336
337 for (ilng = 0; ilng < mlng; ilng++, phip += spt, thetap += spt) {
338 dlng = *phip;
339 sincosd(dlng, &sinlng, &coslng);
340
341 // Compute the native longitude.
342 x = sinlat4 - coslat3*coslng;
343 if (fabs(x) < tol) {
344 // Rearrange formula to reduce roundoff errors.
345 x = -cosd(*latp+eul[1]) + coslat3*(1.0 - coslng);
346 }
347
348 y = -coslat*sinlng;
349 if (x != 0.0 || y != 0.0) {
350 dphi = atan2d(y, x);
351 } else {
352 // Change of origin of longitude.
353 if (eul[1] < 90.0) {
354 dphi = dlng - 180.0;
355 } else {
356 dphi = -dlng;
357 }
358 }
359 *phip = fmod(eul[2] + dphi, 360.0);
360
361 // Normalize the native longitude.
362 if (*phip > 180.0) {
363 *phip -= 360.0;
364 } else if (*phip < -180.0) {
365 *phip += 360.0;
366 }
367
368 // Compute the native latitude.
369 if (fmod(dlng,180.0) == 0.0) {
370 *thetap = *latp + coslng*eul[1];
371 if (*thetap > 90.0) *thetap = 180.0 - *thetap;
372 if (*thetap < -90.0) *thetap = -180.0 - *thetap;
373 } else {
374 z = sinlat3 + coslat4*coslng;
375 if (fabs(z) > 0.99) {
376 // Use an alternative formula for greater accuracy.
377 *thetap = copysign(acosd(sqrt(x*x+y*y)), z);
378 } else {
379 *thetap = asind(z);
380 }
381 }
382 }
383 }
384
385 return 0;
386 }
387
388 //----------------------------------------------------------------------------
389
sphdpa(int nfield,double lng0,double lat0,const double lng[],const double lat[],double dist[],double pa[])390 int sphdpa(
391 int nfield,
392 double lng0,
393 double lat0,
394 const double lng[],
395 const double lat[],
396 double dist[],
397 double pa[])
398
399 {
400 int i;
401 double eul[5];
402
403 // Set the Euler angles for the coordinate transformation.
404 eul[0] = lng0;
405 eul[1] = 90.0 - lat0;
406 eul[2] = 0.0;
407 eul[3] = cosd(eul[1]);
408 eul[4] = sind(eul[1]);
409
410 // Transform field points to the new system.
411 sphs2x(eul, nfield, 0, 1, 1, lng, lat, pa, dist);
412
413 for (i = 0; i < nfield; i++) {
414 // Angular distance is obtained from latitude in the new frame.
415 dist[i] = 90.0 - dist[i];
416
417 // Position angle is obtained from longitude in the new frame.
418 pa[i] = -pa[i];
419 if (pa[i] < -180.0) pa[i] += 360.0;
420 }
421
422 return 0;
423 }
424
425 //----------------------------------------------------------------------------
426
sphpad(int nfield,double lng0,double lat0,const double dist[],const double pa[],double lng[],double lat[])427 int sphpad(
428 int nfield,
429 double lng0,
430 double lat0,
431 const double dist[],
432 const double pa[],
433 double lng[],
434 double lat[])
435
436 {
437 int i;
438 double eul[5];
439
440 // Set the Euler angles for the coordinate transformation.
441 eul[0] = lng0;
442 eul[1] = 90.0 - lat0;
443 eul[2] = 0.0;
444 eul[3] = cosd(eul[1]);
445 eul[4] = sind(eul[1]);
446
447 for (i = 0; i < nfield; i++) {
448 // Latitude in the new frame is obtained from angular distance.
449 lat[i] = 90.0 - dist[i];
450
451 // Longitude in the new frame is obtained from position angle.
452 lng[i] = -pa[i];
453 }
454
455 // Transform field points to the old system.
456 sphx2s(eul, nfield, 0, 1, 1, lng, lat, lng, lat);
457
458 return 0;
459 }
460