← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to prim/general/libsrc/cel.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
prim/general/libsrc/cel.c
 
1
/*=============================================================================
 
2
*
 
3
*   WCSLIB - an implementation of the FITS WCS proposal.
 
4
*   Copyright (C) 1995, Mark Calabretta
 
5
*
 
6
*   This library is free software; you can redistribute it and/or modify it
 
7
*   under the terms of the GNU Library General Public License as published
 
8
*   by the Free Software Foundation; either version 2 of the License, or (at
 
9
*   your option) any later version.
 
10
*
 
11
*   This library is distributed in the hope that it will be useful, but
 
12
*   WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library
 
14
*   General Public License for more details.
 
15
*
 
16
*   You should have received a copy of the GNU Library General Public License
 
17
*   along with this library; if not, write to the Free Software Foundation,
 
18
*   Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
19
*
 
20
*   Correspondence concerning WCSLIB may be directed to:
 
21
*      Internet email: mcalabre@atnf.csiro.au
 
22
*      Postal address: Dr. Mark Calabretta,
 
23
*                      Australia Telescope National Facility,
 
24
*                      P.O. Box 76,
 
25
*                      Epping, NSW, 2121,
 
26
*                      AUSTRALIA
 
27
*
 
28
*.VERSION
 
29
* 070703        last modif
 
30
*=============================================================================
 
31
*
 
32
*   C routines which implement the FITS World Coordinate System (WCS)
 
33
*   convention.
 
34
*
 
35
*   Summary of routines
 
36
*   -------------------
 
37
*   These routines are provided as drivers for the lower level spherical
 
38
*   coordinate transformation and projection routines.  There are separate
 
39
*   driver routines for the forward, celfwd(), and reverse, celrev(),
 
40
*   transformations.
 
41
*
 
42
*   An initialization routine, celset(), computes intermediate values from
 
43
*   the transformation parameters but need not be called explicitly - see the
 
44
*   explanation of cel.flag below.
 
45
*
 
46
*
 
47
*   Initialization routine; celset()
 
48
*   --------------------------------
 
49
*   Initializes members of a celprm data structure which hold intermediate
 
50
*   values.  Note that this routine need not be called directly; it will be
 
51
*   invoked by celfwd() and celrev() if the "flag" structure member is
 
52
*   anything other than a predefined magic value.
 
53
*
 
54
*   Given:
 
55
*      pcode[4] char     WCS projection code (see below).
 
56
*
 
57
*   Given and returned:
 
58
*      cel      celprm*  Spherical coordinate transformation parameters
 
59
*                        (see below).
 
60
*      prj      prjprm*  Projection parameters (usage is described in the
 
61
*                        prologue to "proj.c").
 
62
*
 
63
*   Function return value:
 
64
*               int      Error status
 
65
*                           0: Success.
 
66
*                           1: Invalid coordinate transformation parameters.
 
67
*                           2: Ill-conditioned coordinate transformation
 
68
*                              parameters.
 
69
*
 
70
*   Forward transformation; celfwd()
 
71
*   --------------------------------
 
72
*   Compute (x,y) coordinates in the plane of projection from celestial
 
73
*   coordinates (lng,lat).
 
74
*
 
75
*   Given:
 
76
*      pcode[4] char     WCS projection code (see below).
 
77
*      lng,lat  double   Celestial longitude and latitude of the projected
 
78
*                        point, in degrees.
 
79
*
 
80
*   Given and returned:
 
81
*      cel      celprm*  Spherical coordinate transformation parameters
 
82
*                        (see below).
 
83
*
 
84
*   Returned:
 
85
*      phi,     double*  Longitude and latitude in the native coordinate
 
86
*      theta             system of the projection, in degrees.
 
87
*
 
88
*   Given and returned:
 
89
*      prj      prjprm*  Projection parameters (usage is described in the
 
90
*                        prologue to "proj.c").
 
91
*
 
92
*   Returned:
 
93
*      x,y      double*  Projected coordinates, "degrees".
 
94
*
 
95
*   Function return value:
 
96
*               int      Error status
 
97
*                           0: Success.
 
98
*                           1: Invalid coordinate transformation parameters.
 
99
*                           2: Invalid projection parameters.
 
100
*                           3: Invalid value of (lng,lat).
 
101
*
 
102
*   Reverse transformation; celrev()
 
103
*   --------------------------------
 
104
*   Compute the celestial coordinates (lng,lat) of the point with projected
 
105
*   coordinates (x,y).
 
106
*
 
107
*   Given:
 
108
*      pcode[4] char     WCS projection code (see below).
 
109
*      x,y      double   Projected coordinates, "degrees".
 
110
*
 
111
*   Given and returned:
 
112
*      prj      prjprm*  Projection parameters (usage is described in the
 
113
*                        prologue to "proj.c").
 
114
*
 
115
*   Returned:
 
116
*      phi,     double*  Longitude and latitude in the native coordinate
 
117
*      theta             system of the projection, in degrees.
 
118
*
 
119
*   Given and returned:
 
120
*      cel      celprm*  Spherical coordinate transformation parameters
 
121
*                        (see below).
 
122
*
 
123
*   Returned:
 
124
*      lng,lat  double*  Celestial longitude and latitude of the projected
 
125
*                        point, in degrees.
 
126
*
 
127
*   Function return value:
 
128
*               int      Error status
 
129
*                           0: Success.
 
130
*                           1: Invalid coordinate transformation parameters.
 
131
*                           2: Invalid projection parameters.
 
132
*                           3: Invalid value of (x,y).
 
133
*
 
134
*   Coordinate transformation parameters
 
135
*   ------------------------------------
 
136
*   The celprm struct consists of the following:
 
137
*
 
138
*      int flag
 
139
*         The celprm struct contains pointers to the forward and reverse
 
140
*         projection routines as well as intermediaries computed from the
 
141
*         reference coordinates (see below).  Whenever the projection code
 
142
*         (pcode) or any of ref[4] are set or changed then this flag must be
 
143
*         set to zero to signal the initialization routine, celset(), to
 
144
*         redetermine the function pointers and recompute intermediaries.
 
145
*         Once this has been done pcode itself is ignored.
 
146
*
 
147
*      double ref[4]
 
148
*         The first pair of values should be set to the celestial longitude
 
149
*         and latitude (usually right ascension and declination) of the
 
150
*         reference point of the projection.
 
151
*
 
152
*         The second pair of values are the native longitude and latitude of
 
153
*         the pole of the celestial coordinate system and correspond to the
 
154
*         FITS keywords LONGPOLE and LATPOLE.
 
155
*
 
156
*         LONGPOLE defaults to 0 degrees if the celestial latitude of the
 
157
*         reference point of the projection is greater than the native
 
158
*         latitude, otherwise 180 degrees.  (This is the condition for the
 
159
*         celestial latitude to increase in the same direction as the native
 
160
*         latitude at the reference point.)  ref[2] may be set to 999.0 to
 
161
*         indicate that the correct default should be substituted.
 
162
*
 
163
*         In some circumstances the latitude of the native pole may be
 
164
*         determined by the first three values only to within a sign and
 
165
*         LATPOLE is used to choose between the two solutions.  LATPOLE is
 
166
*         set in ref[3] and the solution closest to this value is used to
 
167
*         reset ref[3].  It is therefore legitimate, for example, to set
 
168
*         ref[3] to 999.0 to choose the more northerly solution - the default
 
169
*         if the LATPOLE card is omitted from the FITS header.  For the
 
170
*         special case where the reference point of the projection is at
 
171
*         native latitude zero, its celestial latitude is zero, and
 
172
*         LONGPOLE = +/- 90 then the native latitude of the pole is not
 
173
*         determined by the first three reference values and LATPOLE
 
174
*         specifies it completely.
 
175
*
 
176
*   The remaining members of the celprm struct are maintained by the
 
177
*   initialization routines and should not be modified.  This is done for the
 
178
*   sake of efficiency and to allow an arbitrary number of contexts to be
 
179
*   maintained simultaneously.
 
180
*
 
181
*      double euler[5]
 
182
*         Euler angles and associated intermediaries derived from the
 
183
*         coordinate reference values.
 
184
*      int (*prjfwd)()
 
185
*      int (*prjrev)()
 
186
*         Pointers to the forward and reverse projection routines.
 
187
*
 
188
*
 
189
*   WCS projection codes
 
190
*   --------------------
 
191
*   Zenithals/azimuthals:
 
192
*      AZP: zenithal/azimuthal perspective
 
193
*      TAN: gnomonic
 
194
*      SIN: synthesis (generalized orthographic)
 
195
*      STG: stereographic
 
196
*      ARC: zenithal/azimuthal equidistant
 
197
*      ZPN: zenithal/azimuthal polynomial
 
198
*      ZEA: zenithal/azimuthal equal area
 
199
*      AIR: Airy
 
200
*
 
201
*   Cylindricals:
 
202
*      CYP: cylindrical perspective
 
203
*      CAR: Cartesian
 
204
*      MER: Mercator
 
205
*      CEA: cylindrical equal area
 
206
*
 
207
*   Conics:
 
208
*      COP: conic perspective
 
209
*      COD: conic equidistant
 
210
*      COE: conic equal area
 
211
*      COO: conic orthomorphic
 
212
*
 
213
*   Polyconics:
 
214
*      BON: Bonne
 
215
*      PCO: polyconic
 
216
*
 
217
*   Pseudo-cylindricals:
 
218
*      GLS: Sanson-Flamsteed (global sinusoidal)
 
219
*      PAR: parabolic
 
220
*      MOL: Mollweide
 
221
*
 
222
*   Conventional:
 
223
*      AIT: Hammer-Aitoff
 
224
*
 
225
*   Quad-cubes:
 
226
*      CSC: COBE quadrilateralized spherical cube
 
227
*      QSC: quadrilateralized spherical cube
 
228
*      TSC: tangential spherical cube
 
229
*
 
230
*   Author: Mark Calabretta, Australia Telescope National Facility
 
231
*===========================================================================*/
 
232
 
 
233
#include <string.h>
 
234
 
 
235
#include <cel.h>
 
236
 
 
237
#ifndef __STDC__
 
238
#ifndef const
 
239
#define const
 
240
#endif
 
241
#endif
 
242
 
 
243
#ifdef SIGNBIT
 
244
#define signbit(X) ((X) < 0.0 ? 0 : 1)
 
245
#endif
 
246
 
 
247
int celset(pcode, cel, prj)
 
248
 
 
249
char pcode[4];
 
250
struct celprm *cel;
 
251
struct prjprm *prj;
 
252
 
 
253
{
 
254
   int dophip;
 
255
   const double tol = 1.0e-10;
 
256
   double clat0, cphip, cthe0, theta0, slat0, sphip, sthe0;
 
257
   double latp, latp1, latp2;
 
258
   double u, v, x, y, z;
 
259
 
 
260
   /* Set pointers to the forward and reverse projection routines. */
 
261
   if (strcmp(pcode, "AZP") == 0) {
 
262
      cel->prjfwd = azpfwd;
 
263
      cel->prjrev = azprev;
 
264
      theta0 = 90.0;
 
265
   } else if (strcmp(pcode, "TAN") == 0) {
 
266
      cel->prjfwd = tanfwd;
 
267
      cel->prjrev = tanrev;
 
268
      theta0 = 90.0;
 
269
   } else if (strcmp(pcode, "SIN") == 0) {
 
270
      cel->prjfwd = sinfwd;
 
271
      cel->prjrev = sinrev;
 
272
      theta0 = 90.0;
 
273
   } else if (strcmp(pcode, "STG") == 0) {
 
274
      cel->prjfwd = stgfwd;
 
275
      cel->prjrev = stgrev;
 
276
      theta0 = 90.0;
 
277
   } else if (strcmp(pcode, "ARC") == 0) {
 
278
      cel->prjfwd = arcfwd;
 
279
      cel->prjrev = arcrev;
 
280
      theta0 = 90.0;
 
281
   } else if (strcmp(pcode, "ZPN") == 0) {
 
282
      cel->prjfwd = zpnfwd;
 
283
      cel->prjrev = zpnrev;
 
284
      theta0 = 90.0;
 
285
   } else if (strcmp(pcode, "ZEA") == 0) {
 
286
      cel->prjfwd = zeafwd;
 
287
      cel->prjrev = zearev;
 
288
      theta0 = 90.0;
 
289
   } else if (strcmp(pcode, "AIR") == 0) {
 
290
      cel->prjfwd = airfwd;
 
291
      cel->prjrev = airrev;
 
292
      theta0 = 90.0;
 
293
   } else if (strcmp(pcode, "CYP") == 0) {
 
294
      cel->prjfwd = cypfwd;
 
295
      cel->prjrev = cyprev;
 
296
      theta0 = 0.0;
 
297
   } else if (strcmp(pcode, "CAR") == 0) {
 
298
      cel->prjfwd = carfwd;
 
299
      cel->prjrev = carrev;
 
300
      theta0 = 0.0;
 
301
   } else if (strcmp(pcode, "MER") == 0) {
 
302
      cel->prjfwd = merfwd;
 
303
      cel->prjrev = merrev;
 
304
      theta0 = 0.0;
 
305
   } else if (strcmp(pcode, "CEA") == 0) {
 
306
      cel->prjfwd = ceafwd;
 
307
      cel->prjrev = cearev;
 
308
      theta0 = 0.0;
 
309
   } else if (strcmp(pcode, "COP") == 0) {
 
310
      cel->prjfwd = copfwd;
 
311
      cel->prjrev = coprev;
 
312
      theta0 = prj->p[1];
 
313
   } else if (strcmp(pcode, "COD") == 0) {
 
314
      cel->prjfwd = codfwd;
 
315
      cel->prjrev = codrev;
 
316
      theta0 = prj->p[1];
 
317
   } else if (strcmp(pcode, "COE") == 0) {
 
318
      cel->prjfwd = coefwd;
 
319
      cel->prjrev = coerev;
 
320
      theta0 = prj->p[1];
 
321
   } else if (strcmp(pcode, "COO") == 0) {
 
322
      cel->prjfwd = coofwd;
 
323
      cel->prjrev = coorev;
 
324
      theta0 = prj->p[1];
 
325
   } else if (strcmp(pcode, "BON") == 0) {
 
326
      cel->prjfwd = bonfwd;
 
327
      cel->prjrev = bonrev;
 
328
      theta0 = 0.0;
 
329
   } else if (strcmp(pcode, "PCO") == 0) {
 
330
      cel->prjfwd = pcofwd;
 
331
      cel->prjrev = pcorev;
 
332
      theta0 = 0.0;
 
333
   } else if (strcmp(pcode, "GLS") == 0) {
 
334
      cel->prjfwd = glsfwd;
 
335
      cel->prjrev = glsrev;
 
336
      theta0 = 0.0;
 
337
   } else if (strcmp(pcode, "PAR") == 0) {
 
338
      cel->prjfwd = parfwd;
 
339
      cel->prjrev = parrev;
 
340
      theta0 = 0.0;
 
341
   } else if (strcmp(pcode, "AIT") == 0) {
 
342
      cel->prjfwd = aitfwd;
 
343
      cel->prjrev = aitrev;
 
344
      theta0 = 0.0;
 
345
   } else if (strcmp(pcode, "MOL") == 0) {
 
346
      cel->prjfwd = molfwd;
 
347
      cel->prjrev = molrev;
 
348
      theta0 = 0.0;
 
349
   } else if (strcmp(pcode, "CSC") == 0) {
 
350
      cel->prjfwd = cscfwd;
 
351
      cel->prjrev = cscrev;
 
352
      theta0 = 0.0;
 
353
   } else if (strcmp(pcode, "QSC") == 0) {
 
354
      cel->prjfwd = qscfwd;
 
355
      cel->prjrev = qscrev;
 
356
      theta0 = 0.0;
 
357
   } else if (strcmp(pcode, "TSC") == 0) {
 
358
      cel->prjfwd = tscfwd;
 
359
      cel->prjrev = tscrev;
 
360
      theta0 = 0.0;
 
361
   } else {
 
362
      /* Unrecognized projection code. */
 
363
      return 1;
 
364
   }
 
365
 
 
366
   /* Set default for native longitude of the celestial pole? */
 
367
   dophip = (cel->ref[2] == 999.0);
 
368
 
 
369
   /* Compute celestial coordinates of the native pole. */
 
370
   if (theta0 == 90.0) {
 
371
      /* Reference point is at the native pole. */
 
372
 
 
373
      if (dophip) {
 
374
         /* Set default for longitude of the celestial pole. */
 
375
         cel->ref[2] = 180.0;
 
376
      }
 
377
 
 
378
      latp = cel->ref[1];
 
379
      cel->ref[3] = latp;
 
380
 
 
381
      cel->euler[0] = cel->ref[0];
 
382
      cel->euler[1] = 90.0 - latp;
 
383
   } else {
 
384
      /* Reference point away from the native pole. */
 
385
 
 
386
      /* Set default for longitude of the celestial pole. */
 
387
      if (dophip) {
 
388
         cel->ref[2] = (cel->ref[1] < theta0) ? 180.0 : 0.0;
 
389
      }
 
390
 
 
391
      clat0 = wcs_cosd(cel->ref[1]);
 
392
      slat0 = wcs_sind(cel->ref[1]);
 
393
      cphip = wcs_cosd(cel->ref[2]);
 
394
      sphip = wcs_sind(cel->ref[2]);
 
395
      cthe0 = wcs_cosd(theta0);
 
396
      sthe0 = wcs_sind(theta0);
 
397
 
 
398
      x = cthe0*cphip;
 
399
      y = sthe0;
 
400
      z = sqrt(x*x + y*y);
 
401
      if (z == 0.0) {
 
402
         if (slat0 != 0.0) {
 
403
            return 1;
 
404
         }
 
405
 
 
406
         /* latp determined by LATPOLE in this case. */
 
407
         latp = cel->ref[3];
 
408
      } else {
 
409
         if (fabs(slat0/z) > 1.0) {
 
410
            return 1;
 
411
         }
 
412
 
 
413
         u = wcs_atan2d(y,x);
 
414
         v = wcs_acosd(slat0/z);
 
415
 
 
416
         latp1 = u + v;
 
417
         if (latp1 > 180.0) {
 
418
            latp1 -= 360.0;
 
419
         } else if (latp1 < -180.0) {
 
420
            latp1 += 360.0;
 
421
         }
 
422
 
 
423
         latp2 = u - v;
 
424
         if (latp2 > 180.0) {
 
425
            latp2 -= 360.0;
 
426
         } else if (latp2 < -180.0) {
 
427
            latp2 += 360.0;
 
428
         }
 
429
 
 
430
         if (fabs(cel->ref[3]-latp1) < fabs(cel->ref[3]-latp2)) {
 
431
            if (fabs(latp1) < 90.0+tol) {
 
432
               latp = latp1;
 
433
            } else {
 
434
               latp = latp2;
 
435
            }
 
436
         } else {
 
437
            if (fabs(latp2) < 90.0+tol) {
 
438
               latp = latp2;
 
439
            } else {
 
440
               latp = latp1;
 
441
            }
 
442
         }
 
443
 
 
444
         cel->ref[3] = latp;
 
445
      }
 
446
 
 
447
      cel->euler[1] = 90.0 - latp;
 
448
 
 
449
      z = wcs_cosd(latp)*clat0;
 
450
      if (fabs(z) < tol) {
 
451
         if (fabs(clat0) < tol) {
 
452
            /* Celestial pole at the reference point. */
 
453
            cel->euler[0] = cel->ref[0];
 
454
            cel->euler[1] = 90.0 - theta0;
 
455
         } else if (latp > 0.0) {
 
456
            /* Celestial pole at the native north pole.*/
 
457
            cel->euler[0] = cel->ref[0] + cel->ref[2] - 180.0;
 
458
            cel->euler[1] = 0.0;
 
459
         } else if (latp < 0.0) {
 
460
            /* Celestial pole at the native south pole. */
 
461
            cel->euler[0] = cel->ref[0] - cel->ref[2];
 
462
            cel->euler[1] = 180.0;
 
463
         }
 
464
      } else {
 
465
         x = (sthe0 - wcs_sind(latp)*slat0)/z;
 
466
         y =  sphip*cthe0/clat0;
 
467
         if (x == 0.0 && y == 0.0) {
 
468
            return 1;
 
469
         }
 
470
         cel->euler[0] = cel->ref[0] - wcs_atan2d(y,x);
 
471
      }
 
472
 
 
473
      /* Make euler[0] the same sign as ref[0]. */
 
474
      if (cel->ref[0] >= 0.0) {
 
475
         if (cel->euler[0] < 0.0) cel->euler[0] += 360.0;
 
476
      } else {
 
477
         if (cel->euler[0] > 0.0) cel->euler[0] -= 360.0;
 
478
      }
 
479
   }
 
480
 
 
481
   cel->euler[2] = cel->ref[2];
 
482
   cel->euler[3] = wcs_cosd(cel->euler[1]);
 
483
   cel->euler[4] = wcs_sind(cel->euler[1]);
 
484
   cel->flag = CELSET;
 
485
 
 
486
   /* Check for ill-conditioned parameters. */
 
487
   if (fabs(latp) > 90.0+tol) {
 
488
      return 2;
 
489
   }
 
490
 
 
491
   return 0;
 
492
}
 
493
 
 
494
/*--------------------------------------------------------------------------*/
 
495
 
 
496
int celfwd(pcode, lng, lat, cel, phi, theta, prj, x, y)
 
497
 
 
498
char   pcode[4];
 
499
double lng, lat;
 
500
struct celprm *cel;
 
501
double *phi, *theta;
 
502
struct prjprm *prj;
 
503
double *x, *y;
 
504
 
 
505
{
 
506
int    err;
 
507
 
 
508
int sphfwd();
 
509
 
 
510
 
 
511
 
 
512
   if (cel->flag != CELSET) {
 
513
      if (celset(pcode, cel, prj)) return 1;
 
514
   }
 
515
 
 
516
   /* Compute native coordinates. */
 
517
   sphfwd(lng, lat, cel->euler, phi, theta);
 
518
 
 
519
   /* Apply forward projection. */
 
520
   if ((err = cel->prjfwd(*phi, *theta, prj, x, y))) {
 
521
      return err == 1 ? 2 : 3;
 
522
   }
 
523
 
 
524
   return 0;
 
525
}
 
526
 
 
527
/*--------------------------------------------------------------------------*/
 
528
 
 
529
int celrev(pcode, x, y, prj, phi, theta, cel, lng, lat)
 
530
 
 
531
char   pcode[4];
 
532
double x, y;
 
533
struct prjprm *prj;
 
534
double *phi, *theta;
 
535
struct celprm *cel;
 
536
double *lng, *lat;
 
537
 
 
538
{
 
539
int    err;
 
540
 
 
541
int sphrev();
 
542
 
 
543
 
 
544
 
 
545
   if (cel->flag != CELSET) {
 
546
      if(celset(pcode, cel, prj)) return 1;
 
547
   }
 
548
 
 
549
   /* Apply reverse projection. */
 
550
   if ((err = cel->prjrev(x, y, prj, phi, theta))) {
 
551
      return err == 1 ? 2 : 3;
 
552
   }
 
553
 
 
554
   /* Compute native coordinates. */
 
555
   sphrev(*phi, *theta, cel->euler, lng, lat);
 
556
 
 
557
   return 0;
 
558
}