diff options
Diffstat (limited to 'funtools/wcs/wcslib.c')
-rw-r--r-- | funtools/wcs/wcslib.c | 1334 |
1 files changed, 1334 insertions, 0 deletions
diff --git a/funtools/wcs/wcslib.c b/funtools/wcs/wcslib.c new file mode 100644 index 0000000..31313e1 --- /dev/null +++ b/funtools/wcs/wcslib.c @@ -0,0 +1,1334 @@ +/*============================================================================= +* +* WCSLIB - an implementation of the FITS WCS proposal. +* Copyright (C) 1995-2002, Mark Calabretta +* +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2 of the License, or (at your option) any later version. +* +* This library is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +* Lesser General Public License for more details. +* +* You should have received a copy of the GNU Lesser General Public +* License along with this library; if not, write to the Free Software +* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +* +* Correspondence concerning WCSLIB may be directed to: +* Internet email: mcalabre@atnf.csiro.au +* Postal address: Dr. Mark Calabretta, +* Australia Telescope National Facility, +* P.O. Box 76, +* Epping, NSW, 2121, +* AUSTRALIA +* +*============================================================================= +* +* C routines which implement the FITS World Coordinate System (WCS) +* convention. +* +* Summary of routines +* ------------------- +* wcsfwd() and wcsrev() are high level driver routines for the WCS linear +* transformation, spherical coordinate transformation, and spherical +* projection routines. +* +* Given either the celestial longitude or latitude plus an element of the +* pixel coordinate a hybrid routine, wcsmix(), iteratively solves for the +* unknown elements. +* +* An initialization routine, wcsset(), computes indices from the ctype +* array but need not be called explicitly - see the explanation of +* wcs.flag below. +* +* +* Initialization routine; wcsset() +* -------------------------------- +* Initializes elements of a wcsprm data structure which holds indices into +* the coordinate arrays. Note that this routine need not be called directly; +* it will be invoked by wcsfwd() and wcsrev() if the "flag" structure member +* is anything other than a predefined magic value. +* +* Given: +* naxis const int +* Number of image axes. +* ctype[][9] +* const char +* Coordinate axis types corresponding to the FITS +* CTYPEn header cards. +* +* Returned: +* wcs wcsprm* Indices for the celestial coordinates obtained +* by parsing the ctype[] array (see below). +* +* Function return value: +* int Error status +* 0: Success. +* 1: Inconsistent or unrecognized coordinate axis +* types. +* +* +* Forward transformation; wcsfwd() +* -------------------------------- +* Compute the pixel coordinate for given world coordinates. +* +* Given: +* ctype[][9] +* const char +* Coordinate axis types corresponding to the FITS +* CTYPEn header cards. +* +* Given or returned: +* wcs wcsprm* Indices for the celestial coordinates obtained +* by parsing the ctype[] array (see below). +* +* Given: +* world const double[] +* World coordinates. world[wcs->lng] and +* world[wcs->lat] are the celestial longitude and +* latitude, in degrees. +* +* Given: +* crval const double[] +* Coordinate reference values corresponding to the FITS +* CRVALn header cards (see note 2). +* +* Given and returned: +* cel celprm* Spherical coordinate transformation parameters (usage +* is described in the prologue to "cel.c"). +* +* Returned: +* phi, double* Longitude and latitude in the native coordinate +* theta system of the projection, in degrees. +* +* Given and returned: +* prj prjprm* Projection parameters (usage is described in the +* prologue to "proj.c"). +* +* Returned: +* imgcrd double[] Image coordinate. imgcrd[wcs->lng] and +* imgcrd[wcs->lat] are the projected x-, and +* y-coordinates, in "degrees". For quadcube +* projections with a CUBEFACE axis the face number is +* also returned in imgcrd[wcs->cubeface]. +* +* Given and returned: +* lin linprm* Linear transformation parameters (usage is described +* in the prologue to "lin.c"). +* +* Returned: +* pixcrd double[] Pixel coordinate. +* +* Function return value: +* int Error status +* 0: Success. +* 1: Invalid coordinate transformation parameters. +* 2: Invalid projection parameters. +* 3: Invalid world coordinate. +* 4: Invalid linear transformation parameters. +* +* +* Reverse transformation; wcsrev() +* -------------------------------- +* Compute world coordinates for a given pixel coordinate. +* +* Given: +* ctype[][9] +* const char +* Coordinate axis types corresponding to the FITS +* CTYPEn header cards. +* +* Given or returned: +* wcs wcsprm* Indices for the celestial coordinates obtained +* by parsing the ctype[] array (see below). +* +* Given: +* pixcrd const double[] +* Pixel coordinate. +* +* Given and returned: +* lin linprm* Linear transformation parameters (usage is described +* in the prologue to "lin.c"). +* +* Returned: +* imgcrd double[] Image coordinate. imgcrd[wcs->lng] and +* imgcrd[wcs->lat] are the projected x-, and +* y-coordinates, in "degrees". +* +* Given and returned: +* prj prjprm* Projection parameters (usage is described in the +* prologue to "proj.c"). +* +* Returned: +* phi, double* Longitude and latitude in the native coordinate +* theta system of the projection, in degrees. +* +* Given: +* crval const double[] +* Coordinate reference values corresponding to the FITS +* CRVALn header cards (see note 2). +* +* Given and returned: +* cel celprm* Spherical coordinate transformation parameters +* (usage is described in the prologue to "cel.c"). +* +* Returned: +* world double[] World coordinates. world[wcs->lng] and +* world[wcs->lat] are the celestial longitude and +* latitude, in degrees. +* +* Function return value: +* int Error status +* 0: Success. +* 1: Invalid coordinate transformation parameters. +* 2: Invalid projection parameters. +* 3: Invalid pixel coordinate. +* 4: Invalid linear transformation parameters. +* +* +* Hybrid transformation; wcsmix() +* ------------------------------- +* Given either the celestial longitude or latitude plus an element of the +* pixel coordinate solve for the remaining elements by iterating on the +* unknown celestial coordinate element using wcsfwd(). +* +* Given: +* ctype[][9] +* const char +* Coordinate axis types corresponding to the FITS +* CTYPEn header cards. +* +* Given or returned: +* wcs wcsprm* Indices for the celestial coordinates obtained +* by parsing the ctype[] array (see below). +* +* Given: +* mixpix const int +* Which element of the pixel coordinate is given. +* mixcel const int +* Which element of the celestial coordinate is +* given: +* 1: Celestial longitude is given in +* world[wcs->lng], latitude returned in +* world[wcs->lat]. +* 2: Celestial latitude is given in +* world[wcs->lat], longitude returned in +* world[wcs->lng]. +* vspan[2] const double +* Solution interval for the celestial coordinate, in +* degrees. The ordering of the two limits is +* irrelevant. Longitude ranges may be specified with +* any convenient normalization, for example [-120,+120] +* is the same as [240,480], except that the solution +* will be returned with the same normalization, i.e. +* lie within the interval specified. +* vstep const double +* Step size for solution search, in degrees. If zero, +* a sensible, although perhaps non-optimal default will +* be used. +* viter int +* If a solution is not found then the step size will be +* halved and the search recommenced. viter controls +* how many times the step size is halved. The allowed +* range is 5 - 10. +* +* Given and returned: +* world double[] World coordinates. world[wcs->lng] and +* world[wcs->lat] are the celestial longitude and +* latitude, in degrees. Which is given and which +* returned depends on the value of mixcel. All other +* elements are given. +* +* Given: +* crval const double[] +* Coordinate reference values corresponding to the FITS +* CRVALn header cards (see note 2). +* +* Given and returned: +* cel celprm* Spherical coordinate transformation parameters +* (usage is described in the prologue to "cel.c"). +* +* Returned: +* phi, double* Longitude and latitude in the native coordinate +* theta system of the projection, in degrees. +* +* Given and returned: +* prj prjprm* Projection parameters (usage is described in the +* prologue to "proj.c"). +* +* Returned: +* imgcrd double[] Image coordinate. imgcrd[wcs->lng] and +* imgcrd[wcs->lat] are the projected x-, and +* y-coordinates, in "degrees". +* +* Given and returned: +* lin linprm* Linear transformation parameters (usage is described +* in the prologue to "lin.c"). +* +* Given and returned: +* pixcrd double[] Pixel coordinate. The element indicated by mixpix is +* given and the remaining elements are returned. +* +* Function return value: +* int Error status +* 0: Success. +* 1: Invalid coordinate transformation parameters. +* 2: Invalid projection parameters. +* 3: Coordinate transformation error. +* 4: Invalid linear transformation parameters. +* 5: No solution found in the specified interval. +* +* +* Notes +* ----- +* 1) The CTYPEn must in be upper case and there must be 0 or 1 pair of +* matched celestial axis types. The ctype[][9] should be padded with +* blanks on the right and null-terminated. +* +* 2) Elements of the crval[] array which correspond to celestial axes are +* ignored, the reference coordinate values in cel->ref[0] and +* cel->ref[1] are the ones used. +* +* 3) These functions recognize the NCP projection and convert it to the +* equivalent SIN projection. +* +* They also recognize GLS as a synonym for SFL. +* +* 4) The quadcube projections (TSC, CSC, QSC) may be represented in FITS in +* either of two ways: +* +* a) The six faces may be laid out in one plane and numbered as +* follows: +* +* 0 +* +* 4 3 2 1 4 3 2 +* +* 5 +* +* Faces 2, 3 and 4 may appear on one side or the other (or both). +* The forward routines map faces 2, 3 and 4 to the left but the +* inverse routines accept them on either side. +* +* b) The "COBE" convention in which the six faces are stored in a +* three-dimensional structure using a "CUBEFACE" axis indexed from +* 0 to 5 as above. +* +* These routines support both methods; wcsset() determines which is +* being used by the presence or absence of a CUBEFACE axis in ctype[]. +* wcsfwd() and wcsrev() translate the CUBEFACE axis representation to +* the single plane representation understood by the lower-level WCSLIB +* projection routines. +* +* +* WCS indexing parameters +* ----------------------- +* The wcsprm struct consists of the following: +* +* int flag +* The wcsprm struct contains indexes and other information derived +* from the CTYPEn. Whenever any of the ctype[] are set or changed +* this flag must be set to zero to signal the initialization routine, +* wcsset() to redetermine the indices. The flag is set to 999 if +* there is no celestial axis pair in the CTYPEn. +* +* char pcode[4] +* The WCS projection code. +* +* char lngtyp[5], lattyp[5] +* WCS celestial axis types. +* +* int lng,lat +* Indices into the imgcrd[], and world[] arrays as described above. +* These may also serve as indices for the celestial longitude and +* latitude axes in the pixcrd[] array provided that the PC matrix +* does not transpose axes. +* +* int cubeface +* Index into the pixcrd[] array for the CUBEFACE axis. This is +* optionally used for the quadcube projections where each cube face is +* stored on a separate axis. +* +* +* wcsmix() algorithm +* ------------------ +* Initially the specified solution interval is checked to see if it's a +* "crossing" interval. If it isn't, a search is made for a crossing +* solution by iterating on the unknown celestial coordinate starting at +* the upper limit of the solution interval and decrementing by the +* specified step size. A crossing is indicated if the trial value of the +* pixel coordinate steps through the value specified. If a crossing +* interval is found then the solution is determined by a modified form of +* "regula falsi" division of the crossing interval. If no crossing +* interval was found within the specified solution interval then a search +* is made for a "non-crossing" solution as may arise from a point of +* tangency. The process is complicated by having to make allowance for +* the discontinuities that occur in all map projections. +* +* Once one solution has been determined others may be found by subsequent +* invokations of wcsmix() with suitably restricted solution intervals. +* +* Note the circumstance which arises when the solution point lies at a +* native pole of a projection in which the pole is represented as a +* finite curve, for example the zenithals and conics. In such cases two +* or more valid solutions may exist but WCSMIX only ever returns one. +* +* Because of its generality wcsmix() is very compute-intensive. For +* compute-limited applications more efficient special-case solvers could +* be written for simple projections, for example non-oblique cylindrical +* projections. +* +* Author: Mark Calabretta, Australia Telescope National Facility +* $Id: wcs.c,v 2.23 2002/04/03 01:25:29 mcalabre Exp $ +*===========================================================================*/ + +#include <stdio.h> +#include <math.h> +#include <string.h> +#include <stdio.h> +#include "wcslib.h" + +/* Map error number to error message for each function. */ +const char *wcsset_errmsg[] = { + 0, + "Inconsistent or unrecognized coordinate axis types"}; + +const char *wcsfwd_errmsg[] = { + 0, + "Invalid coordinate transformation parameters", + "Invalid projection parameters", + "Invalid world coordinate", + "Invalid linear transformation parameters"}; + +const char *wcsrev_errmsg[] = { + 0, + "Invalid coordinate transformation parameters", + "Invalid projection parameters", + "Invalid pixel coordinate", + "Invalid linear transformation parameters"}; + +const char *wcsmix_errmsg[] = { + 0, + "Invalid coordinate transformation parameters", + "Invalid projection parameters", + "Coordinate transformation error", + "Invalid linear transformation parameters", + "No solution found in the specified interval"}; + +#define signb(X) ((X) < 0.0 ? 1 : 0) + +int +wcsset (naxis, ctype, wcs) + +const int naxis; +const char ctype[][9]; +struct wcsprm *wcs; + +{ + int nalias = 2; + char aliases [2][4] = {"NCP", "GLS"}; + + int j, k; + int *ndx = NULL; + char requir[9]; + + strcpy(wcs->pcode, ""); + strcpy(requir, ""); + wcs->lng = -1; + wcs->lat = -1; + wcs->cubeface = -1; + + for (j = 0; j < naxis; j++) { + if (ctype[j][4] != '-') { + if (strcmp(ctype[j], "CUBEFACE") == 0) { + if (wcs->cubeface == -1) { + wcs->cubeface = j; + } else { + /* Multiple CUBEFACE axes! */ + return 1; + } + } + continue; + } + + /* Got an axis qualifier, is it a recognized WCS projection? */ + for (k = 0; k < npcode; k++) { + if (strncmp(&ctype[j][5], pcodes[k], 3) == 0) break; + } + + if (k == npcode) { + /* Maybe it's a projection alias. */ + for (k = 0; k < nalias; k++) { + if (strncmp(&ctype[j][5], aliases[k], 3) == 0) break; + } + + /* Not recognized. */ + if (k == nalias) { + continue; + } + } + + /* Parse the celestial axis type. */ + if (strcmp(wcs->pcode, "") == 0) { + sprintf(wcs->pcode, "%.3s", &ctype[j][5]); + + if (strncmp(ctype[j], "RA--", 4) == 0) { + wcs->lng = j; + strcpy(wcs->lngtyp, "RA"); + strcpy(wcs->lattyp, "DEC"); + ndx = &wcs->lat; + sprintf(requir, "DEC--%s", wcs->pcode); + } else if (strncmp(ctype[j], "DEC-", 4) == 0) { + wcs->lat = j; + strcpy(wcs->lngtyp, "RA"); + strcpy(wcs->lattyp, "DEC"); + ndx = &wcs->lng; + sprintf(requir, "RA---%s", wcs->pcode); + } else if (strncmp(&ctype[j][1], "LON", 3) == 0) { + wcs->lng = j; + sprintf(wcs->lngtyp, "%cLON", ctype[j][0]); + sprintf(wcs->lattyp, "%cLAT", ctype[j][0]); + ndx = &wcs->lat; + sprintf(requir, "%s-%s", wcs->lattyp, wcs->pcode); + } else if (strncmp(&ctype[j][1], "LAT", 3) == 0) { + wcs->lat = j; + sprintf(wcs->lngtyp, "%cLON", ctype[j][0]); + sprintf(wcs->lattyp, "%cLAT", ctype[j][0]); + ndx = &wcs->lng; + sprintf(requir, "%s-%s", wcs->lngtyp, wcs->pcode); + } else if (strncmp(&ctype[j][2], "LN", 2) == 0) { + wcs->lng = j; + sprintf(wcs->lngtyp, "%c%cLN", ctype[j][0], ctype[j][1]); + sprintf(wcs->lattyp, "%c%cLT", ctype[j][0], ctype[j][1]); + ndx = &wcs->lat; + sprintf(requir, "%s-%s", wcs->lattyp, wcs->pcode); + } else if (strncmp(&ctype[j][2], "LT", 2) == 0) { + wcs->lat = j; + sprintf(wcs->lngtyp, "%c%cLN", ctype[j][0], ctype[j][1]); + sprintf(wcs->lattyp, "%c%cLT", ctype[j][0], ctype[j][1]); + ndx = &wcs->lng; + sprintf(requir, "%s-%s", wcs->lngtyp, wcs->pcode); + } else { + /* Unrecognized celestial type. */ + return 1; + } + } else { + if (strncmp(ctype[j], requir, 8) != 0) { + /* Inconsistent projection types. */ + return 1; + } + + if (ndx == NULL) + return 1; + *ndx = j; + strcpy(requir, ""); + } + } + + if (strcmp(requir, "")) { + /* Unmatched celestial axis. */ + return 1; + } + + /* Do simple alias translations. */ + if (strncmp(wcs->pcode, "GLS", 3) == 0) { + strcpy(wcs->pcode, "SFL"); + } + + if (strcmp(wcs->pcode, "")) { + wcs->flag = WCSSET; + } else { + /* Signal for no celestial axis pair. */ + wcs->flag = 999; + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int +wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj, imgcrd, lin, pixcrd) + +const char ctype[][9]; +struct wcsprm* wcs; +const double world[]; +const double crval[]; +struct celprm *cel; +double *phi, *theta; +struct prjprm *prj; +double imgcrd[]; +struct linprm *lin; +double pixcrd[]; + +{ + int err, j; + double offset; + + /* Initialize if required. */ + if (wcs->flag != WCSSET) { + if (wcsset(lin->naxis, ctype, wcs)) return 1; + } + + /* Convert to relative physical coordinates. */ + for (j = 0; j < lin->naxis; j++) { + if (j == wcs->lng) continue; + if (j == wcs->lat) continue; + imgcrd[j] = world[j] - crval[j]; + } + + if (wcs->flag != 999) { + /* Compute projected coordinates. */ + if (strcmp(wcs->pcode, "NCP") == 0) { + /* Convert NCP to SIN. */ + if (cel->ref[1] == 0.0) { + return 2; + } + + strcpy(wcs->pcode, "SIN"); + prj->p[1] = 0.0; + prj->p[2] = cosdeg (cel->ref[1])/sindeg (cel->ref[1]); + prj->flag = (prj->flag < 0) ? -1 : 0; + } + + if ((err = celfwd(wcs->pcode, world[wcs->lng], world[wcs->lat], cel, + phi, theta, prj, &imgcrd[wcs->lng], &imgcrd[wcs->lat]))) { + return err; + } + + /* Do we have a CUBEFACE axis? */ + if (wcs->cubeface != -1) { + /* Separation between faces. */ + if (prj->r0 == 0.0) { + offset = 90.0; + } else { + offset = prj->r0*PI/2.0; + } + + /* Stack faces in a cube. */ + if (imgcrd[wcs->lat] < -0.5*offset) { + imgcrd[wcs->lat] += offset; + imgcrd[wcs->cubeface] = 5.0; + } else if (imgcrd[wcs->lat] > 0.5*offset) { + imgcrd[wcs->lat] -= offset; + imgcrd[wcs->cubeface] = 0.0; + } else if (imgcrd[wcs->lng] > 2.5*offset) { + imgcrd[wcs->lng] -= 3.0*offset; + imgcrd[wcs->cubeface] = 4.0; + } else if (imgcrd[wcs->lng] > 1.5*offset) { + imgcrd[wcs->lng] -= 2.0*offset; + imgcrd[wcs->cubeface] = 3.0; + } else if (imgcrd[wcs->lng] > 0.5*offset) { + imgcrd[wcs->lng] -= offset; + imgcrd[wcs->cubeface] = 2.0; + } else { + imgcrd[wcs->cubeface] = 1.0; + } + } + } + + /* Apply forward linear transformation. */ + if (linfwd(imgcrd, lin, pixcrd)) { + return 4; + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int +wcsrev(ctype, wcs, pixcrd, lin, imgcrd, prj, phi, theta, crval, cel, world) + +const char ctype[][9]; +struct wcsprm *wcs; +const double pixcrd[]; +struct linprm *lin; +double imgcrd[]; +struct prjprm *prj; +double *phi, *theta; +const double crval[]; +struct celprm *cel; +double world[]; + +{ + int err, face, j; + double offset; + + /* Initialize if required. */ + if (wcs->flag != WCSSET) { + if (wcsset(lin->naxis, ctype, wcs)) return 1; + } + + /* Apply reverse linear transformation. */ + if (linrev(pixcrd, lin, imgcrd)) { + return 4; + } + + /* Convert to world coordinates. */ + for (j = 0; j < lin->naxis; j++) { + if (j == wcs->lng) continue; + if (j == wcs->lat) continue; + world[j] = imgcrd[j] + crval[j]; + } + + + if (wcs->flag != 999) { + /* Do we have a CUBEFACE axis? */ + if (wcs->cubeface != -1) { + face = (int)(imgcrd[wcs->cubeface] + 0.5); + if (fabs(imgcrd[wcs->cubeface]-face) > 1e-10) { + return 3; + } + + /* Separation between faces. */ + if (prj->r0 == 0.0) { + offset = 90.0; + } else { + offset = prj->r0*PI/2.0; + } + + /* Lay out faces in a plane. */ + switch (face) { + case 0: + imgcrd[wcs->lat] += offset; + break; + case 1: + break; + case 2: + imgcrd[wcs->lng] += offset; + break; + case 3: + imgcrd[wcs->lng] += offset*2; + break; + case 4: + imgcrd[wcs->lng] += offset*3; + break; + case 5: + imgcrd[wcs->lat] -= offset; + break; + default: + return 3; + } + } + + /* Compute celestial coordinates. */ + if (strcmp(wcs->pcode, "NCP") == 0) { + /* Convert NCP to SIN. */ + if (cel->ref[1] == 0.0) { + return 2; + } + + strcpy(wcs->pcode, "SIN"); + prj->p[1] = 0.0; + prj->p[2] = cosdeg (cel->ref[1])/sindeg (cel->ref[1]); + prj->flag = (prj->flag < 0) ? -1 : 0; + } + + if ((err = celrev(wcs->pcode, imgcrd[wcs->lng], imgcrd[wcs->lat], prj, + phi, theta, cel, &world[wcs->lng], &world[wcs->lat]))) { + return err; + } + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int +wcsmix(ctype, wcs, mixpix, mixcel, vspan, vstep, viter, world, crval, cel, + phi, theta, prj, imgcrd, lin, pixcrd) + +const char ctype[][9]; +struct wcsprm *wcs; +const int mixpix, mixcel; +const double vspan[2], vstep; +int viter; +double world[]; +const double crval[]; +struct celprm *cel; +double *phi, *theta; +struct prjprm *prj; +double imgcrd[]; +struct linprm *lin; +double pixcrd[]; + +{ + const int niter = 60; + int crossed, err, istep, iter, j, k, nstep, retry; + const double tol = 1.0e-10; + const double tol2 = 100.0*tol; + double lambda, span[2], step; + double pixmix; + double dlng, lng, lng0, lng0m, lng1, lng1m; + double dlat, lat, lat0, lat0m, lat1, lat1m; + double d, d0, d0m, d1, d1m; + double dx = 0.0; + double dabs, dmin, lmin; + double dphi, phi0, phi1; + struct celprm cel0; + + /* Initialize if required. */ + if (wcs->flag != WCSSET) { + if (wcsset(lin->naxis, ctype, wcs)) return 1; + } + + /* Check vspan. */ + if (vspan[0] <= vspan[1]) { + span[0] = vspan[0]; + span[1] = vspan[1]; + } else { + /* Swap them. */ + span[0] = vspan[1]; + span[1] = vspan[0]; + } + + /* Check vstep. */ + step = fabs(vstep); + if (step == 0.0) { + step = (span[1] - span[0])/10.0; + if (step > 1.0 || step == 0.0) step = 1.0; + } + + /* Check viter. */ + nstep = viter; + if (nstep < 5) { + nstep = 5; + } else if (nstep > 10) { + nstep = 10; + } + + /* Given pixel element. */ + pixmix = pixcrd[mixpix]; + + /* Iterate on the step size. */ + for (istep = 0; istep <= nstep; istep++) { + if (istep) step /= 2.0; + + /* Iterate on the sky coordinate between the specified range. */ + if (mixcel == 1) { + /* Celestial longitude is given. */ + + /* Check whether the solution interval is a crossing interval. */ + lat0 = span[0]; + world[wcs->lat] = lat0; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj, + imgcrd, lin, pixcrd))) { + return err; + } + d0 = pixcrd[mixpix] - pixmix; + + dabs = fabs(d0); + if (dabs < tol) return 0; + + lat1 = span[1]; + world[wcs->lat] = lat1; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj, + imgcrd, lin, pixcrd))) { + return err; + } + d1 = pixcrd[mixpix] - pixmix; + + dabs = fabs(d1); + if (dabs < tol) return 0; + + lmin = lat1; + dmin = dabs; + + /* Check for a crossing point. */ + if (signb(d0) != signb(d1)) { + crossed = 1; + dx = d1; + } else { + crossed = 0; + lat0 = span[1]; + } + + for (retry = 0; retry < 4; retry++) { + /* Refine the solution interval. */ + while (lat0 > span[0]) { + lat0 -= step; + if (lat0 < span[0]) lat0 = span[0]; + world[wcs->lat] = lat0; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d0 = pixcrd[mixpix] - pixmix; + + /* Check for a solution. */ + dabs = fabs(d0); + if (dabs < tol) return 0; + + /* Record the point of closest approach. */ + if (dabs < dmin) { + lmin = lat0; + dmin = dabs; + } + + /* Check for a crossing point. */ + if (signb(d0) != signb(d1)) { + crossed = 2; + dx = d0; + break; + } + + /* Advance to the next subinterval. */ + lat1 = lat0; + d1 = d0; + } + + if (crossed) { + /* A crossing point was found. */ + for (iter = 0; iter < niter; iter++) { + /* Use regula falsi division of the interval. */ + lambda = d0/(d0-d1); + if (lambda < 0.1) { + lambda = 0.1; + } else if (lambda > 0.9) { + lambda = 0.9; + } + + dlat = lat1 - lat0; + lat = lat0 + lambda*dlat; + world[wcs->lat] = lat; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + + /* Check for a solution. */ + d = pixcrd[mixpix] - pixmix; + dabs = fabs(d); + if (dabs < tol) return 0; + + if (dlat < tol) { + /* An artifact of numerical imprecision. */ + if (dabs < tol2) return 0; + + /* Must be a discontinuity. */ + break; + } + + /* Record the point of closest approach. */ + if (dabs < dmin) { + lmin = lat; + dmin = dabs; + } + + if (signb(d0) == signb(d)) { + lat0 = lat; + d0 = d; + } else { + lat1 = lat; + d1 = d; + } + } + + /* No convergence, must have been a discontinuity. */ + if (crossed == 1) lat0 = span[1]; + lat1 = lat0; + d1 = dx; + crossed = 0; + + } else { + /* No crossing point; look for a tangent point. */ + if (lmin == span[0]) break; + if (lmin == span[1]) break; + + lat = lmin; + lat0 = lat - step; + if (lat0 < span[0]) lat0 = span[0]; + lat1 = lat + step; + if (lat1 > span[1]) lat1 = span[1]; + + world[wcs->lat] = lat0; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d0 = fabs(pixcrd[mixpix] - pixmix); + + d = dmin; + + world[wcs->lat] = lat1; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d1 = fabs(pixcrd[mixpix] - pixmix); + + for (iter = 0; iter < niter; iter++) { + lat0m = (lat0 + lat)/2.0; + world[wcs->lat] = lat0m; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d0m = fabs(pixcrd[mixpix] - pixmix); + + if (d0m < tol) return 0; + + lat1m = (lat1 + lat)/2.0; + world[wcs->lat] = lat1m; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d1m = fabs(pixcrd[mixpix] - pixmix); + + if (d1m < tol) return 0; + + if (d0m < d && d0m <= d1m) { + lat1 = lat; + d1 = d; + lat = lat0m; + d = d0m; + } else if (d1m < d) { + lat0 = lat; + d0 = d; + lat = lat1m; + d = d1m; + } else { + lat0 = lat0m; + d0 = d0m; + lat1 = lat1m; + d1 = d1m; + } + } + } + } + + } else { + /* Celestial latitude is given. */ + + /* Check whether the solution interval is a crossing interval. */ + lng0 = span[0]; + world[wcs->lng] = lng0; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj, + imgcrd, lin, pixcrd))) { + return err; + } + d0 = pixcrd[mixpix] - pixmix; + + dabs = fabs(d0); + if (dabs < tol) return 0; + + lng1 = span[1]; + world[wcs->lng] = lng1; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj, + imgcrd, lin, pixcrd))) { + return err; + } + d1 = pixcrd[mixpix] - pixmix; + + dabs = fabs(d1); + if (dabs < tol) return 0; + lmin = lng1; + dmin = dabs; + + /* Check for a crossing point. */ + if (signb(d0) != signb(d1)) { + crossed = 1; + dx = d1; + } else { + crossed = 0; + lng0 = span[1]; + } + + for (retry = 0; retry < 4; retry++) { + /* Refine the solution interval. */ + while (lng0 > span[0]) { + lng0 -= step; + if (lng0 < span[0]) lng0 = span[0]; + world[wcs->lng] = lng0; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d0 = pixcrd[mixpix] - pixmix; + + /* Check for a solution. */ + dabs = fabs(d0); + if (dabs < tol) return 0; + + /* Record the point of closest approach. */ + if (dabs < dmin) { + lmin = lng0; + dmin = dabs; + } + + /* Check for a crossing point. */ + if (signb(d0) != signb(d1)) { + crossed = 2; + dx = d0; + break; + } + + /* Advance to the next subinterval. */ + lng1 = lng0; + d1 = d0; + } + + if (crossed) { + /* A crossing point was found. */ + for (iter = 0; iter < niter; iter++) { + /* Use regula falsi division of the interval. */ + lambda = d0/(d0-d1); + if (lambda < 0.1) { + lambda = 0.1; + } else if (lambda > 0.9) { + lambda = 0.9; + } + + dlng = lng1 - lng0; + lng = lng0 + lambda*dlng; + world[wcs->lng] = lng; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + + /* Check for a solution. */ + d = pixcrd[mixpix] - pixmix; + dabs = fabs(d); + if (dabs < tol) return 0; + + if (dlng < tol) { + /* An artifact of numerical imprecision. */ + if (dabs < tol2) return 0; + + /* Must be a discontinuity. */ + break; + } + + /* Record the point of closest approach. */ + if (dabs < dmin) { + lmin = lng; + dmin = dabs; + } + + if (signb(d0) == signb(d)) { + lng0 = lng; + d0 = d; + } else { + lng1 = lng; + d1 = d; + } + } + + /* No convergence, must have been a discontinuity. */ + if (crossed == 1) lng0 = span[1]; + lng1 = lng0; + d1 = dx; + crossed = 0; + + } else { + /* No crossing point; look for a tangent point. */ + if (lmin == span[0]) break; + if (lmin == span[1]) break; + + lng = lmin; + lng0 = lng - step; + if (lng0 < span[0]) lng0 = span[0]; + lng1 = lng + step; + if (lng1 > span[1]) lng1 = span[1]; + + world[wcs->lng] = lng0; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d0 = fabs(pixcrd[mixpix] - pixmix); + + d = dmin; + + world[wcs->lng] = lng1; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d1 = fabs(pixcrd[mixpix] - pixmix); + + for (iter = 0; iter < niter; iter++) { + lng0m = (lng0 + lng)/2.0; + world[wcs->lng] = lng0m; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d0m = fabs(pixcrd[mixpix] - pixmix); + + if (d0m < tol) return 0; + + lng1m = (lng1 + lng)/2.0; + world[wcs->lng] = lng1m; + if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, + prj, imgcrd, lin, pixcrd))) { + return err; + } + d1m = fabs(pixcrd[mixpix] - pixmix); + + if (d1m < tol) return 0; + + if (d0m < d && d0m <= d1m) { + lng1 = lng; + d1 = d; + lng = lng0m; + d = d0m; + } else if (d1m < d) { + lng0 = lng; + d0 = d; + lng = lng1m; + d = d1m; + } else { + lng0 = lng0m; + d0 = d0m; + lng1 = lng1m; + d1 = d1m; + } + } + } + } + } + } + + + /* Set cel0 to the unity transformation. */ + cel0.flag = CELSET; + cel0.ref[0] = cel->ref[0]; + cel0.ref[1] = cel->ref[1]; + cel0.ref[2] = cel->ref[2]; + cel0.ref[3] = cel->ref[3]; + cel0.euler[0] = -90.0; + cel0.euler[1] = 0.0; + cel0.euler[2] = 90.0; + cel0.euler[3] = 1.0; + cel0.euler[4] = 0.0; + + /* No convergence, check for aberrant behaviour at a native pole. */ + *theta = -90.0; + for (j = 1; j <= 2; j++) { + /* Could the celestial coordinate element map to a native pole? */ + *theta = -*theta; + err = sphrev(0.0, *theta, cel->euler, &lng, &lat); + + if (mixcel == 1) { + if (fabs(fmod(world[wcs->lng]-lng,360.0)) > tol) continue; + if (lat < span[0]) continue; + if (lat > span[1]) continue; + world[wcs->lat] = lat; + } else { + if (fabs(world[wcs->lat]-lat) > tol) continue; + if (lng < span[0]) lng += 360.0; + if (lng > span[1]) lng -= 360.0; + if (lng < span[0]) continue; + if (lng > span[1]) continue; + world[wcs->lng] = lng; + } + + /* Is there a solution for the given pixel coordinate element? */ + lng = world[wcs->lng]; + lat = world[wcs->lat]; + + /* Feed native coordinates to wcsfwd() with cel0 set to unity. */ + world[wcs->lng] = -180.0; + world[wcs->lat] = *theta; + if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj, + imgcrd, lin, pixcrd))) { + return err; + } + d0 = pixcrd[mixpix] - pixmix; + + /* Check for a solution. */ + if (fabs(d0) < tol) { + /* Recall saved world coordinates. */ + world[wcs->lng] = lng; + world[wcs->lat] = lat; + return 0; + } + + /* Search for a crossing interval. */ + phi0 = -180.0; + for (k = -179; k <= 180; k++) { + phi1 = (double) k; + world[wcs->lng] = phi1; + if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj, + imgcrd, lin, pixcrd))) { + return err; + } + d1 = pixcrd[mixpix] - pixmix; + + /* Check for a solution. */ + dabs = fabs(d1); + if (dabs < tol) { + /* Recall saved world coordinates. */ + world[wcs->lng] = lng; + world[wcs->lat] = lat; + return 0; + } + + /* Is it a crossing interval? */ + if (signb(d0) != signb(d1)) break; + + phi0 = phi1; + d0 = d1; + } + + for (iter = 1; iter <= niter; iter++) { + /* Use regula falsi division of the interval. */ + lambda = d0/(d0-d1); + if (lambda < 0.1) { + lambda = 0.1; + } else if (lambda > 0.9) { + lambda = 0.9; + } + + dphi = phi1 - phi0; + world[wcs->lng] = phi0 + lambda*dphi; + if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj, + imgcrd, lin, pixcrd))) { + return err; + } + + /* Check for a solution. */ + d = pixcrd[mixpix] - pixmix; + dabs = fabs(d); + if (dabs < tol || (dphi < tol && dabs < tol2)) { + /* Recall saved world coordinates. */ + world[wcs->lng] = lng; + world[wcs->lat] = lat; + return 0; + } + + if (signb(d0) == signb(d)) { + phi0 = world[wcs->lng]; + d0 = d; + } else { + phi1 = world[wcs->lng]; + d1 = d; + } + } + } + + + /* No solution. */ + return 5; + +} +/* Dec 20 1999 Doug Mink - Change signbit() to signb() and always define it + * Dec 20 1999 Doug Mink - Include wcslib.h, which includes wcs.h, wcstrig.h + * + * Mar 20 2001 Doug Mink - Include stdio.h for sprintf() + * Mar 20 2001 Doug Mink - Add () around err assignments in if statements + * Sep 19 2001 Doug Mink - Add above changes to WCSLIB-2.7 version + * + * Mar 15 2002 Doug Mink - Add above changes to WCSLIB-2.8.2 + * Apr 3 2002 Mark Calabretta - Fix bug in code checking section + * + * Jun 20 2006 Doug Mink - Initialized uninitialized variables + */ |