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: twcsfix.c,v 7.7 2021/07/12 06:36:49 mcalabre Exp $
23 *=============================================================================
24 *
25 * twcsfix tests the translation routines for non-standard WCS keyvalues, the
26 * wcsfix() suite. It also tests the change of celestial coordinate system
27 * routine, wcsccs(), and the spectral coordinate translation routine,
28 * wcssptr().
29 *
30 *---------------------------------------------------------------------------*/
31
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <string.h>
35
36 #include <wcs.h>
37 #include <wcserr.h>
38 #include <wcsfix.h>
39 #include <wcsprintf.h>
40 #include <wcsunits.h>
41 #include <wcsutil.h>
42
43
44 void parser(struct wcsprm *);
45
46 const int NAXIS = 3;
47 const double CRPIX[3] = {90.0, 90.0, 1.0};
48 const double PC[3][3] = {{ 1.0, 0.0, 0.0},
49 { 0.0, 1.0, 0.0},
50 { 0.0, 0.0, 1.0}};
51 const double CDELT[3] = {-1.0, 1.0, 19.68717093222};
52
53 char CUNIT[3][9] = {"ARCSEC", "ARCSEC", "KM/SEC"};
54
55 // N.B. non-standard spectral axis.
56 char CTYPE[3][9] = {"RA---NCP", "DEC--NCP", "FELO-HEL"};
57
58 // B1950.0 equatorial coordinates of the galactic pole.
59 // CUNITia is set to ARCSEC as an additional test.
60 const double CRVAL[3] = {192.2500*3600.0, 27.4000*3600.0, 5569.27104};
61 const double RESTFRQ = 1.42040575e9;
62 const double RESTWAV = 0.0;
63
64 // N.B. non-standard date-time format.
65 const char DATEOBS[] = "1957/02/15 01:10:00";
66 const char DATEBEG[] = "1957/02/15 01:10:00";
67 const char DATEAVG[] = "1957/02/15 02:10:00";
68 const char DATEEND[] = "1957/02/15 03:10:00";
69
70 const double BEPOCH = 1957.124382563;
71 const double MJDBEG = 35884.048611;
72
73 const double OBSGEO_L = 148.263510;
74 const double OBSGEO_B = -32.998406;
75 const double OBSGEO_H = 411.793;
76
77 const double EQUINOX = 1950.0;
78
79 // For testing spcfix().
80 const int VELREF = 2;
81 const char SPECSYS[] = "BARYCENT";
82
main()83 int main()
84
85 {
86 int stat[NWCSFIX];
87
88 wcsprintf("Testing WCSLIB translator for non-standard usage (twcsfix.c)\n"
89 "------------------------------------------------------------\n\n");
90
91 struct wcsprm wcs0;
92 wcs0.flag = -1;
93 parser(&wcs0);
94
95 // Note: to print the unfixed wcsprm struct using wcsprt() the struct
96 // would first have to be initialized by wcsset(). However, if the struct
97 // contains non-standard keyvalues then wcsset() will either fix them
98 // itself or else fail (e.g. for non-standard units). Thus, in general,
99 // wcsprt() cannot be used to print the unmodified struct.
100
101 // Fix non-standard WCS keyvalues.
102 struct wcserr info[NWCSFIX];
103 wcserr_enable(1);
104 int status = wcsfixi(7, 0, &wcs0, stat, info);
105 wcsprintf("wcsfix status returns: (");
106 for (int i = 0; i < NWCSFIX; i++) {
107 wcsprintf(i ? ", %d" : "%d", stat[i]);
108 }
109 wcsprintf(")\n");
110
111 for (int i = 0; i < NWCSFIX; i++) {
112 if (info[i].status < -1 || 0 < info[i].status) {
113 wcsprintf("\n");
114 wcserr_prt(info+i, 0x0);
115
116 // Free memory used to store the message.
117 if (info[i].msg) wcsdealloc(info[i].msg);
118 }
119 }
120
121 if (status) {
122 wcsprintf("\nwcsfix error %d", status);
123 return 1;
124 }
125
126 // Extract information from the FITS header.
127 if (wcsset(&wcs0)) {
128 wcsprintf("\n");
129 wcsperr(&wcs0, 0x0);
130 return 1;
131 }
132
133 wcsprintf("\n");
134 wcsprt(&wcs0);
135 wcsprintf("\n------------------------------------"
136 "------------------------------------\n");
137
138 // Make a copy of the wcsprm struct.
139 struct wcsprm wcs1;
140 wcs1.flag = -1;
141 if (wcssub(1, &wcs0, 0x0, 0x0, &wcs1)) {
142 wcsperr(&wcs1, 0x0);
143 return 1;
144 }
145
146 // Transform equatorial B1950 to galactic coordinates. The WCS has been
147 // constructed with the galactic pole coincident with the native pole of
148 // the projection in order to test the resolution of an indeterminacy.
149 if (wcsccs(&wcs1, 123.0, 27.4, 192.25, "GLON", "GLAT", 0x0, 0.0, "G")) {
150 wcsperr(&wcs1, 0x0);
151 return 1;
152 }
153
154 // Should now have a 'VOPT-F2W' axis, translate it to frequency.
155 char ctypeS[9];
156 strcpy(ctypeS, "FREQ-???");
157 int i = -1;
158 if (wcssptr(&wcs1, &i, ctypeS)) {
159 wcsperr(&wcs1, 0x0);
160 return 1;
161 }
162
163 if (wcsset(&wcs1)) {
164 wcsperr(&wcs1, 0x0);
165 return 1;
166 }
167
168 wcsprt(&wcs1);
169
170 // Print before-and-afters.
171 printf("\nOriginal and new coordinates of reference point "
172 "(%4.1f, %4.1f, %3.1f), lonpole, and latpole:\n",
173 CRPIX[0], CRPIX[1], CRPIX[2]);
174 printf("%14.6f, %14.6f, %14.2f, %14.6f, %14.6f\n",
175 wcs0.crval[0], wcs0.crval[1], wcs0.crval[2], wcs0.lonpole, wcs0.latpole);
176 printf("%14.6f, %14.6f, %14.2f, %14.6f, %14.6f\n",
177 wcs1.crval[0], wcs1.crval[1], wcs1.crval[2], wcs1.lonpole, wcs1.latpole);
178
179 // Compute B1950 coordinates of a field point.
180 double pixcrd[3] = {1000.0, 1000.0, 1.0};
181 printf("\nOriginal and new coordinates of field point "
182 "(%5.1f, %5.1f, %3.1f):\n", pixcrd[0], pixcrd[1], pixcrd[2]);
183
184 double imgcrd[3], phi, theta, world[3];
185 if (wcsp2s(&wcs0, 1, 3, pixcrd, imgcrd, &phi, &theta, world, stat)) {
186 wcsperr(&wcs0, 0x0);
187 return 1;
188 }
189
190 printf("%14.6f, %14.6f, %14.2f\n", world[0], world[1], world[2]);
191
192 // Compute galactic coordinates of the same field point.
193 if (wcsp2s(&wcs1, 1, 3, pixcrd, imgcrd, &phi, &theta, world, stat)) {
194 wcsperr(&wcs1, 0x0);
195 return 1;
196 }
197
198 printf("%14.6f, %14.6f, %14.2f\n", world[0], world[1], world[2]);
199
200 wcsfree(&wcs0);
201 wcsfree(&wcs1);
202
203 return 0;
204 }
205
206 //----------------------------------------------------------------------------
207
parser(wcs)208 void parser(wcs)
209
210 struct wcsprm *wcs;
211
212 {
213 int i, j;
214 double *pcij;
215
216 // In practice a parser would read the FITS header until it encountered
217 // the NAXIS keyword which must occur near the start, before any of the
218 // WCS keywords. It would then use wcsini() to allocate memory for
219 // arrays in the wcsprm struct and set default values. In this
220 // simulation the header keyvalues are set as global variables.
221 wcsnpv(2);
222 wcsini(1, NAXIS, wcs);
223
224
225 // Now the parser scans the FITS header, identifying WCS keywords and
226 // loading their values into the appropriate elements of the wcsprm
227 // struct.
228
229 for (j = 0; j < NAXIS; j++) {
230 wcs->crpix[j] = CRPIX[j];
231 }
232
233 pcij = wcs->pc;
234 for (i = 0; i < NAXIS; i++) {
235 for (j = 0; j < NAXIS; j++) {
236 *(pcij++) = PC[i][j];
237 }
238 }
239
240 for (i = 0; i < NAXIS; i++) {
241 wcs->cdelt[i] = CDELT[i];
242 }
243
244 for (i = 0; i < NAXIS; i++) {
245 strcpy(wcs->cunit[i], &CUNIT[i][0]);
246 }
247
248 for (i = 0; i < NAXIS; i++) {
249 strcpy(wcs->ctype[i], &CTYPE[i][0]);
250 }
251
252 for (i = 0; i < NAXIS; i++) {
253 wcs->crval[i] = CRVAL[i];
254 }
255
256 wcs->restfrq = RESTFRQ;
257 wcs->restwav = RESTWAV;
258
259 wcs->pv[0].i = -1;
260 wcs->pv[0].m = -1;
261 wcs->pv[0].value = -1.0;
262 wcs->npv = 1;
263
264 wcs->velref = VELREF;
265
266 // dateobs and datebeg will be set by datfix().
267 strcpy(wcs->dateavg, DATEAVG);
268 strcpy(wcs->dateend, DATEEND);
269 wcs->bepoch = BEPOCH;
270 wcs->mjdbeg = MJDBEG;
271
272 wcs->obsgeo[3] = OBSGEO_L;
273 wcs->obsgeo[4] = OBSGEO_B;
274 wcs->obsgeo[5] = OBSGEO_H;
275
276 wcs->equinox = EQUINOX;
277
278 strcpy(wcs->specsys, SPECSYS);
279
280 return;
281 }
282