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