Skip to content
GitLab
Commit 3b348a94 authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Changed trunk directory name

parents
Hide whitespace changes
Inline Side-by-side
src/types.h 0 → 100644
View file @ 3b348a94
/*
types.h
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: SExtractor
*
* Author: E.BERTIN (IAP)
*
* Contents: global type definitions.
*
* Last modify: 12/01/2006
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#include <stdio.h>
#ifndef _FITSCAT_H_
#include "fits/fitscat.h"
#endif
/*-------------------------------- flags ------------------------------------*/
#define OBJ_CROWDED 0x0001
#define OBJ_MERGED 0x0002
#define OBJ_SATUR 0x0004
#define OBJ_TRUNC 0x0008
#define OBJ_APERT_PB 0x0010
#define OBJ_ISO_PB 0x0020
#define OBJ_DOVERFLOW 0x0040
#define OBJ_OVERFLOW 0x0080
/*----------------------------- weight flags --------------------------------*/
#define OBJ_WEIGHTZERO 0x0001
#define OBJ_DWEIGHTZERO 0x0002
/*---------------------------- preanalyse flags -----------------------------*/
#define ANALYSE_FAST 0
#define ANALYSE_FULL 1
#define ANALYSE_ROBUST 2
/*--------------------------------- typedefs --------------------------------*/
typedef unsigned char BYTE; /* a byte */
typedef unsigned short USHORT; /* 0 to 65535 integers */
typedef unsigned int FLAGTYPE; /* flag type */
typedef char pliststruct; /* Dummy type for plist */
typedef int LONG;
typedef unsigned int ULONG;
typedef enum {BACK_RELATIVE, BACK_ABSOLUTE}
backenum; /* BACK_TYPE */
typedef enum {CHECK_NONE, CHECK_IDENTICAL, CHECK_BACKGROUND,
CHECK_BACKRMS, CHECK_MINIBACKGROUND, CHECK_MINIBACKRMS,
CHECK_SUBTRACTED, CHECK_FILTERED, CHECK_OBJECTS, CHECK_APERTURES,
CHECK_SEGMENTATION, CHECK_ASSOC, CHECK_SUBOBJECTS,
CHECK_SUBPSFPROTOS, CHECK_PSFPROTOS,
CHECK_SUBPCPROTOS, CHECK_PCPROTOS, CHECK_PCOPROTOS,
CHECK_MAPSOM} checkenum;
/* CHECK_IMAGE type */
typedef enum {WEIGHT_NONE, WEIGHT_FROMBACK, WEIGHT_FROMRMSMAP,
WEIGHT_FROMVARMAP, WEIGHT_FROMWEIGHTMAP, WEIGHT_FROMINTERP}
weightenum; /* WEIGHT_IMAGE type */
/*--------------------------------- objects ---------------------------------*/
/* I: "PIXEL" parameters */
typedef struct
{
/* ---- basic parameters */
int number; /* ID */
int fdnpix; /* nb of extracted pix */
int dnpix; /* nb of pix above thresh */
int npix; /* "" in measured frame */
float fdflux; /* integrated ext. flux */
float dflux; /* integrated det. flux */
float flux; /* integrated mes. flux */
float fluxerr; /* integrated variance */
float flux_prof; /* PROFILE flux*/
float fluxerr_prof; /* PROFILE flux variance */
PIXTYPE fdpeak; /* peak intensity (ADU) */
PIXTYPE dpeak; /* peak intensity (ADU) */
PIXTYPE peak; /* peak intensity (ADU) */
/* ---- astrometric data */
int peakx,peaky; /* pos of brightest pix */
double mx, my; /* barycenter */
double poserr_mx2, poserr_my2,
poserr_mxy; /* Error ellips moments */
/* ---- morphological data */
int xmin,xmax,ymin,ymax,ycmin,ycmax;/* x,y limits */
PIXTYPE *blank, *dblank; /* BLANKing sub-images */
int *submap; /* Pixel-index sub-map */
int subx,suby, subw,subh; /* sub-image pos. and size */
short flag; /* extraction flags */
BYTE wflag; /* weighted extraction flags */
FLAGTYPE imaflag[MAXFLAG]; /* flags from FLAG-images */
BYTE singuflag; /* flags for singularities */
int imanflag[MAXFLAG]; /* number of MOST flags */
double mx2,my2,mxy; /* variances and covariance */
float a, b, theta, abcor; /* moments and angle */
float cxx,cyy,cxy; /* ellipse parameters */
int firstpix; /* ptr to first pixel */
int lastpix; /* ptr to last pixel */
float bkg, dbkg, sigbkg; /* Background stats (ADU) */
float thresh; /* measur. threshold (ADU) */
float dthresh; /* detect. threshold (ADU) */
float mthresh; /* max. threshold (ADU) */
int iso[NISO]; /* isophotal areas */
float fwhm; /* IMAGE FWHM */
} objstruct;
/* II: "BLIND" parameters */
typedef struct
{
/* ---- photometric data */
float flux_iso; /* ISO integrated flux */
float fluxerr_iso; /* RMS error on ISO flux */
float mag_iso; /* ISO mag */
float magerr_iso; /* ISO mag uncertainty */
float flux_isocor; /* ISOCOR integrated flux */
float fluxerr_isocor; /* RMS error on ISOCOR flux */
float mag_isocor; /* ISOCOR mag */
float magerr_isocor; /* ISOCOR mag uncertainty */
float kronfactor; /* kron parameter */
float flux_auto; /* AUTO integrated flux */
float fluxerr_auto; /* RMS error on AUTO flux */
float mag_auto; /* AUTO mag */
float magerr_auto; /* AUTO mag uncertainty */
float petrofactor; /* kron parameter */
float flux_petro; /* AUTO integrated flux */
float fluxerr_petro; /* RMS error on AUTO flux */
float mag_petro; /* AUTO mag */
float magerr_petro; /* AUTO mag uncertainty */
float flux_best; /* BEST integrated flux */
float fluxerr_best; /* RMS error on BEST flux */
float mag_best; /* BEST mag */
float magerr_best; /* BEST mag uncertainty */
float *flux_aper; /* APER flux vector */
float *fluxerr_aper; /* APER flux error vector */
float *mag_aper; /* APER magnitude vector */
float *magerr_aper; /* APER mag error vector */
float flux_prof; /* PROFILE flux*/
float fluxerr_prof; /* PROFILE flux error */
float mag_prof; /* PROFILE magnitude */
float magerr_prof; /* PROFILE magnitude error */
float flux_win; /* WINdowed flux*/
float fluxerr_win; /* WINdowed flux error */
float mag_win; /* WINdowed magnitude */
float magerr_win; /* WINdowed magnitude error */
/* ---- astrometric data */
double posx,posy; /* "FITS" pos. in pixels */
double mamaposx,mamaposy; /* "MAMA" pos. in pixels */
float sposx,sposy; /* single precision pos. */
float poserr_a, poserr_b,
poserr_theta; /* Error ellips parameters */
float poserr_cxx, poserr_cyy,
poserr_cxy; /* pos. error ellipse */
double poserr_mx2w, poserr_my2w,
poserr_mxyw; /* WORLD error moments */
float poserr_aw, poserr_bw,
poserr_thetaw; /* WORLD error parameters */
float poserr_thetas; /* native error pos. angle */
float poserr_theta2000; /* J2000 error pos. angle */
float poserr_theta1950; /* B1950 error pos. angle */
float poserr_cxxw, poserr_cyyw,
poserr_cxyw; /* WORLD error ellipse */
double mx2w,my2w,mxyw; /* WORLD var. and covar. */
double peakxw, peakyw; /* WORLD of brightest pix */
double mxw, myw; /* WORLD barycenters */
double alphas, deltas; /* native alpha, delta */
float thetas; /* native position angle E/N*/
double peakalphas, peakdeltas; /* native for brightest pix */
double peakalpha2000, peakdelta2000; /* J2000 for brightest pix */
double peakalpha1950, peakdelta1950; /* B1950 for brightest pix */
double alpha2000, delta2000; /* J2000 alpha, delta */
float theta2000; /* J2000 position angle E/N */
double alpha1950, delta1950; /* B1950 alpha, delta */
float theta1950; /* B1950 position angle E/N */
float aw, bw; /* WORLD ellipse size */
float thetaw; /* WORLD position angle */
float cxxw,cyyw,cxyw; /* WORLD ellipse parameters */
float npixw, fdnpixw; /* WORLD isophotal areas */
float threshmu; /* det. surface brightnees */
float maxmu; /* max. surface brightnees */
float elong; /* elongation */
float ellip; /* ellipticity */
float polar; /* Kaiser's "polarization" */
float polarw; /* WORLD "polarization" */
float sprob; /* Stellarity index */
float fwhmw; /* WORLD FWHM */
float *assoc; /* ASSOCiated data */
int assoc_number; /* nb of ASSOCiated objects */
float *vignet; /* Pixel data */
float *vigshift; /* (Shifted) pixel data */
/* Windowed measurements */
double winpos_x,winpos_y; /* Windowed barycenter */
double winposerr_mx2, winposerr_my2,
winposerr_mxy; /* Error ellips moments */
float winposerr_a, winposerr_b,
winposerr_theta; /* Error ellips parameters */
float winposerr_cxx, winposerr_cyy,
winposerr_cxy; /* pos. error ellipse */
double winposerr_mx2w, winposerr_my2w,
winposerr_mxyw; /* WORLD error moments */
float winposerr_aw, winposerr_bw,
winposerr_thetaw; /* WORLD error parameters */
float winposerr_thetas; /* native error pos. angle */
float winposerr_theta2000; /* J2000 error pos. angle */
float winposerr_theta1950; /* B1950 error pos. angle */
float winposerr_cxxw, winposerr_cyyw,
winposerr_cxyw; /* WORLD error ellipse */
double win_mx2, win_my2,
win_mxy; /* Windowed moments */
float win_a, win_b,
win_theta; /* Windowed ellipse parameters*/
float win_polar; /* Windowed "polarization" */
float win_cxx, win_cyy,
win_cxy; /* Windowed ellipse parameters*/
double win_mx2w, win_my2w,
win_mxyw; /* WORLD windowed moments */
float win_aw, win_bw,
win_thetaw; /* WORLD ellipse parameters */
float win_polarw; /* WORLD WIN "polarization" */
float win_thetas; /* native error pos. angle */
float win_theta2000; /* J2000 error pos. angle */
float win_theta1950; /* B1950 error pos. angle */
float win_cxxw, win_cyyw,
win_cxyw; /* WORLD ellipse parameters */
double winpos_xw, winpos_yw; /* WORLD coordinates */
double winpos_alphas, winpos_deltas; /* native alpha, delta */
double winpos_alpha2000, winpos_delta2000; /* J2000 alpha, delta */
double winpos_alpha1950, winpos_delta1950; /* B1950 alpha, delta */
short winpos_niter; /* Number of WIN iterations */
short win_flag; /* 1:x2<0 2:xy=x2 4:flux<0 */
/* ---- SOM fitting */
float flux_somfit; /* Fitted amplitude */
float fluxerr_somfit; /* RMS error on SOM flux */
float mag_somfit; /* Magnitude from SOM fit */
float magerr_somfit; /* Mag. err. from SOM fit */
float stderr_somfit; /* Fitting reduced error */
float *vector_somfit; /* SOM fit vector position */
/* ---- Growth curves and stuff */
float *flux_growth; /* Cumulated growth_curve */
float flux_growthstep; /* Growth-curve step */
float *mag_growth; /* Cumulated growth_curve */
float mag_growthstep; /* Growth-curve step */
float *flux_radius; /* f-light-radii */
float hl_radius; /* Scalar half-light radius */
/* ---- PSF-fitting */
float *flux_psf; /* Flux from PSF-fitting */
float *fluxerr_psf; /* RMS error on PSF flux */
float *mag_psf; /* Mag from PSF-fitting */
float *magerr_psf; /* RMS mag from PSF-fitting */
float *x_psf, *y_psf; /* Coords from PSF-fitting */
short niter_psf; /* # of PSF-fitting iterat. */
short npsf; /* # of fitted PSFs */
float chi2_psf; /* Red. chi2 of PSF-fitting */
double xw_psf, yw_psf; /* WORLD coords */
double alphas_psf, deltas_psf; /* native alpha, delta */
double alpha2000_psf, delta2000_psf; /* J2000 alpha, delta */
double alpha1950_psf, delta1950_psf; /* B1950 alpha, delta */
double poserrmx2_psf, poserrmy2_psf,
poserrmxy_psf; /* Error ellips moments */
float poserra_psf, poserrb_psf,
poserrtheta_psf; /* Error ellips parameters */
float poserrcxx_psf, poserrcyy_psf,
poserrcxy_psf; /* pos. error ellipse */
double poserrmx2w_psf, poserrmy2w_psf,
poserrmxyw_psf; /* WORLD error moments */
float poserraw_psf, poserrbw_psf,
poserrthetaw_psf; /* WORLD error parameters */
float poserrthetas_psf; /* native error pos. angle */
float poserrtheta2000_psf; /* J2000 error pos. angle */
float poserrtheta1950_psf; /* B1950 error pos. angle */
float poserrcxxw_psf, poserrcyyw_psf,
poserrcxyw_psf; /* WORLD error ellipse */
/* ---- PC-fitting */
double mx2_pc,my2_pc,mxy_pc; /* PC 2nd-order parameters */
float a_pc,b_pc,theta_pc; /* PC shape parameters */
float *vector_pc; /* Principal components */
float gdposang; /* Gal. disk position angle */
float gdscale; /* Gal. disk scalelength */
float gdaspect; /* Gal. disk aspect-ratio */
float gde1,gde2; /* Gal. disk ellipticities */
float gbratio; /* Galaxy B/T */
float gbposang; /* Gal. bulge position angle */
float gbscale; /* Gal. bulge scalelength */
float gbaspect; /* Gal. bulge aspect-ratio */
float gbe1,gbe2; /* Gal. bulge ellipticities */
float flux_galfit; /* Galaxy tot. flux from fit */
float fluxerr_galfit; /* RMS error on galfit flux */
float mag_galfit; /* Galaxy tot. mag from fit */
float magerr_galfit; /* RMS error on galfit mag */
/* ---- MEF */
short ext_number; /* FITS extension number */
} obj2struct;
/*----------------------------- lists of objects ----------------------------*/
typedef struct
{
int nobj; /* number of objects in list */
objstruct *obj; /* pointer to the object array */
int npix; /* number of pixels in pixel-list */
pliststruct *plist; /* pointer to the pixel-list */
PIXTYPE dthresh; /* detection threshold */
PIXTYPE thresh; /* analysis threshold */
} objliststruct;
/*----------------------------- image parameters ----------------------------*/
typedef struct pic
{
char filename[MAXCHAR]; /* pointer to the image filename */
char *rfilename; /* pointer to the reduced image name */
char ident[MAXCHAR]; /* field identifier (read from FITS)*/
char rident[MAXCHAR]; /* field identifier (relative) */
FILE *file; /* pointer the image file structure */
char *fitshead; /* pointer to the FITS header */
int fitsheadsize; /* FITS header size */
/* ---- main image parameters */
int bitpix, bytepix; /* nb of bits and bytes per pixel */
int bitsgn; /* non-zero if signed integer data */
int width, height; /* x,y size of the field */
KINGSIZE_T npix; /* total number of pixels */
double bscale, bzero; /* FITS scale and offset */
double ngamma; /* normalized photo gamma */
int nlevels; /* nb of quantification levels */
float pixmin, pixmax; /* min and max values in frame */
int y; /* y current position in field */
int ymin; /* y limit (lowest accessible) */
int ymax; /* y limit (highest accessible+1) */
int yblank; /* y blanking limit (highest+1) */
PIXTYPE *strip; /* pointer to the image buffer */
FLAGTYPE *fstrip; /* pointer to the FLAG buffer */
int stripheight; /* height of a strip (in lines) */
int stripmargin; /* number of lines in margin */
int stripstep; /* number of lines at each read */
int stripy; /* y position in buffer */
int stripylim; /* y limit in buffer */
int stripysclim; /* y scroll limit in buffer */
/* ---- image (de-)compression */
enum {ICOMPRESS_NONE, ICOMPRESS_BASEBYTE, ICOMPRESS_PREVPIX}
compress_type; /* image compression type */
char *compress_buf; /* de-compression buffer */
char *compress_bufptr; /* present pixel in buffer */
int compress_curval; /* current pixel or checksum value */
int compress_checkval; /* foreseen pixel or checksum value */
int compress_npix; /* remaining pixels in buffer */
/* ---- basic astrometric parameters */
double pixscale; /* pixel size in arcsec.pix-1 */
double epoch; /* epoch of coordinates */
/* ---- background parameters */
float *back; /* ptr to the background map in mem */
float *dback; /* ptr to the background deriv. map */
float *sigma; /* ptr to the sigma map */
float *dsigma; /* Ptr to the sigma deriv. map */
int backw, backh; /* x,y size of a bkgnd mesh */
int nbackp; /* total nb of pixels per bkgnd mesh */
int nbackx, nbacky; /* x,y number of bkgnd meshes */
int nback; /* total number of bkgnd meshes */
int nbackfx, nbackfy; /* x,y size of bkgnd filtering mask */
float backmean; /* median bkgnd value in image */
float backsig; /* median bkgnd rms in image */
float sigfac; /* scaling RMS factor (for WEIGHTs) */
PIXTYPE *backline; /* current interpolated bkgnd line */
PIXTYPE dthresh; /* detection threshold */
PIXTYPE thresh; /* analysis threshold */
backenum back_type; /* Background type */
/* ---- astrometric parameters */
struct structastrom *astrom; /* astrometric data */
struct structassoc *assoc; /* ptr to the assoc-list */
int flags; /* flags defining the field type */
/* ---- image interpolation */
int interp_flag; /* interpolation for this field? */
PIXTYPE *interp_backup; /* backup line for interpolation */
PIXTYPE weight_thresh; /* interpolation threshold */
int *interp_ytimeoutbuf; /* interpolation timeout line buffer */
int interp_xtimeout; /* interpolation timeout value in x */
int interp_ytimeout; /* interpolation timeout value in y */
struct pic *reffield; /* pointer to a reference field */
OFF_T mefpos; /* Position in a MEF file */
} picstruct;
/*-------------------------------- catalog ---------------------------------*/
typedef struct
{
int ndetect; /* nb of detections */
int ntotal; /* Total object nb */
int nparam; /* Nb of parameters */
/*----- Misc. strings defining the extraction */
char prefs_name[MAXCHAR]; /* Prefs filename*/
char image_name[MAXCHAR]; /* image filename*/
char psf_name[MAXCHAR]; /* PSF filename*/
char nnw_name[MAXCHAR]; /* NNW name */
char filter_name[MAXCHAR]; /* Filter name */
char soft_name[MAXCHAR]; /* Sextractor version*/
/*----- time */
char ext_date[16],ext_time[16]; /* date and time */
double ext_elapsed; /* processing time */
/*----- MEF */
int currext; /* current extension */
int next; /* Nb of extensions */
} sexcatstruct;
src/wcs/Makefile.am 0 → 100644
View file @ 3b348a94
# Program Makefile for the WCS library
# Copyright (C) 2002 Emmanuel Bertin.
noinst_LIBRARIES = libwcs_c.a
libwcs_c_a_SOURCES = cel.c lin.c poly.c proj.c sph.c tnx.c wcs.c \
wcstrig.c \
cel.h lin.h poly.h proj.h sph.h tnx.h wcs.h \
wcsmath.h wcstrig.h
src/wcs/Makefile.in 0 → 100644
View file @ 3b348a94
# Makefile.in generated by automake 1.9.4 from Makefile.am.
# @configure_input@
# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
# 2003, 2004 Free Software Foundation, Inc.
# This Makefile.in is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY, to the extent permitted by law; without
# even the implied warranty of MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE.
@SET_MAKE@
SOURCES = $(libwcs_c_a_SOURCES)
srcdir = @srcdir@
top_srcdir = @top_srcdir@
VPATH = @srcdir@
pkgdatadir = $(datadir)/@PACKAGE@
pkglibdir = $(libdir)/@PACKAGE@
pkgincludedir = $(includedir)/@PACKAGE@
top_builddir = ../..
am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd
INSTALL = @INSTALL@
install_sh_DATA = $(install_sh) -c -m 644
install_sh_PROGRAM = $(install_sh) -c
install_sh_SCRIPT = $(install_sh) -c
INSTALL_HEADER = $(INSTALL_DATA)
transform = $(program_transform_name)
NORMAL_INSTALL = :
PRE_INSTALL = :
POST_INSTALL = :
NORMAL_UNINSTALL = :
PRE_UNINSTALL = :
POST_UNINSTALL = :
subdir = src/wcs
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/acx_urbi_resolve_dir.m4 \
$(top_srcdir)/acx_prog_cc_optim.m4 \
$(top_srcdir)/acx_prog_cc_optim.m4 \
$(top_srcdir)/acx_urbi_resolve_dir.m4 \
$(top_srcdir)/configure.ac
am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \
$(ACLOCAL_M4)
mkinstalldirs = $(SHELL) $(top_srcdir)/autoconf/mkinstalldirs
CONFIG_HEADER = $(top_builddir)/config.h
CONFIG_CLEAN_FILES =
LIBRARIES = $(noinst_LIBRARIES)
AR = ar
ARFLAGS = cru
libwcs_c_a_AR = $(AR) $(ARFLAGS)
libwcs_c_a_LIBADD =
am_libwcs_c_a_OBJECTS = cel.$(OBJEXT) lin.$(OBJEXT) poly.$(OBJEXT) \
proj.$(OBJEXT) sph.$(OBJEXT) tnx.$(OBJEXT) wcs.$(OBJEXT) \
wcstrig.$(OBJEXT)
libwcs_c_a_OBJECTS = $(am_libwcs_c_a_OBJECTS)
DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir)
depcomp = $(SHELL) $(top_srcdir)/autoconf/depcomp
am__depfiles_maybe = depfiles
COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS)
CCLD = $(CC)
LINK = $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@
SOURCES = $(libwcs_c_a_SOURCES)
DIST_SOURCES = $(libwcs_c_a_SOURCES)
ETAGS = etags
CTAGS = ctags
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
ACLOCAL = @ACLOCAL@
AMDEP_FALSE = @AMDEP_FALSE@
AMDEP_TRUE = @AMDEP_TRUE@
AMTAR = @AMTAR@
AUTOCONF = @AUTOCONF@
AUTOHEADER = @AUTOHEADER@
AUTOMAKE = @AUTOMAKE@
AWK = @AWK@
CC = @CC@
CCDEPMODE = @CCDEPMODE@
CFLAGS = @CFLAGS@
CPP = @CPP@
CPPFLAGS = @CPPFLAGS@
CYGPATH_W = @CYGPATH_W@
DATE2 = @DATE2@
DATE3 = @DATE3@
DEFS = @DEFS@
DEPDIR = @DEPDIR@
ECHO_C = @ECHO_C@
ECHO_N = @ECHO_N@
ECHO_T = @ECHO_T@
EGREP = @EGREP@
EXEEXT = @EXEEXT@
INSTALL_DATA = @INSTALL_DATA@
INSTALL_PROGRAM = @INSTALL_PROGRAM@
INSTALL_SCRIPT = @INSTALL_SCRIPT@
INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@
LDFLAGS = @LDFLAGS@
LIBOBJS = @LIBOBJS@
LIBS = @LIBS@
LTLIBOBJS = @LTLIBOBJS@
MAKEINFO = @MAKEINFO@
OBJEXT = @OBJEXT@
PACKAGE = @PACKAGE@
PACKAGER = @PACKAGER@
PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@
PACKAGE_NAME = @PACKAGE_NAME@
PACKAGE_STRING = @PACKAGE_STRING@
PACKAGE_TARNAME = @PACKAGE_TARNAME@
PACKAGE_VERSION = @PACKAGE_VERSION@
PATH_SEPARATOR = @PATH_SEPARATOR@
RANLIB = @RANLIB@
SET_MAKE = @SET_MAKE@
SHELL = @SHELL@
STRIP = @STRIP@
VERSION = @VERSION@
ac_ct_CC = @ac_ct_CC@
ac_ct_RANLIB = @ac_ct_RANLIB@
ac_ct_STRIP = @ac_ct_STRIP@
am__fastdepCC_FALSE = @am__fastdepCC_FALSE@
am__fastdepCC_TRUE = @am__fastdepCC_TRUE@
am__include = @am__include@
am__leading_dot = @am__leading_dot@
am__quote = @am__quote@
am__tar = @am__tar@
am__untar = @am__untar@
bindir = @bindir@
build_alias = @build_alias@
datadir = @datadir@
exec_prefix = @exec_prefix@
host_alias = @host_alias@
includedir = @includedir@
infodir = @infodir@
install_sh = @install_sh@
libdir = @libdir@
libexecdir = @libexecdir@
localstatedir = @localstatedir@
mandir = @mandir@
mkdir_p = @mkdir_p@
oldincludedir = @oldincludedir@
prefix = @prefix@
program_transform_name = @program_transform_name@
sbindir = @sbindir@
sharedstatedir = @sharedstatedir@
sysconfdir = @sysconfdir@
target_alias = @target_alias@
# Program Makefile for the WCS library
# Copyright (C) 2002 Emmanuel Bertin.
noinst_LIBRARIES = libwcs_c.a
libwcs_c_a_SOURCES = cel.c lin.c poly.c proj.c sph.c tnx.c wcs.c \
wcstrig.c \
cel.h lin.h poly.h proj.h sph.h tnx.h wcs.h \
wcsmath.h wcstrig.h
all: all-am
.SUFFIXES:
.SUFFIXES: .c .o .obj
$(srcdir)/Makefile.in: $(srcdir)/Makefile.am $(am__configure_deps)
@for dep in $?; do \
case '$(am__configure_deps)' in \
*$$dep*) \
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh \
&& exit 0; \
exit 1;; \
esac; \
done; \
echo ' cd $(top_srcdir) && $(AUTOMAKE) --foreign src/wcs/Makefile'; \
cd $(top_srcdir) && \
$(AUTOMAKE) --foreign src/wcs/Makefile
.PRECIOUS: Makefile
Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status
@case '$?' in \
*config.status*) \
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \
*) \
echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \
cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \
esac;
$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES)
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
$(top_srcdir)/configure: $(am__configure_deps)
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
$(ACLOCAL_M4): $(am__aclocal_m4_deps)
cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh
clean-noinstLIBRARIES:
-test -z "$(noinst_LIBRARIES)" || rm -f $(noinst_LIBRARIES)
libwcs_c.a: $(libwcs_c_a_OBJECTS) $(libwcs_c_a_DEPENDENCIES)
-rm -f libwcs_c.a
$(libwcs_c_a_AR) libwcs_c.a $(libwcs_c_a_OBJECTS) $(libwcs_c_a_LIBADD)
$(RANLIB) libwcs_c.a
mostlyclean-compile:
-rm -f *.$(OBJEXT)
distclean-compile:
-rm -f *.tab.c
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cel.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lin.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/poly.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/proj.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/sph.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/tnx.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/wcs.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/wcstrig.Po@am__quote@
.c.o:
@am__fastdepCC_TRUE@ if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ $<; \
@am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/$*.Tpo" "$(DEPDIR)/$*.Po"; else rm -f "$(DEPDIR)/$*.Tpo"; exit 1; fi
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(COMPILE) -c $<
.c.obj:
@am__fastdepCC_TRUE@ if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ `$(CYGPATH_W) '$<'`; \
@am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/$*.Tpo" "$(DEPDIR)/$*.Po"; else rm -f "$(DEPDIR)/$*.Tpo"; exit 1; fi
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(COMPILE) -c `$(CYGPATH_W) '$<'`
uninstall-info-am:
ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES)
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) ' { files[$$0] = 1; } \
END { for (i in files) print i; }'`; \
mkid -fID $$unique
tags: TAGS
TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
$(TAGS_FILES) $(LISP)
tags=; \
here=`pwd`; \
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) ' { files[$$0] = 1; } \
END { for (i in files) print i; }'`; \
if test -z "$(ETAGS_ARGS)$$tags$$unique"; then :; else \
test -n "$$unique" || unique=$$empty_fix; \
$(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \
$$tags $$unique; \
fi
ctags: CTAGS
CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
$(TAGS_FILES) $(LISP)
tags=; \
here=`pwd`; \
list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
unique=`for i in $$list; do \
if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
done | \
$(AWK) ' { files[$$0] = 1; } \
END { for (i in files) print i; }'`; \
test -z "$(CTAGS_ARGS)$$tags$$unique" \
|| $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \
$$tags $$unique
GTAGS:
here=`$(am__cd) $(top_builddir) && pwd` \
&& cd $(top_srcdir) \
&& gtags -i $(GTAGS_ARGS) $$here
distclean-tags:
-rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags
distdir: $(DISTFILES)
@srcdirstrip=`echo "$(srcdir)" | sed 's|.|.|g'`; \
topsrcdirstrip=`echo "$(top_srcdir)" | sed 's|.|.|g'`; \
list='$(DISTFILES)'; for file in $$list; do \
case $$file in \
$(srcdir)/*) file=`echo "$$file" | sed "s|^$$srcdirstrip/||"`;; \
$(top_srcdir)/*) file=`echo "$$file" | sed "s|^$$topsrcdirstrip/|$(top_builddir)/|"`;; \
esac; \
if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \
dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \
if test "$$dir" != "$$file" && test "$$dir" != "."; then \
dir="/$$dir"; \
$(mkdir_p) "$(distdir)$$dir"; \
else \
dir=''; \
fi; \
if test -d $$d/$$file; then \
if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \
cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \
fi; \
cp -pR $$d/$$file $(distdir)$$dir || exit 1; \
else \
test -f $(distdir)/$$file \
|| cp -p $$d/$$file $(distdir)/$$file \
|| exit 1; \
fi; \
done
check-am: all-am
check: check-am
all-am: Makefile $(LIBRARIES)
installdirs:
install: install-am
install-exec: install-exec-am
install-data: install-data-am
uninstall: uninstall-am
install-am: all-am
@$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am
installcheck: installcheck-am
install-strip:
$(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \
install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \
`test -z '$(STRIP)' || \
echo "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'"` install
mostlyclean-generic:
clean-generic:
distclean-generic:
-test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES)
maintainer-clean-generic:
@echo "This command is intended for maintainers to use"
@echo "it deletes files that may require special tools to rebuild."
clean: clean-am
clean-am: clean-generic clean-noinstLIBRARIES mostlyclean-am
distclean: distclean-am
-rm -rf ./$(DEPDIR)
-rm -f Makefile
distclean-am: clean-am distclean-compile distclean-generic \
distclean-tags
dvi: dvi-am
dvi-am:
html: html-am
info: info-am
info-am:
install-data-am:
install-exec-am:
install-info: install-info-am
install-man:
installcheck-am:
maintainer-clean: maintainer-clean-am
-rm -rf ./$(DEPDIR)
-rm -f Makefile
maintainer-clean-am: distclean-am maintainer-clean-generic
mostlyclean: mostlyclean-am
mostlyclean-am: mostlyclean-compile mostlyclean-generic
pdf: pdf-am
pdf-am:
ps: ps-am
ps-am:
uninstall-am: uninstall-info-am
.PHONY: CTAGS GTAGS all all-am check check-am clean clean-generic \
clean-noinstLIBRARIES ctags distclean distclean-compile \
distclean-generic distclean-tags distdir dvi dvi-am html \
html-am info info-am install install-am install-data \
install-data-am install-exec install-exec-am install-info \
install-info-am install-man install-strip installcheck \
installcheck-am installdirs maintainer-clean \
maintainer-clean-generic mostlyclean mostlyclean-compile \
mostlyclean-generic pdf pdf-am ps ps-am tags uninstall \
uninstall-am uninstall-info-am
# Tell versions [3.59,3.63) of GNU make to not export all variables.
# Otherwise a system limit (for SysV at least) may be exceeded.
.NOEXPORT:
src/wcs/cel.c 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
* -------------------
* These routines are provided as drivers for the lower level spherical
* coordinate transformation and projection routines. There are separate
* driver routines for the forward, celfwd(), and reverse, celrev(),
* transformations.
*
* An initialization routine, celset(), computes intermediate values from
* the transformation parameters but need not be called explicitly - see the
* explanation of cel.flag below.
*
*
* Initialization routine; celset()
* --------------------------------
* Initializes members of a celprm data structure which hold intermediate
* values. Note that this routine need not be called directly; it will be
* invoked by celfwd() and celrev() if the "flag" structure member is
* anything other than a predefined magic value.
*
* Given:
* pcode[4] const char
* WCS projection code (see below).
*
* Given and returned:
* cel celprm* Spherical coordinate transformation parameters
* (see below).
* prj prjprm* Projection parameters (usage is described in the
* prologue to "proj.c").
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid coordinate transformation parameters.
* 2: Ill-conditioned coordinate transformation
* parameters.
*
* Forward transformation; celfwd()
* --------------------------------
* Compute (x,y) coordinates in the plane of projection from celestial
* coordinates (lng,lat).
*
* Given:
* pcode[4] const char
* WCS projection code (see below).
* lng,lat const double
* Celestial longitude and latitude of the projected
* point, in degrees.
*
* Given and returned:
* cel celprm* Spherical coordinate transformation parameters
* (see below).
*
* 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:
* x,y double* Projected coordinates, "degrees".
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid coordinate transformation parameters.
* 2: Invalid projection parameters.
* 3: Invalid value of (lng,lat).
*
* Reverse transformation; celrev()
* --------------------------------
* Compute the celestial coordinates (lng,lat) of the point with projected
* coordinates (x,y).
*
* Given:
* pcode[4] const char
* WCS projection code (see below).
* x,y const double
* Projected coordinates, "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 and returned:
* cel celprm* Spherical coordinate transformation parameters
* (see below).
*
* Returned:
* lng,lat double* Celestial longitude and latitude of the projected
* point, in degrees.
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid coordinate transformation parameters.
* 2: Invalid projection parameters.
* 3: Invalid value of (x,y).
*
* Coordinate transformation parameters
* ------------------------------------
* The celprm struct consists of the following:
*
* int flag
* The celprm struct contains pointers to the forward and reverse
* projection routines as well as intermediaries computed from the
* reference coordinates (see below). Whenever the projection code
* (pcode) or any of ref[4] are set or changed then this flag must be
* set to zero to signal the initialization routine, celset(), to
* redetermine the function pointers and recompute intermediaries.
* Once this has been done pcode itself is ignored.
*
* double ref[4]
* The first pair of values should be set to the celestial longitude
* and latitude (usually right ascension and declination) of the
* reference point of the projection.
*
* The second pair of values are the native longitude and latitude of
* the pole of the celestial coordinate system and correspond to the
* FITS keywords LONGPOLE and LATPOLE.
*
* LONGPOLE defaults to 0 degrees if the celestial latitude of the
* reference point of the projection is greater than the native
* latitude, otherwise 180 degrees. (This is the condition for the
* celestial latitude to increase in the same direction as the native
* latitude at the reference point.) ref[2] may be set to 999.0 to
* indicate that the correct default should be substituted.
*
* In some circumstances the latitude of the native pole may be
* determined by the first three values only to within a sign and
* LATPOLE is used to choose between the two solutions. LATPOLE is
* set in ref[3] and the solution closest to this value is used to
* reset ref[3]. It is therefore legitimate, for example, to set
* ref[3] to 999.0 to choose the more northerly solution - the default
* if the LATPOLE card is omitted from the FITS header. For the
* special case where the reference point of the projection is at
* native latitude zero, its celestial latitude is zero, and
* LONGPOLE = +/- 90 then the native latitude of the pole is not
* determined by the first three reference values and LATPOLE
* specifies it completely.
*
* The remaining members of the celprm struct are maintained by the
* initialization routines and should not be modified. This is done for the
* sake of efficiency and to allow an arbitrary number of contexts to be
* maintained simultaneously.
*
* double euler[5]
* Euler angles and associated intermediaries derived from the
* coordinate reference values.
* int (*prjfwd)()
* int (*prjrev)()
* Pointers to the forward and reverse projection routines.
*
*
* WCS projection codes
* --------------------
* Zenithals/azimuthals:
* AZP: zenithal/azimuthal perspective
* TAN: gnomonic
* SIN: synthesis (generalized orthographic)
* STG: stereographic
* ARC: zenithal/azimuthal equidistant
* ZPN: zenithal/azimuthal polynomial
* ZEA: zenithal/azimuthal equal area
* AIR: Airy
* TNX: IRAF's polynomial correction to TAN
*
* Cylindricals:
* CYP: cylindrical perspective
* CAR: Cartesian
* MER: Mercator
* CEA: cylindrical equal area
*
* Conics:
* COP: conic perspective
* COD: conic equidistant
* COE: conic equal area
* COO: conic orthomorphic
*
* Polyconics:
* BON: Bonne
* PCO: polyconic
*
* Pseudo-cylindricals:
* GLS: Sanson-Flamsteed (global sinusoidal)
* PAR: parabolic
* MOL: Mollweide
*
* Conventional:
* AIT: Hammer-Aitoff
*
* Quad-cubes:
* CSC: COBE quadrilateralized spherical cube
* QSC: quadrilateralized spherical cube
* TSC: tangential spherical cube
*
* Author: Mark Calabretta, Australia Telescope National Facility
* IRAF's TNX added by E.Bertin 2000/03/28
* Filtering of abs(phi)>180 and abs(theta)>90 added by E.Bertin 2000/11/11
* $Id: cel.c,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include <string.h>
#include "wcstrig.h"
#include "cel.h"
#include "sph.h"
#include "tnx.h"
int npcode = 26;
char pcodes[26][4] =
{"AZP", "TAN", "SIN", "STG", "ARC", "ZPN", "ZEA", "AIR", "CYP", "CAR",
"MER", "CEA", "COP", "COD", "COE", "COO", "BON", "PCO", "GLS", "PAR",
"AIT", "MOL", "CSC", "QSC", "TSC", "TNX"};
/* Map error number to error message for each function. */
const char *celset_errmsg[] = {
0,
"Invalid coordinate transformation parameters",
"Ill-conditioned coordinate transformation parameters"};
const char *celfwd_errmsg[] = {
0,
"Invalid coordinate transformation parameters",
"Invalid projection parameters",
"Invalid value of (lng,lat)"};
const char *celrev_errmsg[] = {
0,
"Invalid coordinate transformation parameters",
"Invalid projection parameters",
"Invalid value of (x,y)"};
int celset(pcode, cel, prj)
const char pcode[4];
struct celprm *cel;
struct prjprm *prj;
{
int dophip;
const double tol = 1.0e-10;
double clat0, cphip, cthe0, theta0, slat0, sphip, sthe0;
double latp, latp1, latp2;
double u, v, x, y, z;
/* Set pointers to the forward and reverse projection routines. */
if (strcmp(pcode, "AZP") == 0) {
cel->prjfwd = azpfwd;
cel->prjrev = azprev;
theta0 = 90.0;
} else if (strcmp(pcode, "TAN") == 0) {
cel->prjfwd = tanfwd;
cel->prjrev = tanrev;
theta0 = 90.0;
} else if (strcmp(pcode, "SIN") == 0) {
cel->prjfwd = sinfwd;
cel->prjrev = sinrev;
theta0 = 90.0;
} else if (strcmp(pcode, "STG") == 0) {
cel->prjfwd = stgfwd;
cel->prjrev = stgrev;
theta0 = 90.0;
} else if (strcmp(pcode, "ARC") == 0) {
cel->prjfwd = arcfwd;
cel->prjrev = arcrev;
theta0 = 90.0;
} else if (strcmp(pcode, "ZPN") == 0) {
cel->prjfwd = zpnfwd;
cel->prjrev = zpnrev;
theta0 = 90.0;
} else if (strcmp(pcode, "ZEA") == 0) {
cel->prjfwd = zeafwd;
cel->prjrev = zearev;
theta0 = 90.0;
} else if (strcmp(pcode, "AIR") == 0) {
cel->prjfwd = airfwd;
cel->prjrev = airrev;
theta0 = 90.0;
} else if (strcmp(pcode, "CYP") == 0) {
cel->prjfwd = cypfwd;
cel->prjrev = cyprev;
theta0 = 0.0;
} else if (strcmp(pcode, "CAR") == 0) {
cel->prjfwd = carfwd;
cel->prjrev = carrev;
theta0 = 0.0;
} else if (strcmp(pcode, "MER") == 0) {
cel->prjfwd = merfwd;
cel->prjrev = merrev;
theta0 = 0.0;
} else if (strcmp(pcode, "CEA") == 0) {
cel->prjfwd = ceafwd;
cel->prjrev = cearev;
theta0 = 0.0;
} else if (strcmp(pcode, "COP") == 0) {
cel->prjfwd = copfwd;
cel->prjrev = coprev;
theta0 = prj->p[1];
} else if (strcmp(pcode, "COD") == 0) {
cel->prjfwd = codfwd;
cel->prjrev = codrev;
theta0 = prj->p[1];
} else if (strcmp(pcode, "COE") == 0) {
cel->prjfwd = coefwd;
cel->prjrev = coerev;
theta0 = prj->p[1];
} else if (strcmp(pcode, "COO") == 0) {
cel->prjfwd = coofwd;
cel->prjrev = coorev;
theta0 = prj->p[1];
} else if (strcmp(pcode, "BON") == 0) {
cel->prjfwd = bonfwd;
cel->prjrev = bonrev;
theta0 = 0.0;
} else if (strcmp(pcode, "PCO") == 0) {
cel->prjfwd = pcofwd;
cel->prjrev = pcorev;
theta0 = 0.0;
} else if (strcmp(pcode, "GLS") == 0) {
cel->prjfwd = glsfwd;
cel->prjrev = glsrev;
theta0 = 0.0;
} else if (strcmp(pcode, "PAR") == 0) {
cel->prjfwd = parfwd;
cel->prjrev = parrev;
theta0 = 0.0;
} else if (strcmp(pcode, "AIT") == 0) {
cel->prjfwd = aitfwd;
cel->prjrev = aitrev;
theta0 = 0.0;
} else if (strcmp(pcode, "MOL") == 0) {
cel->prjfwd = molfwd;
cel->prjrev = molrev;
theta0 = 0.0;
} else if (strcmp(pcode, "CSC") == 0) {
cel->prjfwd = cscfwd;
cel->prjrev = cscrev;
theta0 = 0.0;
} else if (strcmp(pcode, "QSC") == 0) {
cel->prjfwd = qscfwd;
cel->prjrev = qscrev;
theta0 = 0.0;
} else if (strcmp(pcode, "TSC") == 0) {
cel->prjfwd = tscfwd;
cel->prjrev = tscrev;
theta0 = 0.0;
} else if (strcmp(pcode, "TNX") == 0) {
cel->prjfwd = tnxfwd;
cel->prjrev = tnxrev;
theta0 = 90.0;
} else {
/* Unrecognized projection code. */
return 1;
}
/* Set default for native longitude of the celestial pole? */
dophip = (cel->ref[2] == 999.0);
/* Compute celestial coordinates of the native pole. */
if (theta0 == 90.0) {
/* Reference point is at the native pole. */
if (dophip) {
/* Set default for longitude of the celestial pole. */
cel->ref[2] = 180.0;
}
latp = cel->ref[1];
cel->ref[3] = latp;
cel->euler[0] = cel->ref[0];
cel->euler[1] = 90.0 - latp;
} else {
/* Reference point away from the native pole. */
/* Set default for longitude of the celestial pole. */
if (dophip) {
cel->ref[2] = (cel->ref[1] < theta0) ? 180.0 : 0.0;
}
clat0 = wcs_cosd(cel->ref[1]);
slat0 = wcs_sind(cel->ref[1]);
cphip = wcs_cosd(cel->ref[2]);
sphip = wcs_sind(cel->ref[2]);
cthe0 = wcs_cosd(theta0);
sthe0 = wcs_sind(theta0);
x = cthe0*cphip;
y = sthe0;
z = sqrt(x*x + y*y);
if (z == 0.0) {
if (slat0 != 0.0) {
return 1;
}
/* latp determined by LATPOLE in this case. */
latp = cel->ref[3];
} else {
if (fabs(slat0/z) > 1.0) {
return 1;
}
u = wcs_atan2d(y,x);
v = wcs_acosd(slat0/z);
latp1 = u + v;
if (latp1 > 180.0) {
latp1 -= 360.0;
} else if (latp1 < -180.0) {
latp1 += 360.0;
}
latp2 = u - v;
if (latp2 > 180.0) {
latp2 -= 360.0;
} else if (latp2 < -180.0) {
latp2 += 360.0;
}
if (fabs(cel->ref[3]-latp1) < fabs(cel->ref[3]-latp2)) {
if (fabs(latp1) < 90.0+tol) {
latp = latp1;
} else {
latp = latp2;
}
} else {
if (fabs(latp2) < 90.0+tol) {
latp = latp2;
} else {
latp = latp1;
}
}
cel->ref[3] = latp;
}
cel->euler[1] = 90.0 - latp;
z = wcs_cosd(latp)*clat0;
if (fabs(z) < tol) {
if (fabs(clat0) < tol) {
/* Celestial pole at the reference point. */
cel->euler[0] = cel->ref[0];
cel->euler[1] = 90.0 - theta0;
} else if (latp > 0.0) {
/* Celestial pole at the native north pole.*/
cel->euler[0] = cel->ref[0] + cel->ref[2] - 180.0;
cel->euler[1] = 0.0;
} else if (latp < 0.0) {
/* Celestial pole at the native south pole. */
cel->euler[0] = cel->ref[0] - cel->ref[2];
cel->euler[1] = 180.0;
}
} else {
x = (sthe0 - wcs_sind(latp)*slat0)/z;
y = sphip*cthe0/clat0;
if (x == 0.0 && y == 0.0) {
return 1;
}
cel->euler[0] = cel->ref[0] - wcs_atan2d(y,x);
}
/* Make euler[0] the same sign as ref[0]. */
if (cel->ref[0] >= 0.0) {
if (cel->euler[0] < 0.0) cel->euler[0] += 360.0;
} else {
if (cel->euler[0] > 0.0) cel->euler[0] -= 360.0;
}
}
cel->euler[2] = cel->ref[2];
cel->euler[3] = wcs_cosd(cel->euler[1]);
cel->euler[4] = wcs_sind(cel->euler[1]);
cel->flag = CELSET;
/* Check for ill-conditioned parameters. */
if (fabs(latp) > 90.0+tol) {
return 2;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int celfwd(pcode, lng, lat, cel, phi, theta, prj, x, y)
const char pcode[4];
const double lng, lat;
struct celprm *cel;
double *phi, *theta;
struct prjprm *prj;
double *x, *y;
{
int err;
if (cel->flag != CELSET) {
if (celset(pcode, cel, prj)) return 1;
}
/* Compute native coordinates. */
sphfwd(lng, lat, cel->euler, phi, theta);
/* Apply forward projection. */
if ((err = cel->prjfwd(*phi, *theta, prj, x, y))) {
return err == 1 ? 2 : 3;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int celrev(pcode, x, y, prj, phi, theta, cel, lng, lat)
const char pcode[4];
const double x, y;
struct prjprm *prj;
double *phi, *theta;
struct celprm *cel;
double *lng, *lat;
{
int err;
if (cel->flag != CELSET) {
if(celset(pcode, cel, prj)) return 1;
}
/* Apply reverse projection. */
if ((err = cel->prjrev(x, y, prj, phi, theta))) {
return err == 1 ? 2 : 3;
}
if (fabs(*phi)>180.0 || fabs(*theta)>90.0)
return 2;
/* Compute native coordinates. */
sphrev(*phi, *theta, cel->euler, lng, lat);
return 0;
}
src/wcs/cel.h 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: cel.h,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifndef WCSLIB_CEL
#define WCSLIB_CEL
#include "proj.h"
#ifdef __cplusplus
extern "C" {
#endif
extern int npcode;
extern char pcodes[26][4];
struct celprm {
int flag;
double ref[4];
double euler[5];
#if __STDC__ || defined(__cplusplus)
int (*prjfwd)(const double, const double,
struct prjprm *,
double *, double *);
int (*prjrev)(const double, const double,
struct prjprm *,
double *, double *);
#else
int (*prjfwd)();
int (*prjrev)();
#endif
};
#if __STDC__ || defined(__cplusplus)
int celset(const char *, struct celprm *, struct prjprm *);
int celfwd(const char *,
const double, const double,
struct celprm *,
double *, double *,
struct prjprm *,
double *, double *);
int celrev(const char *,
const double, const double,
struct prjprm *,
double *, double *,
struct celprm *,
double *, double *);
#else
int celset(), celfwd(), celrev();
#endif
extern const char *celset_errmsg[];
extern const char *celfwd_errmsg[];
extern const char *celrev_errmsg[];
#define CELSET 137
#ifdef __cplusplus
}
#endif
#endif /* WCSLIB_CEL */
src/wcs/lin.c 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
* -------------------
* These utility routines apply the linear transformation defined by the WCS
* FITS header cards. There are separate routines for the image-to-pixel,
* linfwd(), and pixel-to-image, linrev(), transformations.
*
* An initialization routine, linset(), computes intermediate values from
* the transformation parameters but need not be called explicitly - see the
* explanation of lin.flag below.
*
* An auxiliary matrix inversion routine, matinv(), is included. It uses
* LU-triangular factorization with scaled partial pivoting.
*
*
* Initialization routine; linset()
* --------------------------------
* Initializes members of a linprm data structure which hold intermediate
* values. Note that this routine need not be called directly; it will be
* invoked by linfwd() and linrev() if the "flag" structure member is
* anything other than a predefined magic value.
*
* Given and/or returned:
* lin linprm* Linear transformation parameters (see below).
*
* Function return value:
* int Error status
* 0: Success.
* 1: Memory allocation error.
* 2: PC matrix is singular.
*
* Forward transformation; linfwd()
* --------------------------------
* Compute pixel coordinates from image coordinates. Note that where
* celestial coordinate systems are concerned the image coordinates
* correspond to (x,y) in the plane of projection, not celestial (lng,lat).
*
* Given:
* imgcrd const double[]
* Image (world) coordinate.
*
* Given and returned:
* lin linprm* Linear transformation parameters (see below).
*
* Returned:
* pixcrd d[] Pixel coordinate.
*
* Function return value:
* int Error status
* 0: Success.
* 1: The transformation is not invertible.
*
* Reverse transformation; linrev()
* --------------------------------
* Compute image coordinates from pixel coordinates. Note that where
* celestial coordinate systems are concerned the image coordinates
* correspond to (x,y) in the plane of projection, not celestial (lng,lat).
*
* Given:
* pixcrd const double[]
* Pixel coordinate.
*
* Given and/or returned:
* lin linprm* Linear transformation parameters (see below).
*
* Returned:
* imgcrd d[] Image (world) coordinate.
*
* Function return value:
* int Error status
* 0: Success.
* 1: Error.
*
* Linear transformation parameters
* --------------------------------
* The linprm struct consists of the following:
*
* int flag
* This flag must be set to zero whenever any of the following members
* are set or modified. This signals the initialization routine,
* linset(), to recompute intermediaries.
* int naxis
* Number of image axes.
* double *crpix
* Pointer to the first element of an array of double containing the
* coordinate reference pixel, CRPIXn.
* double *pc
* Pointer to the first element of the PC (pixel coordinate)
* transformation matrix.
* double *cdelt
* Pointer to the first element of an array of double containing the
* coordinate increments, CDELTn.
*
* The remaining members of the linprm struct are maintained by the
* initialization routine and should not be modified.
*
* double *piximg
* Pointer to the first element of the matrix containing the product
* of the CDELTn diagonal matrix and the PC matrix.
* double *imgpix
* Pointer to the first element of the inverse of the piximg matrix.
*
* linset allocates storage for the above arrays using malloc(). Note,
* however, that these routines do not free this storage so if a linprm
* variable has itself been malloc'd then these structure members must be
* explicitly freed before the linprm variable is free'd otherwise a memory
* leak will result.
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: lin.c,v 1.1.1.1 2003/07/09 19:06:26 bertin Exp $
*===========================================================================*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include <stdlib.h>
#include "lin.h"
/* Map error number to error message for each function. */
const char *linset_errmsg[] = {
0,
"Memory allocation error",
"PC matrix is singular"};
const char *linfwd_errmsg[] = {
0,
"Memory allocation error",
"PC matrix is singular"};
const char *linrev_errmsg[] = {
0,
"Memory allocation error",
"PC matrix is singular"};
int linset(lin)
struct linprm *lin;
{
int i, ij, j, mem, n;
n = lin->naxis;
/* Allocate memory for internal arrays. */
mem = n * n * sizeof(double);
lin->piximg = (double*)malloc(mem);
if (lin->piximg == (double*)0) return 1;
lin->imgpix = (double*)malloc(mem);
if (lin->imgpix == (double*)0) {
free(lin->piximg);
return 1;
}
/* Compute the pixel-to-image transformation matrix. */
for (i = 0, ij = 0; i < n; i++) {
for (j = 0; j < n; j++, ij++) {
lin->piximg[ij] = lin->cdelt[i] * lin->pc[ij];
}
}
/* Compute the image-to-pixel transformation matrix. */
if (matinv(n, lin->piximg, lin->imgpix)) return 2;
lin->flag = LINSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int linfwd(imgcrd, lin, pixcrd)
const double imgcrd[];
struct linprm *lin;
double pixcrd[];
{
int i, ij, j, n;
n = lin->naxis;
if (lin->flag != LINSET) {
if (linset(lin)) return 1;
}
for (i = 0, ij = 0; i < n; i++) {
pixcrd[i] = 0.0;
for (j = 0; j < n; j++, ij++) {
pixcrd[i] += lin->imgpix[ij] * imgcrd[j];
}
}
for (j = 0; j < n; j++) {
pixcrd[j] += lin->crpix[j];
}
return 0;
}
/*--------------------------------------------------------------------------*/
int linrev(pixcrd, lin, imgcrd)
const double pixcrd[];
struct linprm *lin;
double imgcrd[];
{
int i, ij, j, n;
double temp;
n = lin->naxis;
if (lin->flag != LINSET) {
if (linset(lin)) return 1;
}
for (i = 0; i < n; i++) {
imgcrd[i] = 0.0;
}
for (j = 0; j < n; j++) {
temp = pixcrd[j] - lin->crpix[j];
for (i = 0, ij = j; i < n; i++, ij+=n) {
imgcrd[i] += lin->piximg[ij] * temp;
}
}
return 0;
}
/*--------------------------------------------------------------------------*/
int matinv(n, mat, inv)
const int n;
const double mat[];
double inv[];
{
register int i, ij, ik, j, k, kj, pj;
int itemp, mem, *mxl, *lxm, pivot;
double colmax, *lu, *rowmax, dtemp;
/* Allocate memory for internal arrays. */
mem = n * sizeof(int);
if ((mxl = (int*)malloc(mem)) == (int*)0) return 1;
if ((lxm = (int*)malloc(mem)) == (int*)0) {
free(mxl);
return 1;
}
mem = n * sizeof(double);
if ((rowmax = (double*)malloc(mem)) == (double*)0) {
free(mxl);
free(lxm);
return 1;
}
mem *= n;
if ((lu = (double*)malloc(mem)) == (double*)0) {
free(mxl);
free(lxm);
free(rowmax);
return 1;
}
/* Initialize arrays. */
for (i = 0, ij = 0; i < n; i++) {
/* Vector which records row interchanges. */
mxl[i] = i;
rowmax[i] = 0.0;
for (j = 0; j < n; j++, ij++) {
dtemp = fabs(mat[ij]);
if (dtemp > rowmax[i]) rowmax[i] = dtemp;
lu[ij] = mat[ij];
}
/* A row of zeroes indicates a singular matrix. */
if (rowmax[i] == 0.0) {
free(mxl);
free(lxm);
free(rowmax);
free(lu);
return 2;
}
}
/* Form the LU triangular factorization using scaled partial pivoting. */
for (k = 0; k < n; k++) {
/* Decide whether to pivot. */
colmax = fabs(lu[k*n+k]) / rowmax[k];
pivot = k;
for (i = k+1; i < n; i++) {
ik = i*n + k;
dtemp = fabs(lu[ik]) / rowmax[i];
if (dtemp > colmax) {
colmax = dtemp;
pivot = i;
}
}
if (pivot > k) {
/* We must pivot, interchange the rows of the design matrix. */
for (j = 0, pj = pivot*n, kj = k*n; j < n; j++, pj++, kj++) {
dtemp = lu[pj];
lu[pj] = lu[kj];
lu[kj] = dtemp;
}
/* Amend the vector of row maxima. */
dtemp = rowmax[pivot];
rowmax[pivot] = rowmax[k];
rowmax[k] = dtemp;
/* Record the interchange for later use. */
itemp = mxl[pivot];
mxl[pivot] = mxl[k];
mxl[k] = itemp;
}
/* Gaussian elimination. */
for (i = k+1; i < n; i++) {
ik = i*n + k;
/* Nothing to do if lu[ik] is zero. */
if (lu[ik] != 0.0) {
/* Save the scaling factor. */
lu[ik] /= lu[k*n+k];
/* Subtract rows. */
for (j = k+1; j < n; j++) {
lu[i*n+j] -= lu[ik]*lu[k*n+j];
}
}
}
}
/* mxl[i] records which row of mat corresponds to row i of lu. */
/* lxm[i] records which row of lu corresponds to row i of mat. */
for (i = 0; i < n; i++) {
lxm[mxl[i]] = i;
}
/* Determine the inverse matrix. */
for (i = 0, ij = 0; i < n; i++) {
for (j = 0; j < n; j++, ij++) {
inv[ij] = 0.0;
}
}
for (k = 0; k < n; k++) {
inv[lxm[k]*n+k] = 1.0;
/* Forward substitution. */
for (i = lxm[k]+1; i < n; i++) {
for (j = lxm[k]; j < i; j++) {
inv[i*n+k] -= lu[i*n+j]*inv[j*n+k];
}
}
/* Backward substitution. */
for (i = n-1; i >= 0; i--) {
for (j = i+1; j < n; j++) {
inv[i*n+k] -= lu[i*n+j]*inv[j*n+k];
}
inv[i*n+k] /= lu[i*n+i];
}
}
free(mxl);
free(lxm);
free(rowmax);
free(lu);
return 0;
}
src/wcs/lin.h 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: lin.h,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifndef WCSLIB_LIN
#define WCSLIB_LIN
#ifdef __cplusplus
extern "C" {
#endif
#if !defined(__STDC__) && !defined(__cplusplus)
#ifndef const
#define const
#endif
#endif
struct linprm {
int flag;
int naxis;
double *crpix;
double *pc;
double *cdelt;
/* Intermediates. */
double *piximg;
double *imgpix;
};
#if __STDC__ || defined(__cplusplus)
int linset(struct linprm *);
int linfwd(const double[], struct linprm *, double[]);
int linrev(const double[], struct linprm *, double[]);
int matinv(const int, const double [], double []);
#else
int linset(), linfwd(), linrev(), matinv();
#endif
extern const char *linset_errmsg[];
extern const char *linfwd_errmsg[];
extern const char *linrev_errmsg[];
#define LINSET 137
#ifdef __cplusplus
};
#endif
#endif /* WCSLIB_LIN */
src/wcs/poly.c 0 → 100644
View file @ 3b348a94
/*
poly.c
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: A program using Polynomials
*
* Author: E.BERTIN (IAP)
*
* Contents: Polynomial fitting
*
* Last modify: 08/03/2005
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "poly.h"
#define QCALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)calloc((size_t)(nel),sizeof(typ)))) \
qerror("Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QMALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)malloc((size_t)(nel)*sizeof(typ)))) \
qerror("Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
/********************************* qerror ************************************/
/*
I hope it will never be used!
*/
void qerror(char *msg1, char *msg2)
{
fprintf(stderr, "\n> %s%s\n\n",msg1,msg2);
exit(-1);
}
/****** poly_init ************************************************************
PROTO polystruct *poly_init(int *group, int ndim, int *degree, int ngroup)
PURPOSE Allocate and initialize a polynom structure.
INPUT 1D array containing the group for each parameter,
number of dimensions (parameters),
1D array with the polynomial degree for each group,
number of groups.
OUTPUT polystruct pointer.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 08/03/2003
***/
polystruct *poly_init(int *group, int ndim, int *degree, int ngroup)
{
void qerror(char *msg1, char *msg2);
polystruct *poly;
char str[512];
int nd[POLY_MAXDIM];
int *groupt,
d,g,n,num,den;
QCALLOC(poly, polystruct, 1);
if ((poly->ndim=ndim) > POLY_MAXDIM)
{
sprintf(str, "The dimensionality of the polynom (%d) exceeds the maximum\n"
"allowed one (%d)", ndim, POLY_MAXDIM);
qerror("*Error*: ", str);
}
if (ndim)
QMALLOC(poly->group, int, poly->ndim);
for (groupt=poly->group, d=ndim; d--;)
*(groupt++) = *(group++)-1;
poly->ngroup = ngroup;
if (ngroup)
{
group = poly->group; /* Forget the original *group */
QMALLOC(poly->degree, int, poly->ngroup);
/*-- Compute the number of context parameters for each group */
memset(nd, 0, ngroup*sizeof(int));
for (d=0; d<ndim; d++)
{
if ((g=group[d])>=ngroup)
qerror("*Error*: polynomial GROUP out of range", "");
nd[g]++;
}
}
/* Compute the total number of coefficients */
poly->ncoeff = 1;
for (g=0; g<ngroup; g++)
{
if ((d=poly->degree[g]=*(degree++))>POLY_MAXDEGREE)
{
sprintf(str, "The degree of the polynom (%d) exceeds the maximum\n"
"allowed one (%d)", poly->degree[g], POLY_MAXDEGREE);
qerror("*Error*: ", str);
}
/*-- There are (n+d)!/(n!d!) coeffs per group, that is Prod_(i<=d) (n+i)/i */
for (num=den=1, n=nd[g]; d; num*=(n+d), den*=d--);
poly->ncoeff *= num/den;
}
QMALLOC(poly->basis, double, poly->ncoeff);
QCALLOC(poly->coeff, double, poly->ncoeff);
return poly;
}
/****** poly_end *************************************************************
PROTO void poly_end(polystruct *poly)
PURPOSE Free a polynom structure and everything it contains.
INPUT polystruct pointer.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 09/04/2000
***/
void poly_end(polystruct *poly)
{
if (poly)
{
free(poly->coeff);
free(poly->basis);
free(poly->degree);
free(poly->group);
free(poly);
}
}
/****** poly_func ************************************************************
PROTO double poly_func(polystruct *poly, double *pos)
PURPOSE Evaluate a multidimensional polynom.
INPUT polystruct pointer,
pointer to the 1D array of input vector data.
OUTPUT Polynom value.
NOTES Values of the basis functions are updated in poly->basis.
AUTHOR E. Bertin (IAP)
VERSION 03/03/2004
***/
double poly_func(polystruct *poly, double *pos)
{
double xpol[POLY_MAXDIM+1];
double *post, *xpolt, *basis, *coeff, xval;
long double val;
int expo[POLY_MAXDIM+1], gexpo[POLY_MAXDIM+1];
int *expot, *degree,*degreet, *group,*groupt, *gexpot,
d,g,t, ndim;
/* Prepare the vectors and counters */
ndim = poly->ndim;
basis = poly->basis;
coeff = poly->coeff;
group = poly->group;
degree = poly->degree;
if (ndim)
{
for (xpolt=xpol, expot=expo, post=pos, d=ndim; --d;)
{
*(++xpolt) = 1.0;
*(++expot) = 0;
}
for (gexpot=gexpo, degreet=degree, g=poly->ngroup; g--;)
*(gexpot++) = *(degreet++);
if (gexpo[*group])
gexpo[*group]--;
}
/* The constant term is handled separately */
val = *(coeff++);
*(basis++) = 1.0;
*expo = 1;
*xpol = *pos;
/* Compute the rest of the polynom */
for (t=poly->ncoeff; --t; )
{
/*-- xpol[0] contains the current product of the x^n's */
val += (*(basis++)=*xpol)**(coeff++);
/*-- A complex recursion between terms of the polynom speeds up computations */
/*-- Not too good for roundoff errors (prefer Horner's), but much easier for */
/*-- multivariate polynomials: this is why we use a long double accumulator */
post = pos;
groupt = group;
expot = expo;
xpolt = xpol;
for (d=0; d<ndim; d++, groupt++)
if (gexpo[*groupt]--)
{
++*(expot++);
xval = (*(xpolt--) *= *post);
while (d--)
*(xpolt--) = xval;
break;
}
else
{
gexpo[*groupt] = *expot;
*(expot++) = 0;
*(xpolt++) = 1.0;
post++;
}
}
return (double)val;
}
/****** poly_fit *************************************************************
PROTO double poly_fit(polystruct *poly, double *x, double *y, double *w,
int ndata, double *extbasis)
PURPOSE Least-Square fit of a multidimensional polynom to weighted data.
INPUT polystruct pointer,
pointer to the (pseudo)2D array of inputs to basis functions,
pointer to the 1D array of data values,
pointer to the 1D array of data weights,
number of data points,
pointer to a (pseudo)2D array of computed basis function values.
OUTPUT Chi2 of the fit.
NOTES If different from NULL, extbasis can be provided to store the
values of the basis functions. If x==NULL and extbasis!=NULL, the
precomputed basis functions stored in extbasis are used (which saves
CPU). If w is NULL, all points are given identical weight.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 08/03/2005
***/
void poly_fit(polystruct *poly, double *x, double *y, double *w, int ndata,
double *extbasis)
{
void qerror(char *msg1, char *msg2);
double /*offset[POLY_MAXDIM],*/x2[POLY_MAXDIM],
*alpha,*alphat, *beta,*betat, *basis,*basis1,*basis2, *coeff,
*extbasist,*xt,
val,wval,yval;
int ncoeff, ndim, matsize,
d,i,j,n;
if (!x && !extbasis)
qerror("*Internal Error*: One of x or extbasis should be "
"different from NULL\nin ", "poly_func()");
ncoeff = poly->ncoeff;
ndim = poly->ndim;
matsize = ncoeff*ncoeff;
basis = poly->basis;
extbasist = extbasis;
QCALLOC(alpha, double, matsize);
QCALLOC(beta, double, ncoeff);
/* Subtract an average offset to maintain precision (droped for now ) */
/*
if (x)
{
for (d=0; d<ndim; d++)
offset[d] = 0.0;
xt = x;
for (n=ndata; n--;)
for (d=0; d<ndim; d++)
offset[d] += *(xt++);
for (d=0; d<ndim; d++)
offset[d] /= (double)ndata;
}
*/
/* Build the covariance matrix */
xt = x;
for (n=ndata; n--;)
{
if (x)
{
/*---- If x!=NULL, compute the basis functions */
for (d=0; d<ndim; d++)
x2[d] = *(xt++)/* - offset[d]*/;
poly_func(poly, x2);
/*---- If, in addition, extbasis is provided, then fill it */
if (extbasis)
for (basis1=basis,j=ncoeff; j--;)
*(extbasist++) = *(basis1++);
}
else
/*---- If x==NULL, then rely on pre-computed basis functions */
for (basis1=basis,j=ncoeff; j--;)
*(basis1++) = *(extbasist++);
basis1 = basis;
wval = w? *(w++) : 1.0;
yval = *(y++);
betat = beta;
alphat = alpha;
for (j=ncoeff; j--;)
{
val = *(basis1++)*wval;
*(betat++) += val*yval;
for (basis2=basis,i=ncoeff; i--;)
*(alphat++) += val**(basis2++);
}
}
/* Solve the system */
poly_solve(alpha,beta,ncoeff);
free(alpha);
/* Now fill the coeff array with the result of the fit */
betat = beta;
coeff = poly->coeff;
for (j=ncoeff; j--;)
*(coeff++) = *(betat++);
/*
poly_addcste(poly, offset);
*/
free(beta);
return;
}
/****** poly_addcste *********************************************************
PROTO void poly_addcste(polystruct *poly, double *cste)
PURPOSE Modify matrix coefficients to mimick the effect of adding a cst to
the input of a polynomial.
INPUT Pointer to the polynomial structure,
Pointer to the vector of cst.
OUTPUT -.
NOTES Requires quadruple-precision. **For the time beeing, this function
returns completely wrong results!!**
AUTHOR E. Bertin (IAP)
VERSION 03/03/2004
***/
void poly_addcste(polystruct *poly, double *cste)
{
long double *acoeff;
double *coeff,*mcoeff,*mcoefft,
val;
int *mpowers,*powers,*powerst,*powerst2,
i,j,n,p, denum, flag, maxdegree, ncoeff, ndim;
ncoeff = poly->ncoeff;
ndim = poly->ndim;
maxdegree = 0;
for (j=0; j<poly->ngroup; j++)
if (maxdegree < poly->degree[j])
maxdegree = poly->degree[j];
maxdegree++; /* Actually we need maxdegree+1 terms */
QCALLOC(acoeff, long double, ncoeff);
QCALLOC(mcoeff, double, ndim*maxdegree);
QCALLOC(mpowers, int, ndim);
mcoefft = mcoeff; /* To avoid gcc -Wall warnings */
powerst = powers = poly_powers(poly);
coeff = poly->coeff;
for (i=0; i<ncoeff; i++)
{
for (j=0; j<ndim; j++)
{
mpowers[j] = n = *(powerst++);
mcoefft = mcoeff+j*maxdegree+n;
denum = 1;
val = 1.0;
for (p=n+1; p--;)
{
*(mcoefft--) = val;
val *= (cste[j]*(n--))/(denum++); /* This is C_n^p X^(n-p) */
}
}
/*-- Update all valid coefficients */
powerst2 = powers;
for (p=0; p<ncoeff; p++)
{
/*---- Check that this combination of powers is included in the series above */
flag = 0;
for (j=0; j<ndim; j++)
if (mpowers[j] < powerst2[j])
{
flag = 1;
powerst2 += ndim;
break;
}
if (flag == 1)
continue;
val = 1.0;
mcoefft = mcoeff;
for (j=ndim; j--; mcoefft += maxdegree)
val *= mcoefft[*(powerst2++)];
acoeff[i] += val*coeff[p];
/*
printf("%g \n", val);
*/
}
}
/* Add the new coefficients to the previous ones */
for (i=0; i<ncoeff; i++)
{
/*
printf("%g %g\n", coeff[i], (double)acoeff[i]);
*/
coeff[i] = (double)acoeff[i];
}
free(acoeff);
free(mcoeff);
free(mpowers);
free(powers);
return;
}
/****** poly_solve ************************************************************
PROTO void poly_solve(double *a, double *b, int n)
PURPOSE Solve a system of linear equations, using Cholesky decomposition or
SVD (if the former fails due to hidden correlation between variables).
INPUT Pointer to the (pseudo 2D) matrix of coefficients,
pointer to the 1D column vector,
matrix size.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 21/09/2004
***/
void poly_solve(double *a, double *b, int n)
{
double *vmat,*wmat;
if (cholsolve(a,b,n))
{
QMALLOC(vmat, double, n*n);
QMALLOC(wmat, double, n);
svdsolve(a, b, n,n, vmat,wmat);
free(vmat);
free(wmat);
}
return;
}
/****** cholsolve *************************************************************
PROTO void cholsolve(double *a, double *b, int n)
PURPOSE Solve a system of linear equations, using Cholesky decomposition.
INPUT Pointer to the (pseudo 2D) matrix of coefficients,
pointer to the 1D column vector,
matrix size.
OUTPUT -1 if the matrix is not positive-definite, 0 otherwise.
NOTES Based on Numerical Recipes, 2nd ed. (Chap 2.9). The matrix of
coefficients must be symmetric and positive definite.
AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
VERSION 28/10/2003
***/
int cholsolve(double *a, double *b, int n)
{
void qerror(char *msg1, char *msg2);
double *p, *x, sum;
int i,j,k;
/* Allocate memory to store the diagonal elements */
QMALLOC(p, double, n);
/* Cholesky decomposition */
for (i=0; i<n; i++)
for (j=i; j<n; j++)
{
for (sum=a[i*n+j],k=i-1; k>=0; k--)
sum -= a[i*n+k]*a[j*n+k];
if (i==j)
{
if (sum <= 0.0)
{
free(p);
return -1;
}
p[i] = sqrt(sum);
}
else
a[j*n+i] = sum/p[i];
}
/* Solve the system */
x = b; /* Just to save memory: the solution replaces b */
for (i=0; i<n; i++)
{
for (sum=b[i],k=i-1; k>=0; k--)
sum -= a[i*n+k]*x[k];
x[i] = sum/p[i];
}
for (i=n-1; i>=0; i--)
{
for (sum=x[i],k=i+1; k<n; k++)
sum -= a[k*n+i]*x[k];
x[i] = sum/p[i];
}
free(p);
return 0;
}
/****** svdsolve *************************************************************
PROTO void svdsolve(double *a, double *b, int m, int n, double *vmat,
double *wmat)
PURPOSE General least-square fit A.x = b, based on Singular Value
Decomposition (SVD).
Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671).
INPUT Pointer to the (pseudo 2D) matrix of coefficients,
pointer to the 1D column vector (replaced by solution in output),
number of matrix rows,
number of matrix columns,
pointer to the (pseudo 2D) SVD matrix,
pointer to the diagonal (1D) matrix of singular values.
OUTPUT -.
NOTES Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671). The a
and v matrices are transposed with respect to the N.R. convention.
AUTHOR E. Bertin (IAP)
VERSION 26/12/2003
***/
void svdsolve(double *a, double *b, int m, int n, double *vmat, double *wmat)
{
#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
(maxarg1) : (maxarg2))
#define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
(ct=bt/at,at*sqrt(1.0+ct*ct)) \
: (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define TOL 1.0e-11
void qerror(char *msg1, char *msg2);
int flag,i,its,j,jj,k,l,mmi,nm, nml;
double *w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
*bp,*tmpp, *rv1,*tmp, *sol,
c,f,h,s,x,y,z,
anorm, g, scale,
at,bt,ct,maxarg1,maxarg2,
thresh, wmax;
anorm = g = scale = 0.0;
if (m < n)
qerror("*Error*: Not enough rows for solving the system ", "in svdfit()");
sol = b; /* The solution overwrites the input column matrix */
QMALLOC(rv1, double, n);
QMALLOC(tmp, double, n);
l = nm = nml = 0; /* To avoid gcc -Wall warnings */
for (i=0;i<n;i++)
{
l = i+1;
nml = n-l;
rv1[i] = scale*g;
g = s = scale = 0.0;
if ((mmi = m - i) > 0)
{
ap = ap0 = a+i*(m+1);
for (k=mmi;k--;)
scale += fabs(*(ap++));
if (scale)
{
for (ap=ap0,k=mmi; k--; ap++)
{
*ap /= scale;
s += *ap**ap;
}
f = *ap0;
g = -SIGN(sqrt(s),f);
h = f*g-s;
*ap0 = f-g;
ap10 = a+l*m+i;
for (j=nml;j--; ap10+=m)
{
for (s=0.0,ap=ap0,ap1=ap10,k=mmi; k--;)
s += *(ap1++)**(ap++);
f = s/h;
for (ap=ap0,ap1=ap10,k=mmi; k--;)
*(ap1++) += f**(ap++);
}
for (ap=ap0,k=mmi; k--;)
*(ap++) *= scale;
}
}
wmat[i] = scale*g;
g = s = scale = 0.0;
if (i < m && i+1 != n)
{
ap = ap0 = a+i+m*l;
for (k=nml;k--; ap+=m)
scale += fabs(*ap);
if (scale)
{
for (ap=ap0,k=nml;k--; ap+=m)
{
*ap /= scale;
s += *ap**ap;
}
f=*ap0;
g = -SIGN(sqrt(s),f);
h=f*g-s;
*ap0=f-g;
rv1p = rv1+l;
for (ap=ap0,k=nml;k--; ap+=m)
*(rv1p++) = *ap/h;
ap10 = a+l+m*l;
for (j=m-l; j--; ap10++)
{
for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
s += *ap1**ap;
rv1p = rv1+l;
for (ap1=ap10,k=nml;k--; ap1+=m)
*ap1 += s**(rv1p++);
}
for (ap=ap0,k=nml;k--; ap+=m)
*ap *= scale;
}
}
anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
}
for (i=n-1;i>=0;i--)
{
if (i < n-1)
{
if (g)
{
ap0 = a+l*m+i;
vp0 = vmat+i*n+l;
vp10 = vmat+l*n+l;
g *= *ap0;
for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
*(vp++) = *ap/g;
for (j=nml; j--; vp10+=n)
{
for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
s += *ap**(vp1++);
for (vp=vp0,vp1=vp10,k=nml; k--;)
*(vp1++) += s**(vp++);
}
}
vp = vmat+l*n+i;
vp1 = vmat+i*n+l;
for (j=nml; j--; vp+=n)
*vp = *(vp1++) = 0.0;
}
vmat[i*n+i]=1.0;
g=rv1[i];
l=i;
nml = n-l;
}
for (i=(m<n?m:n); --i>=0;)
{
l=i+1;
nml = n-l;
mmi=m-i;
g=wmat[i];
ap0 = a+i*m+i;
ap10 = ap0 + m;
for (ap=ap10,j=nml;j--;ap+=m)
*ap=0.0;
if (g)
{
g=1.0/g;
for (j=nml;j--; ap10+=m)
{
for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
s += *(++ap)**(++ap1);
f = (s/(*ap0))*g;
for (ap=ap0,ap1=ap10,k=mmi;k--;)
*(ap1++) += f**(ap++);
}
for (ap=ap0,j=mmi;j--;)
*(ap++) *= g;
}
else
for (ap=ap0,j=mmi;j--;)
*(ap++)=0.0;
++(*ap0);
}
for (k=n; --k>=0;)
{
for (its=0;its<100;its++)
{
flag=1;
for (l=k;l>=0;l--)
{
nm=l-1;
if (fabs(rv1[l])+anorm == anorm)
{
flag=0;
break;
}
if (fabs(wmat[nm])+anorm == anorm)
break;
}
if (flag)
{
c=0.0;
s=1.0;
ap0 = a+nm*m;
ap10 = a+l*m;
for (i=l; i<=k; i++,ap10+=m)
{
f=s*rv1[i];
if (fabs(f)+anorm == anorm)
break;
g=wmat[i];
h=PYTHAG(f,g);
wmat[i]=h;
h=1.0/h;
c=g*h;
s=(-f*h);
for (ap=ap0,ap1=ap10,j=m; j--;)
{
z = *ap1;
y = *ap;
*(ap++) = y*c+z*s;
*(ap1++) = z*c-y*s;
}
}
}
z=wmat[k];
if (l == k)
{
if (z < 0.0)
{
wmat[k] = -z;
vp = vmat+k*n;
for (j=n; j--; vp++)
*vp = (-*vp);
}
break;
}
if (its == 99)
qerror("*Error*: No convergence in 100 SVD iterations ",
"in svdfit()");
x=wmat[l];
nm=k-1;
y=wmat[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=PYTHAG(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
c=s=1.0;
ap10 = a+l*m;
vp10 = vmat+l*n;
for (j=l;j<=nm;j++,ap10+=m,vp10+=n)
{
i=j+1;
g=rv1[i];
y=wmat[i];
h=s*g;
g=c*g;
z=PYTHAG(f,h);
rv1[j]=z;
c=f/z;
s=h/z;
f=x*c+g*s;
g=g*c-x*s;
h=y*s;
y=y*c;
for (vp=(vp1=vp10)+n,jj=n; jj--;)
{
z = *vp;
x = *vp1;
*(vp1++) = x*c+z*s;
*(vp++) = z*c-x*s;
}
z=PYTHAG(f,h);
wmat[j]=z;
if (z)
{
z=1.0/z;
c=f*z;
s=h*z;
}
f=c*g+s*y;
x=c*y-s*g;
for (ap=(ap1=ap10)+m,jj=m; jj--;)
{
z = *ap;
y = *ap1;
*(ap1++) = y*c+z*s;
*(ap++) = z*c-y*s;
}
}
rv1[l]=0.0;
rv1[k]=f;
wmat[k]=x;
}
}
wmax=0.0;
w = wmat;
for (j=n;j--; w++)
if (*w > wmax)
wmax=*w;
thresh=TOL*wmax;
w = wmat;
for (j=n;j--; w++)
if (*w < thresh)
*w = 0.0;
w = wmat;
ap = a;
tmpp = tmp;
for (j=n; j--; w++)
{
s=0.0;
if (*w)
{
bp = b;
for (i=m; i--;)
s += *(ap++)**(bp++);
s /= *w;
}
else
ap += m;
*(tmpp++) = s;
}
vp0 = vmat;
for (j=0; j<n; j++,vp0++)
{
s=0.0;
tmpp = tmp;
for (vp=vp0,jj=n; jj--; vp+=n)
s += *vp**(tmpp++);
sol[j]=s;
}
/* Free temporary arrays */
free(tmp);
free(rv1);
return;
}
#undef SIGN
#undef MAX
#undef PYTHAG
#undef TOL
/****** poly_powers ***********************************************************
PROTO int *poly_powers(polystruct *poly)
PURPOSE Return an array of powers of polynom terms
INPUT polystruct pointer,
OUTPUT Pointer to an array of polynom powers (int *), (ncoeff*ndim numbers).
NOTES The returned pointer is mallocated.
AUTHOR E. Bertin (IAP)
VERSION 23/10/2003
***/
int *poly_powers(polystruct *poly)
{
int expo[POLY_MAXDIM+1], gexpo[POLY_MAXDIM+1];
int *expot, *degree,*degreet, *group,*groupt, *gexpot,
*powers, *powerst,
d,g,t, ndim;
/* Prepare the vectors and counters */
ndim = poly->ndim;
group = poly->group;
degree = poly->degree;
QMALLOC(powers, int, ndim*poly->ncoeff);
if (ndim)
{
for (expot=expo, d=ndim; --d;)
*(++expot) = 0;
for (gexpot=gexpo, degreet=degree, g=poly->ngroup; g--;)
*(gexpot++) = *(degreet++);
if (gexpo[*group])
gexpo[*group]--;
}
/* The constant term is handled separately */
powerst = powers;
for (d=0; d<ndim; d++)
*(powerst++) = 0;
*expo = 1;
/* Compute the rest of the polynom */
for (t=poly->ncoeff; --t; )
{
for (d=0; d<ndim; d++)
*(powerst++) = expo[d];
/*-- A complex recursion between terms of the polynom speeds up computations */
groupt = group;
expot = expo;
for (d=0; d<ndim; d++, groupt++)
if (gexpo[*groupt]--)
{
++*(expot++);
break;
}
else
{
gexpo[*groupt] = *expot;
*(expot++) = 0;
}
}
return powers;
}
src/wcs/poly.h 0 → 100644
View file @ 3b348a94
/*
poly.h
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: A program using polynomial fits
*
* Author: E.BERTIN (IAP)
*
* Contents: Include for poly.c
*
* Last modify: 03/03/2004
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifndef _POLY_H_
#define _POLY_H_
/*--------------------------------- constants -------------------------------*/
#define POLY_MAXDIM 4 /* Max dimensionality of polynom */
#define POLY_MAXDEGREE 10 /* Max degree of the polynom */
/*---------------------------------- macros ---------------------------------*/
/*--------------------------- structure definitions -------------------------*/
typedef struct poly
{
double *basis; /* Current values of the basis functions */
double *coeff; /* Polynom coefficients */
int ncoeff; /* Number of coefficients */
int *group; /* Groups */
int ndim; /* dimensionality of the polynom */
int *degree; /* Degree in each group */
int ngroup; /* Number of different groups */
} polystruct;
/*---------------------------------- protos --------------------------------*/
extern polystruct *poly_init(int *group,int ndim,int *degree,int ngroup);
extern double poly_func(polystruct *poly, double *pos);
extern int cholsolve(double *a, double *b, int n),
*poly_powers(polystruct *poly);
extern void poly_addcste(polystruct *poly, double *cste),
poly_end(polystruct *poly),
poly_fit(polystruct *poly, double *x, double *y,
double *w, int ndata, double *extbasis),
poly_solve(double *a, double *b, int n),
svdsolve(double *a, double *b, int m, int n,
double *vmat, double *wmat);
#endif
src/wcs/proj.c 0 → 100644
View file @ 3b348a94
/*============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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 implementation of the spherical map projections recognized by the FITS
* "World Coordinate System" (WCS) convention.
*
* Summary of routines
* -------------------
* Each projection is implemented via separate functions for the forward,
* *fwd(), and reverse, *rev(), transformation.
*
* Initialization routines, *set(), compute intermediate values from the
* projection parameters but need not be called explicitly - see the
* explanation of prj.flag below.
*
* azpset azpfwd azprev AZP: zenithal/azimuthal perspective
* tanset tanfwd tanrev TAN: gnomonic
* sinset sinfwd sinrev SIN: orthographic/synthesis
* stgset stgfwd stgrev STG: stereographic
* arcset arcfwd arcrev ARC: zenithal/azimuthal equidistant
* zpnset zpnfwd zpnrev ZPN: zenithal/azimuthal polynomial
* zeaset zeafwd zearev ZEA: zenithal/azimuthal equal area
* airset airfwd airrev AIR: Airy
* cypset cypfwd cyprev CYP: cylindrical perspective
* carset carfwd carrev CAR: Cartesian
* merset merfwd merrev MER: Mercator
* ceaset ceafwd cearev CEA: cylindrical equal area
* copset copfwd coprev COP: conic perspective
* codset codfwd codrev COD: conic equidistant
* coeset coefwd coerev COE: conic equal area
* cooset coofwd coorev COO: conic orthomorphic
* bonset bonfwd bonrev BON: Bonne
* pcoset pcofwd pcorev PCO: polyconic
* glsset glsfwd glsrev GLS: Sanson-Flamsteed (global sinusoidal)
* parset parfwd parrev PAR: parabolic
* aitset aitfwd aitrev AIT: Hammer-Aitoff
* molset molfwd molrev MOL: Mollweide
* cscset cscfwd cscrev CSC: COBE quadrilateralized spherical cube
* qscset qscfwd qscrev QSC: quadrilateralized spherical cube
* tscset tscfwd tscrev TSC: tangential spherical cube
* tnxset tnxfwd tnxrev TNX: IRAF's gnomonic
*
*
* Initialization routine; *set()
* ------------------------------
* Initializes members of a prjprm data structure which hold intermediate
* values. Note that this routine need not be called directly; it will be
* invoked by prjfwd() and prjrev() if the "flag" structure member is
* anything other than a predefined magic value.
*
* Given and/or returned:
* prj prjprm* Projection parameters (see below).
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid projection parameters.
*
* Forward transformation; *fwd()
* -----------------------------
* Compute (x,y) coordinates in the plane of projection from native spherical
* coordinates (phi,theta).
*
* Given:
* phi, const double
* theta Longitude and latitude of the projected point in
* native spherical coordinates, in degrees.
*
* Given and returned:
* prj prjprm* Projection parameters (see below).
*
* Returned:
* x,y double* Projected coordinates.
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid projection parameters.
* 2: Invalid value of (phi,theta).
*
* Reverse transformation; *rev()
* -----------------------------
* Compute native spherical coordinates (phi,theta) from (x,y) coordinates in
* the plane of projection.
*
* Given:
* x,y const double
* Projected coordinates.
*
* Given and returned:
* prj prjprm* Projection parameters (see below).
*
* Returned:
* phi, double* Longitude and latitude of the projected point in
* theta native spherical coordinates, in degrees.
*
* Function return value:
* int Error status
* 0: Success.
* 1: Invalid projection parameters.
* 2: Invalid value of (x,y).
*
* Projection parameters
* ---------------------
* The prjprm struct consists of the following:
*
* int flag
* This flag must be set to zero whenever any of p[10] or r0 are set
* or changed. This signals the initialization routine to recompute
* intermediaries. flag may also be set to -1 to disable strict bounds
* checking for the AZP, TAN, SIN, ZPN, and COP projections.
* double r0
* r0; The radius of the generating sphere for the projection, a linear
* scaling parameter. If this is zero, it will be reset to the default
* value of 180/pi (the value for FITS WCS).
* double p[10]
* The first 10 elements contain projection parameters which correspond
* to the PROJPn keywords in FITS, so p[0] is PROJP0, and p[9] is
* PROJP9. Many projections use p[1] (PROJP1) and some also use p[2]
* (PROJP2). ZPN is the only projection which uses any of the others.
*
* The remaining members of the prjprm struct are maintained by the
* initialization routines and should not be modified. This is done for the
* sake of efficiency and to allow an arbitrary number of contexts to be
* maintained simultaneously.
*
* int n
* double w[10]
* Intermediate values derived from the projection parameters.
*
* Usage of the p[] array as it applies to each projection is described in
* the prologue to each trio of projection routines.
*
* Argument checking
* -----------------
* Forward routines:
*
* The values of phi and theta (the native longitude and latitude)
* normally lie in the range [-180,180] for phi, and [-90,90] for theta.
* However, all forward projections will accept any value of phi and will
* not normalize it.
*
* The forward projection routines do not explicitly check that theta lies
* within the range [-90,90]. They do check for any value of theta which
* produces an invalid argument to the projection equations (e.g. leading
* to division by zero). The forward routines for AZP, TAN, SIN, ZPN, and
* COP also return error 2 if (phi,theta) corresponds to the overlapped
* (far) side of the projection but also return the corresponding value of
* (x,y). This strict bounds checking may be relaxed by setting prj->flag
* to -1 (rather than 0) when these projections are initialized.
*
* Reverse routines:
*
* Error checking on the projected coordinates (x,y) is limited to that
* required to ascertain whether a solution exists. Where a solution does
* exist no check is made that the value of phi and theta obtained lie
* within the ranges [-180,180] for phi, and [-90,90] for theta.
*
* Accuracy
* --------
* Closure to a precision of at least 1E-10 degree of longitude and latitude
* has been verified for typical projection parameters on the 1 degree grid
* of native longitude and latitude (to within 5 degrees of any latitude
* where the projection may diverge).
*
* Author: Mark Calabretta, Australia Telescope National Facility
* IRAF's TNX added by E.Bertin 2000/08/23
* $Id: proj.c,v 1.1.1.1 2004/01/04 21:33:26 bertin Exp $
*===========================================================================*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include "poly.h"
#include "proj.h"
#include "tnx.h"
#include "wcsmath.h"
#include "wcstrig.h"
/* Map error number to error message for each function. */
const char *prjset_errmsg[] = {
0,
"Invalid projection parameters"};
const char *prjfwd_errmsg[] = {
0,
"Invalid projection parameters",
"Invalid value of (phi,theta)"};
const char *prjrev_errmsg[] = {
0,
"Invalid projection parameters",
"Invalid value of (x,y)"};
#define wcs_copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
/*============================================================================
* AZP: zenithal/azimuthal perspective projection.
*
* Given:
* prj->p[1] AZP distance parameter, mu in units of r0.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(mu+1)
* prj->w[1] 1/prj->w[0]
* prj->w[2] Boundary parameter, -mu for |mu| <= 1,
* -1/mu for |mu| >= 1.
*===========================================================================*/
int azpset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) prj->r0 = R2D;
prj->w[0] = prj->r0*(prj->p[1] + 1.0);
if (prj->w[0] == 0.0) {
return 1;
}
prj->w[1] = 1.0/prj->w[0];
if (fabs(prj->p[1]) <= 1.0) {
prj->w[2] = -prj->p[1];
} else {
prj->w[2] = -1.0/prj->p[1];
}
if (prj->flag == -1) {
prj->flag = -PRJSET;
} else {
prj->flag = PRJSET;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int azpfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double r, s, sthe;
if (abs(prj->flag) != PRJSET) {
if (azpset(prj)) return 1;
}
sthe = wcs_sind(theta);
s = prj->p[1] + sthe;
if (s == 0.0) {
return 2;
}
r = prj->w[0]*wcs_cosd(theta)/s;
*x = r*wcs_sind(phi);
*y = -r*wcs_cosd(phi);
if (prj->flag == PRJSET && sthe < prj->w[2]) {
return 2;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int azprev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double r, rho, s;
const double tol = 1.0e-13;
if (abs(prj->flag) != PRJSET) {
if (azpset(prj)) return 1;
}
r = sqrt(x*x + y*y);
if (r == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(x, -y);
}
rho = r*prj->w[1];
s = rho*prj->p[1]/sqrt(rho*rho+1.0);
if (fabs(s) > 1.0) {
if (fabs(s) > 1.0+tol) {
return 2;
}
*theta = wcs_atan2d(1.0,rho) - wcs_copysign(90.0,s);
} else {
*theta = wcs_atan2d(1.0,rho) - wcs_asind(s);
}
return 0;
}
/*============================================================================
* TAN: gnomonic projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
*===========================================================================*/
int tanset(prj)
struct prjprm *prj;
{
int k;
if (prj->r0 == 0.0) prj->r0 = R2D;
if (prj->flag == -1) {
prj->flag = -PRJSET;
} else {
prj->flag = PRJSET;
}
for (k = 99; k >= 0 && prj->p[k] == 0.0 && prj->p[k+100] == 0.0; k--);
if (k < 0)
k = 0;
prj->n = k;
return 0;
}
/*--------------------------------------------------------------------------*/
int tanfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double r, s, xp[2];
if (abs(prj->flag) != PRJSET) {
if(tanset(prj)) return 1;
}
s = wcs_sind(theta);
if (s == 0.0) return 2;
r = prj->r0*wcs_cosd(theta)/s;
xp[0] = r*wcs_sind(phi);
xp[1] = -r*wcs_cosd(phi);
*x = prj->inv_x? poly_func(prj->inv_x, xp) : xp[0];
*y = prj->inv_y? poly_func(prj->inv_y, xp) : xp[1];
if (prj->flag == PRJSET && s < 0.0) {
return 2;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int tanrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double xp,yp, rp;
if (abs(prj->flag) != PRJSET) {
if (tanset(prj)) return 1;
}
if (prj->n)
raw_to_pv(prj, x,y, &xp, &yp);
else
{
xp = x;
yp = y;
}
rp = sqrt(xp*xp+yp*yp);
if (rp == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(xp, -yp);
}
*theta = wcs_atan2d(prj->r0, rp);
return 0;
}
/*============================================================================
* SIN: orthographic/synthesis projection.
*
* Given:
* prj->p[1:2] SIN obliqueness parameters, alpha and beta.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] 1/r0
* prj->w[1] alpha**2 + beta**2
* prj->w[2] 2*(alpha**2 + beta**2)
* prj->w[3] 2*(alpha**2 + beta**2 + 1)
* prj->w[4] alpha**2 + beta**2 - 1
*===========================================================================*/
int sinset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) prj->r0 = R2D;
prj->w[0] = 1.0/prj->r0;
prj->w[1] = prj->p[1]*prj->p[1] + prj->p[2]*prj->p[2];
prj->w[2] = 2.0*prj->w[1];
prj->w[3] = prj->w[2] + 2.0;
prj->w[4] = prj->w[1] - 1.0;
if (prj->flag == -1) {
prj->flag = -PRJSET;
} else {
prj->flag = PRJSET;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int sinfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double cphi, cthe, sphi, t, z;
if (abs(prj->flag) != PRJSET) {
if (sinset(prj)) return 1;
}
t = (90.0 - fabs(theta))*D2R;
if (t < 1.0e-5) {
if (theta > 0.0) {
z = -t*t/2.0;
} else {
z = -2.0 + t*t/2.0;
}
cthe = t;
} else {
z = wcs_sind(theta) - 1.0;
cthe = wcs_cosd(theta);
}
cphi = wcs_cosd(phi);
sphi = wcs_sind(phi);
*x = prj->r0*(cthe*sphi + prj->p[1]*z);
*y = -prj->r0*(cthe*cphi + prj->p[2]*z);
/* Validate this solution. */
if (prj->flag == PRJSET) {
if (prj->w[1] == 0.0) {
/* Orthographic projection. */
if (theta < 0.0) {
return 2;
}
} else {
/* "Synthesis" projection. */
t = wcs_atand(prj->p[1]*sphi + prj->p[2]*cphi);
if (theta < t) {
return 2;
}
}
}
return 0;
}
/*--------------------------------------------------------------------------*/
int sinrev (x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
const double tol = 1.0e-13;
double a, b, c, d, r2, sth, sth1, sth2, sxy, x0, xp, y0, yp, z;
if (abs(prj->flag) != PRJSET) {
if (sinset(prj)) return 1;
}
/* Compute intermediaries. */
x0 = x*prj->w[0];
y0 = y*prj->w[0];
r2 = x0*x0 + y0*y0;
if (prj->w[1] == 0.0) {
/* Orthographic projection. */
if (r2 != 0.0) {
*phi = wcs_atan2d(x0, -y0);
} else {
*phi = 0.0;
}
if (r2 < 0.5) {
*theta = wcs_acosd(sqrt(r2));
} else if (r2 <= 1.0) {
*theta = wcs_asind(sqrt(1.0 - r2));
} else {
return 2;
}
} else {
/* "Synthesis" projection. */
if (r2 < 1.0e-10) {
/* Use small angle formula. */
z = -r2/2.0;
*theta = 90.0 - R2D*sqrt(r2/(1.0 - x0*prj->p[1] + y0*prj->p[2]));
} else {
sxy = 2.0*(prj->p[1]*x0 - prj->p[2]*y0);
a = prj->w[3];
b = -(sxy + prj->w[2]);
c = r2 + sxy + prj->w[4];
d = b*b - 2.0*a*c;
/* Check for a solution. */
if (d < 0.0) {
return 2;
}
d = sqrt(d);
/* Choose solution closest to pole. */
sth1 = (-b + d)/a;
sth2 = (-b - d)/a;
sth = (sth1>sth2) ? sth1 : sth2;
if (sth > 1.0) {
if (sth-1.0 < tol) {
sth = 1.0;
} else {
sth = (sth1<sth2) ? sth1 : sth2;
}
}
if (sth > 1.0 || sth < -1.0) {
return 2;
}
*theta = wcs_asind(sth);
z = sth - 1.0;
}
xp = -y0 - prj->p[2]*z;
yp = x0 - prj->p[1]*z;
if (xp == 0.0 && yp == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(yp,xp);
}
}
return 0;
}
/*============================================================================
* STG: stereographic projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] 2*r0
* prj->w[1] 1/(2*r0)
*===========================================================================*/
int stgset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 360.0/PI;
prj->w[1] = PI/360.0;
} else {
prj->w[0] = 2.0*prj->r0;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int stgfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double r, s;
if (prj->flag != PRJSET) {
if (stgset(prj)) return 1;
}
s = 1.0 + wcs_sind(theta);
if (s == 0.0) {
return 2;
}
r = prj->w[0]*wcs_cosd(theta)/s;
*x = r*wcs_sind(phi);
*y = -r*wcs_cosd(phi);
return 0;
}
/*--------------------------------------------------------------------------*/
int stgrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double r;
if (prj->flag != PRJSET) {
if (stgset(prj)) return 1;
}
r = sqrt(x*x + y*y);
if (r == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(x, -y);
}
*theta = 90.0 - 2.0*wcs_atand(r*prj->w[1]);
return 0;
}
/*============================================================================
* ARC: zenithal/azimuthal equidistant projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/180)
* prj->w[1] (180/pi)/r0
*===========================================================================*/
int arcset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 1.0;
prj->w[1] = 1.0;
} else {
prj->w[0] = prj->r0*D2R;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int arcfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double r;
if (prj->flag != PRJSET) {
if (arcset(prj)) return 1;
}
r = prj->w[0]*(90.0 - theta);
*x = r*wcs_sind(phi);
*y = -r*wcs_cosd(phi);
return 0;
}
/*--------------------------------------------------------------------------*/
int arcrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double r;
if (prj->flag != PRJSET) {
if (arcset(prj)) return 1;
}
r = sqrt(x*x + y*y);
if (r == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(x, -y);
}
*theta = 90.0 - r*prj->w[1];
return 0;
}
/*============================================================================
* ZPN: zenithal/azimuthal polynomial projection.
*
* Given:
* prj->p[0:99] Polynomial coefficients.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->n Degree of the polynomial, N.
* prj->w[0] Co-latitude of the first point of inflection (N > 2).
* prj->w[1] Radius of the first point of inflection (N > 2).
*===========================================================================*/
int zpnset(prj)
struct prjprm *prj;
{
int i, j, k;
double d, d1, d2, r, zd, zd1, zd2;
const double tol = 1.0e-13;
if (prj->r0 == 0.0) prj->r0 = R2D;
/* Find the highest non-zero coefficient. */
for (k = 99; k >= 0 && prj->p[k] == 0.0; k--);
if (k < 0) return 1;
prj->n = k;
if (k >= 3) {
/* Find the point of inflection closest to the pole. */
zd1 = 0.0;
d1 = prj->p[1];
if (d1 <= 0.0) {
return 1;
}
/* Find the point where the derivative first goes negative. */
for (i = 0; i < 180; i++) {
zd2 = i*D2R;
d2 = 0.0;
for (j = k; j > 0; j--) {
d2 = d2*zd2 + j*prj->p[j];
}
if (d2 <= 0.0) break;
zd1 = zd2;
d1 = d2;
}
if (i == 180) {
/* No negative derivative -> no point of inflection. */
zd = PI;
} else {
/* Find where the derivative is zero. */
for (i = 1; i <= 10; i++) {
zd = zd1 - d1*(zd2-zd1)/(d2-d1);
d = 0.0;
for (j = k; j > 0; j--) {
d = d*zd + j*prj->p[j];
}
if (fabs(d) < tol) break;
if (d < 0.0) {
zd2 = zd;
d2 = d;
} else {
zd1 = zd;
d1 = d;
}
}
}
r = 0.0;
for (j = k; j >= 0; j--) {
r = r*zd + prj->p[j];
}
prj->w[0] = zd;
prj->w[1] = r;
}
if (prj->flag == -1) {
prj->flag = -PRJSET;
} else {
prj->flag = PRJSET;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int zpnfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
int j;
double r, s;
if (abs(prj->flag) != PRJSET) {
if (zpnset(prj)) return 1;
}
s = (90.0 - theta)*D2R;
r = 0.0;
for (j = prj->n; j >= 0; j--) {
r = r*s + prj->p[j];
}
r = prj->r0*r;
*x = r*wcs_sind(phi);
*y = -r*wcs_cosd(phi);
if (prj->flag == PRJSET && s > prj->w[0]) {
return 2;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int zpnrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
int i, j, k;
double a, b, c, d, lambda, r, r1, r2, rt, zd, zd1, zd2;
const double tol = 1.0e-13;
if (abs(prj->flag) != PRJSET) {
if (zpnset(prj)) return 1;
}
k = prj->n;
r = sqrt(x*x + y*y)/prj->r0;
if (k < 1) {
/* Constant - no solution. */
return 1;
} else if (k == 1) {
/* Linear. */
zd = (r - prj->p[0])/prj->p[1];
} else if (k == 2) {
/* Quadratic. */
a = prj->p[2];
b = prj->p[1];
c = prj->p[0] - r;
d = b*b - 4.0*a*c;
if (d < 0.0) {
return 2;
}
d = sqrt(d);
/* Choose solution closest to pole. */
zd1 = (-b + d)/(2.0*a);
zd2 = (-b - d)/(2.0*a);
zd = (zd1<zd2) ? zd1 : zd2;
if (zd < -tol) zd = (zd1>zd2) ? zd1 : zd2;
if (zd < 0.0) {
if (zd < -tol) {
return 2;
}
zd = 0.0;
} else if (zd > PI) {
if (zd > PI+tol) {
return 2;
}
zd = PI;
}
} else {
/* Higher order - solve iteratively. */
zd1 = 0.0;
r1 = prj->p[0];
zd2 = prj->w[0];
r2 = prj->w[1];
if (r < r1) {
if (r < r1-tol) {
return 2;
}
zd = zd1;
} else if (r > r2) {
if (r > r2+tol) {
return 2;
}
zd = zd2;
} else {
/* Disect the interval. */
for (j = 0; j < 100; j++) {
lambda = (r2 - r)/(r2 - r1);
if (lambda < 0.1) {
lambda = 0.1;
} else if (lambda > 0.9) {
lambda = 0.9;
}
zd = zd2 - lambda*(zd2 - zd1);
rt = 0.0;
for (i = k; i >= 0; i--) {
rt = (rt * zd) + prj->p[i];
}
if (rt < r) {
if (r-rt < tol) break;
r1 = rt;
zd1 = zd;
} else {
if (rt-r < tol) break;
r2 = rt;
zd2 = zd;
}
if (fabs(zd2-zd1) < tol) break;
}
}
}
if (r == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(x, -y);
}
*theta = 90.0 - zd*R2D;
return 0;
}
/*============================================================================
* ZEA: zenithal/azimuthal equal area projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] 2*r0
* prj->w[1] 1/(2*r0)
*===========================================================================*/
int zeaset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 360.0/PI;
prj->w[1] = PI/360.0;
} else {
prj->w[0] = 2.0*prj->r0;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int zeafwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double r;
if (prj->flag != PRJSET) {
if (zeaset(prj)) return 1;
}
r = prj->w[0]*wcs_sind((90.0 - theta)/2.0);
*x = r*wcs_sind(phi);
*y = -r*wcs_cosd(phi);
return 0;
}
/*--------------------------------------------------------------------------*/
int zearev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double r, s;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (zeaset(prj)) return 1;
}
r = sqrt(x*x + y*y);
if (r == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(x, -y);
}
s = r*prj->w[1];
if (fabs(s) > 1.0) {
if (fabs(r - prj->w[0]) < tol) {
*theta = -90.0;
} else {
return 2;
}
} else {
*theta = 90.0 - 2.0*wcs_asind(s);
}
return 0;
}
/*============================================================================
* AIR: Airy's projection.
*
* Given:
* prj->p[1] Latitude theta_b within which the error is minimized,
* in degrees.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] 2*r0
* prj->w[1] ln(cos(xi_b))/tan(xi_b)**2, where xi_b = (90-theta_b)/2
* prj->w[2] 1/2 - prj->w[1]
* prj->w[3] 2*r0*prj->w[2]
* prj->w[4] tol, cutoff for using small angle approximation, in
* radians.
* prj->w[5] prj->w[2]*tol
* prj->w[6] (180/pi)/prj->w[2]
*===========================================================================*/
int airset(prj)
struct prjprm *prj;
{
const double tol = 1.0e-4;
double cxi;
if (prj->r0 == 0.0) prj->r0 = R2D;
prj->w[0] = 2.0*prj->r0;
if (prj->p[1] == 90.0) {
prj->w[1] = -0.5;
prj->w[2] = 1.0;
} else if (prj->p[1] > -90.0) {
cxi = wcs_cosd((90.0 - prj->p[1])/2.0);
prj->w[1] = log(cxi)*(cxi*cxi)/(1.0-cxi*cxi);
prj->w[2] = 0.5 - prj->w[1];
} else {
return 1;
}
prj->w[3] = prj->w[0] * prj->w[2];
prj->w[4] = tol;
prj->w[5] = prj->w[2]*tol;
prj->w[6] = R2D/prj->w[2];
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int airfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double cxi, r, txi, xi;
if (prj->flag != PRJSET) {
if (airset(prj)) return 1;
}
if (theta == 90.0) {
r = 0.0;
} else if (theta > -90.0) {
xi = D2R*(90.0 - theta)/2.0;
if (xi < prj->w[4]) {
r = xi*prj->w[3];
} else {
cxi = wcs_cosd((90.0 - theta)/2.0);
txi = sqrt(1.0-cxi*cxi)/cxi;
r = -prj->w[0]*(log(cxi)/txi + prj->w[1]*txi);
}
} else {
return 2;
}
*x = r*wcs_sind(phi);
*y = -r*wcs_cosd(phi);
return 0;
}
/*--------------------------------------------------------------------------*/
int airrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
int j;
double cxi, lambda, r, r1, r2, rt, txi, x1, x2, xi;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (airset(prj)) return 1;
}
r = sqrt(x*x + y*y)/prj->w[0];
if (r == 0.0) {
xi = 0.0;
} else if (r < prj->w[5]) {
xi = r*prj->w[6];
} else {
/* Find a solution interval. */
x1 = 1.0;
r1 = 0.0;
for (j = 0; j < 30; j++) {
x2 = x1/2.0;
txi = sqrt(1.0-x2*x2)/x2;
r2 = -(log(x2)/txi + prj->w[1]*txi);
if (r2 >= r) break;
x1 = x2;
r1 = r2;
}
if (j == 30) return 2;
for (j = 0; j < 100; j++) {
/* Weighted division of the interval. */
lambda = (r2-r)/(r2-r1);
if (lambda < 0.1) {
lambda = 0.1;
} else if (lambda > 0.9) {
lambda = 0.9;
}
cxi = x2 - lambda*(x2-x1);
txi = sqrt(1.0-cxi*cxi)/cxi;
rt = -(log(cxi)/txi + prj->w[1]*txi);
if (rt < r) {
if (r-rt < tol) break;
r1 = rt;
x1 = cxi;
} else {
if (rt-r < tol) break;
r2 = rt;
x2 = cxi;
}
}
if (j == 100) return 2;
xi = wcs_acosd(cxi);
}
if (r == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(x, -y);
}
*theta = 90.0 - 2.0*xi;
return 0;
}
/*============================================================================
* CYP: cylindrical perspective projection.
*
* Given:
* prj->p[1] Distance of point of projection from the centre of the
* generating sphere, mu, in units of r0.
* prj->p[2] Radius of the cylinder of projection, lambda, in units
* of r0.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*lambda*(pi/180)
* prj->w[1] (180/pi)/(r0*lambda)
* prj->w[2] r0*(mu + lambda)
* prj->w[3] 1/(r0*(mu + lambda))
*===========================================================================*/
int cypset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = prj->p[2];
if (prj->w[0] == 0.0) {
return 1;
}
prj->w[1] = 1.0/prj->w[0];
prj->w[2] = R2D*(prj->p[1] + prj->p[2]);
if (prj->w[2] == 0.0) {
return 1;
}
prj->w[3] = 1.0/prj->w[2];
} else {
prj->w[0] = prj->r0*prj->p[2]*D2R;
if (prj->w[0] == 0.0) {
return 1;
}
prj->w[1] = 1.0/prj->w[0];
prj->w[2] = prj->r0*(prj->p[1] + prj->p[2]);
if (prj->w[2] == 0.0) {
return 1;
}
prj->w[3] = 1.0/prj->w[2];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int cypfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double s;
if (prj->flag != PRJSET) {
if (cypset(prj)) return 1;
}
s = prj->p[1] + wcs_cosd(theta);
if (s == 0.0) {
return 2;
}
*x = prj->w[0]*phi;
*y = prj->w[2]*wcs_sind(theta)/s;
return 0;
}
/*--------------------------------------------------------------------------*/
int cyprev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double eta;
if (prj->flag != PRJSET) {
if (cypset(prj)) return 1;
}
*phi = x*prj->w[1];
eta = y*prj->w[3];
*theta = wcs_atan2d(eta,1.0) + wcs_asind(eta*prj->p[1]/sqrt(eta*eta+1.0));
return 0;
}
/*============================================================================
* CAR: Cartesian projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/180)
* prj->w[1] (180/pi)/r0
*===========================================================================*/
int carset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 1.0;
prj->w[1] = 1.0;
} else {
prj->w[0] = prj->r0*D2R;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int carfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
if (prj->flag != PRJSET) {
if (carset(prj)) return 1;
}
*x = prj->w[0]*phi;
*y = prj->w[0]*theta;
return 0;
}
/*--------------------------------------------------------------------------*/
int carrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
if (prj->flag != PRJSET) {
if (carset(prj)) return 1;
}
*phi = prj->w[1]*x;
*theta = prj->w[1]*y;
return 0;
}
/*============================================================================
* MER: Mercator's projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/180)
* prj->w[1] (180/pi)/r0
*===========================================================================*/
int merset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 1.0;
prj->w[1] = 1.0;
} else {
prj->w[0] = prj->r0*D2R;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int merfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
if (prj->flag != PRJSET) {
if (merset(prj)) return 1;
}
if (theta <= -90.0 || theta >= 90.0) {
return 2;
}
*x = prj->w[0]*phi;
*y = prj->r0*log(wcs_tand((90.0+theta)/2.0));
return 0;
}
/*--------------------------------------------------------------------------*/
int merrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
if (prj->flag != PRJSET) {
if (merset(prj)) return 1;
}
*phi = x*prj->w[1];
*theta = 2.0*wcs_atand(exp(y/prj->r0)) - 90.0;
return 0;
}
/*============================================================================
* CEA: cylindrical equal area projection.
*
* Given:
* prj->p[1] Square of the cosine of the latitude at which the
* projection is conformal, lambda.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/180)
* prj->w[1] (180/pi)/r0
* prj->w[2] r0/lambda
* prj->w[3] lambda/r0
*===========================================================================*/
int ceaset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 1.0;
prj->w[1] = 1.0;
if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) {
return 1;
}
prj->w[2] = prj->r0/prj->p[1];
prj->w[3] = prj->p[1]/prj->r0;
} else {
prj->w[0] = prj->r0*D2R;
prj->w[1] = R2D/prj->r0;
if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) {
return 1;
}
prj->w[2] = prj->r0/prj->p[1];
prj->w[3] = prj->p[1]/prj->r0;
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int ceafwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
if (prj->flag != PRJSET) {
if (ceaset(prj)) return 1;
}
*x = prj->w[0]*phi;
*y = prj->w[2]*wcs_sind(theta);
return 0;
}
/*--------------------------------------------------------------------------*/
int cearev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double s;
const double tol = 1.0e-13;
if (prj->flag != PRJSET) {
if (ceaset(prj)) return 1;
}
s = y*prj->w[3];
if (fabs(s) > 1.0) {
if (fabs(s) > 1.0+tol) {
return 2;
}
s = copysign(1.0,s);
}
*phi = x*prj->w[1];
*theta = wcs_asind(s);
return 0;
}
/*============================================================================
* COP: conic perspective projection.
*
* Given:
* prj->p[1] sigma = (theta2+theta1)/2
* prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
* latitudes of the standard parallels, in degrees.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] C = sin(sigma)
* prj->w[1] 1/C
* prj->w[2] Y0 = r0*cos(delta)*cot(sigma)
* prj->w[3] r0*cos(delta)
* prj->w[4] 1/(r0*cos(delta)
* prj->w[5] cot(sigma)
*===========================================================================*/
int copset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) prj->r0 = R2D;
prj->w[0] = wcs_sind(prj->p[1]);
if (prj->w[0] == 0.0) {
return 1;
}
prj->w[1] = 1.0/prj->w[0];
prj->w[3] = prj->r0*wcs_cosd(prj->p[2]);
if (prj->w[3] == 0.0) {
return 1;
}
prj->w[4] = 1.0/prj->w[3];
prj->w[5] = 1.0/wcs_tand(prj->p[1]);
prj->w[2] = prj->w[3]*prj->w[5];
if (prj->flag == -1) {
prj->flag = -PRJSET;
} else {
prj->flag = PRJSET;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int copfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double a, r, s, t;
if (abs(prj->flag) != PRJSET) {
if (copset(prj)) return 1;
}
t = theta - prj->p[1];
s = wcs_cosd(t);
if (s == 0.0) {
return 2;
}
a = prj->w[0]*phi;
r = prj->w[2] - prj->w[3]*wcs_sind(t)/s;
*x = r*wcs_sind(a);
*y = prj->w[2] - r*wcs_cosd(a);
if (prj->flag == PRJSET && r*prj->w[0] < 0.0) {
return 2;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int coprev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double a, dy, r;
if (abs(prj->flag) != PRJSET) {
if (copset(prj)) return 1;
}
dy = prj->w[2] - y;
r = sqrt(x*x + dy*dy);
if (prj->p[1] < 0.0) r = -r;
if (r == 0.0) {
a = 0.0;
} else {
a = wcs_atan2d(x/r, dy/r);
}
*phi = a*prj->w[1];
*theta = prj->p[1] + wcs_atand(prj->w[5] - r*prj->w[4]);
return 0;
}
/*============================================================================
* COD: conic equidistant projection.
*
* Given:
* prj->p[1] sigma = (theta2+theta1)/2
* prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
* latitudes of the standard parallels, in degrees.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] C = r0*sin(sigma)*sin(delta)/delta
* prj->w[1] 1/C
* prj->w[2] Y0 = delta*cot(delta)*cot(sigma)
* prj->w[3] Y0 + sigma
*===========================================================================*/
int codset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) prj->r0 = R2D;
if (prj->p[2] == 0.0) {
prj->w[0] = prj->r0*wcs_sind(prj->p[1])*D2R;
} else {
prj->w[0] = prj->r0*wcs_sind(prj->p[1])*wcs_sind(prj->p[2])/prj->p[2];
}
if (prj->w[0] == 0.0) {
return 1;
}
prj->w[1] = 1.0/prj->w[0];
prj->w[2] = prj->r0*wcs_cosd(prj->p[2])*wcs_cosd(prj->p[1])/prj->w[0];
prj->w[3] = prj->w[2] + prj->p[1];
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int codfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double a, r;
if (prj->flag != PRJSET) {
if (codset(prj)) return 1;
}
a = prj->w[0]*phi;
r = prj->w[3] - theta;
*x = r*wcs_sind(a);
*y = prj->w[2] - r*wcs_cosd(a);
return 0;
}
/*--------------------------------------------------------------------------*/
int codrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double a, dy, r;
if (prj->flag != PRJSET) {
if (codset(prj)) return 1;
}
dy = prj->w[2] - y;
r = sqrt(x*x + dy*dy);
if (prj->p[1] < 0.0) r = -r;
if (r == 0.0) {
a = 0.0;
} else {
a = wcs_atan2d(x/r, dy/r);
}
*phi = a*prj->w[1];
*theta = prj->w[3] - r;
return 0;
}
/*============================================================================
* COE: conic equal area projection.
*
* Given:
* prj->p[1] sigma = (theta2+theta1)/2
* prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
* latitudes of the standard parallels, in degrees.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] C = (sin(theta1) + sin(theta2))/2
* prj->w[1] 1/C
* prj->w[2] Y0 = chi*sqrt(psi - 2C*wcs_sind(sigma))
* prj->w[3] chi = r0/C
* prj->w[4] psi = 1 + sin(theta1)*sin(theta2)
* prj->w[5] 2C
* prj->w[6] (1 + sin(theta1)*sin(theta2))*(r0/C)**2
* prj->w[7] C/(2*r0**2)
* prj->w[8] chi*sqrt(psi + 2C)
*===========================================================================*/
int coeset(prj)
struct prjprm *prj;
{
double theta1, theta2;
if (prj->r0 == 0.0) prj->r0 = R2D;
theta1 = prj->p[1] - prj->p[2];
theta2 = prj->p[1] + prj->p[2];
prj->w[0] = (wcs_sind(theta1) + wcs_sind(theta2))/2.0;
if (prj->w[0] == 0.0) {
return 1;
}
prj->w[1] = 1.0/prj->w[0];
prj->w[3] = prj->r0/prj->w[0];
prj->w[4] = 1.0 + wcs_sind(theta1)*wcs_sind(theta2);
prj->w[5] = 2.0*prj->w[0];
prj->w[6] = prj->w[3]*prj->w[3]*prj->w[4];
prj->w[7] = 1.0/(2.0*prj->r0*prj->w[3]);
prj->w[8] = prj->w[3]*sqrt(prj->w[4] + prj->w[5]);
prj->w[2] = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*wcs_sind(prj->p[1]));
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int coefwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double a, r;
if (prj->flag != PRJSET) {
if (coeset(prj)) return 1;
}
a = phi*prj->w[0];
if (theta == -90.0) {
r = prj->w[8];
} else {
r = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*wcs_sind(theta));
}
*x = r*wcs_sind(a);
*y = prj->w[2] - r*wcs_cosd(a);
return 0;
}
/*--------------------------------------------------------------------------*/
int coerev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double a, dy, r, w;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (coeset(prj)) return 1;
}
dy = prj->w[2] - y;
r = sqrt(x*x + dy*dy);
if (prj->p[1] < 0.0) r = -r;
if (r == 0.0) {
a = 0.0;
} else {
a = wcs_atan2d(x/r, dy/r);
}
*phi = a*prj->w[1];
if (fabs(r - prj->w[8]) < tol) {
*theta = -90.0;
} else {
w = (prj->w[6] - r*r)*prj->w[7];
if (fabs(w) > 1.0) {
if (fabs(w-1.0) < tol) {
*theta = 90.0;
} else if (fabs(w+1.0) < tol) {
*theta = -90.0;
} else {
return 2;
}
} else {
*theta = wcs_asind(w);
}
}
return 0;
}
/*============================================================================
* COO: conic orthomorphic projection.
*
* Given:
* prj->p[1] sigma = (theta2+theta1)/2
* prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
* latitudes of the standard parallels, in degrees.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] C = ln(cos(theta2)/cos(theta1))/ln(tan(tau2)/tan(tau1))
* where tau1 = (90 - theta1)/2
* tau2 = (90 - theta2)/2
* prj->w[1] 1/C
* prj->w[2] Y0 = psi*tan((90-sigma)/2)**C
* prj->w[3] psi = (r0*cos(theta1)/C)/tan(tau1)**C
* prj->w[4] 1/psi
*===========================================================================*/
int cooset(prj)
struct prjprm *prj;
{
double cos1, cos2, tan1, tan2, theta1, theta2;
if (prj->r0 == 0.0) prj->r0 = R2D;
theta1 = prj->p[1] - prj->p[2];
theta2 = prj->p[1] + prj->p[2];
tan1 = wcs_tand((90.0 - theta1)/2.0);
cos1 = wcs_cosd(theta1);
if (theta1 == theta2) {
prj->w[0] = wcs_sind(theta1);
} else {
tan2 = wcs_tand((90.0 - theta2)/2.0);
cos2 = wcs_cosd(theta2);
prj->w[0] = log(cos2/cos1)/log(tan2/tan1);
}
if (prj->w[0] == 0.0) {
return 1;
}
prj->w[1] = 1.0/prj->w[0];
prj->w[3] = prj->r0*(cos1/prj->w[0])/pow(tan1,prj->w[0]);
if (prj->w[3] == 0.0) {
return 1;
}
prj->w[2] = prj->w[3]*pow(wcs_tand((90.0 - prj->p[1])/2.0),prj->w[0]);
prj->w[4] = 1.0/prj->w[3];
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int coofwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double a, r;
if (prj->flag != PRJSET) {
if (cooset(prj)) return 1;
}
a = prj->w[0]*phi;
if (theta == -90.0) {
if (prj->w[0] < 0.0) {
r = 0.0;
} else {
return 2;
}
} else {
r = prj->w[3]*pow(wcs_tand((90.0 - theta)/2.0),prj->w[0]);
}
*x = r*wcs_sind(a);
*y = prj->w[2] - r*wcs_cosd(a);
return 0;
}
/*--------------------------------------------------------------------------*/
int coorev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double a, dy, r;
if (prj->flag != PRJSET) {
if (cooset(prj)) return 1;
}
dy = prj->w[2] - y;
r = sqrt(x*x + dy*dy);
if (prj->p[1] < 0.0) r = -r;
if (r == 0.0) {
a = 0.0;
} else {
a = wcs_atan2d(x/r, dy/r);
}
*phi = a*prj->w[1];
if (r == 0.0) {
if (prj->w[0] < 0.0) {
*theta = -90.0;
} else {
return 2;
}
} else {
*theta = 90.0 - 2.0*wcs_atand(pow(r*prj->w[4],prj->w[1]));
}
return 0;
}
/*============================================================================
* BON: Bonne's projection.
*
* Given:
* prj->p[1] Bonne conformal latitude, theta1, in degrees.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[1] r0*pi/180
* prj->w[2] Y0 = r0*cot(theta1) + theta1*pi/180)
*===========================================================================*/
int bonset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[1] = 1.0;
prj->w[2] = prj->r0*wcs_cosd(prj->p[1])/wcs_sind(prj->p[1]) + prj->p[1];
} else {
prj->w[1] = prj->r0*D2R;
prj->w[2] = prj->r0*(wcs_cosd(prj->p[1])/wcs_sind(prj->p[1]) + prj->p[1]*D2R);
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int bonfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double a, r;
if (prj->p[1] == 0.0) {
/* Sanson-Flamsteed. */
return glsfwd(phi, theta, prj, x, y);
}
if (prj->flag != PRJSET) {
if (bonset(prj)) return 1;
}
r = prj->w[2] - theta*prj->w[1];
a = prj->r0*phi*wcs_cosd(theta)/r;
*x = r*wcs_sind(a);
*y = prj->w[2] - r*wcs_cosd(a);
return 0;
}
/*--------------------------------------------------------------------------*/
int bonrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double a, dy, costhe, r;
if (prj->p[1] == 0.0) {
/* Sanson-Flamsteed. */
return glsrev(x, y, prj, phi, theta);
}
if (prj->flag != PRJSET) {
if (bonset(prj)) return 1;
}
dy = prj->w[2] - y;
r = sqrt(x*x + dy*dy);
if (prj->p[1] < 0.0) r = -r;
if (r == 0.0) {
a = 0.0;
} else {
a = wcs_atan2d(x/r, dy/r);
}
*theta = (prj->w[2] - r)/prj->w[1];
costhe = wcs_cosd(*theta);
if (costhe == 0.0) {
*phi = 0.0;
} else {
*phi = a*(r/prj->r0)/wcs_cosd(*theta);
}
return 0;
}
/*============================================================================
* PCO: polyconic projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/180)
* prj->w[1] 1/r0
* prj->w[2] 2*r0
*===========================================================================*/
int pcoset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 1.0;
prj->w[1] = 1.0;
prj->w[2] = 360.0/PI;
} else {
prj->w[0] = prj->r0*D2R;
prj->w[1] = 1.0/prj->w[0];
prj->w[2] = 2.0*prj->r0;
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int pcofwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double a, costhe, cotthe, sinthe;
if (prj->flag != PRJSET) {
if (pcoset(prj)) return 1;
}
costhe = wcs_cosd(theta);
sinthe = wcs_sind(theta);
a = phi*sinthe;
if (sinthe == 0.0) {
*x = prj->w[0]*phi;
*y = 0.0;
} else {
cotthe = costhe/sinthe;
*x = prj->r0*cotthe*wcs_sind(a);
*y = prj->r0*(cotthe*(1.0 - wcs_cosd(a)) + theta*D2R);
}
return 0;
}
/*--------------------------------------------------------------------------*/
int pcorev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
int j;
double f, fneg, fpos, lambda, tanthe, theneg, thepos, w, xp, xx, ymthe, yp;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (pcoset(prj)) return 1;
}
w = fabs(y*prj->w[1]);
if (w < tol) {
*phi = x*prj->w[1];
*theta = 0.0;
} else if (fabs(w-90.0) < tol) {
*phi = 0.0;
*theta = wcs_copysign(90.0,y);
} else {
/* Iterative solution using weighted division of the interval. */
if (y > 0.0) {
thepos = 90.0;
} else {
thepos = -90.0;
}
theneg = 0.0;
xx = x*x;
ymthe = y - prj->w[0]*thepos;
fpos = xx + ymthe*ymthe;
fneg = -999.0;
for (j = 0; j < 64; j++) {
if (fneg < -100.0) {
/* Equal division of the interval. */
*theta = (thepos+theneg)/2.0;
} else {
/* Weighted division of the interval. */
lambda = fpos/(fpos-fneg);
if (lambda < 0.1) {
lambda = 0.1;
} else if (lambda > 0.9) {
lambda = 0.9;
}
*theta = thepos - lambda*(thepos-theneg);
}
/* Compute the residue. */
ymthe = y - prj->w[0]*(*theta);
tanthe = wcs_tand(*theta);
f = xx + ymthe*(ymthe - prj->w[2]/tanthe);
/* Check for convergence. */
if (fabs(f) < tol) break;
if (fabs(thepos-theneg) < tol) break;
/* Redefine the interval. */
if (f > 0.0) {
thepos = *theta;
fpos = f;
} else {
theneg = *theta;
fneg = f;
}
}
xp = prj->r0 - ymthe*tanthe;
yp = x*tanthe;
if (xp == 0.0 && yp == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(yp, xp)/wcs_sind(*theta);
}
}
return 0;
}
/*============================================================================
* GLS: Sanson-Flamsteed ("global sinusoid") projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/180)
* prj->w[1] (180/pi)/r0
*===========================================================================*/
int glsset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 1.0;
prj->w[1] = 1.0;
} else {
prj->w[0] = prj->r0*D2R;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int glsfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
if (prj->flag != PRJSET) {
if (glsset(prj)) return 1;
}
*x = prj->w[0]*phi*wcs_cosd(theta);
*y = prj->w[0]*theta;
return 0;
}
/*--------------------------------------------------------------------------*/
int glsrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double w;
if (prj->flag != PRJSET) {
if (glsset(prj)) return 1;
}
w = cos(y/prj->r0);
if (w == 0.0) {
*phi = 0.0;
} else {
*phi = x*prj->w[1]/cos(y/prj->r0);
}
*theta = y*prj->w[1];
return 0;
}
/*============================================================================
* PAR: parabolic projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/180)
* prj->w[1] (180/pi)/r0
* prj->w[2] pi*r0
* prj->w[3] 1/(pi*r0)
*===========================================================================*/
int parset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 1.0;
prj->w[1] = 1.0;
prj->w[2] = 180.0;
prj->w[3] = 1.0/prj->w[2];
} else {
prj->w[0] = prj->r0*D2R;
prj->w[1] = 1.0/prj->w[0];
prj->w[2] = PI*prj->r0;
prj->w[3] = 1.0/prj->w[2];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int parfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double s;
if (prj->flag != PRJSET) {
if (parset(prj)) return 1;
}
s = wcs_sind(theta/3.0);
*x = prj->w[0]*phi*(1.0 - 4.0*s*s);
*y = prj->w[2]*s;
return 0;
}
/*--------------------------------------------------------------------------*/
int parrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double s, t;
if (prj->flag != PRJSET) {
if (parset(prj)) return 1;
}
s = y*prj->w[3];
if (s > 1.0 || s < -1.0) {
return 2;
}
t = 1.0 - 4.0*s*s;
if (t == 0.0) {
if (x == 0.0) {
*phi = 0.0;
} else {
return 2;
}
} else {
*phi = prj->w[1]*x/t;
}
*theta = 3.0*wcs_asind(s);
return 0;
}
/*============================================================================
* AIT: Hammer-Aitoff projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] 2*r0**2
* prj->w[1] 1/(2*r0)**2
* prj->w[2] 1/(4*r0)**2
* prj->w[3] 1/(2*r0)
*===========================================================================*/
int aitset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) prj->r0 = R2D;
prj->w[0] = 2.0*prj->r0*prj->r0;
prj->w[1] = 1.0/(2.0*prj->w[0]);
prj->w[2] = prj->w[1]/4.0;
prj->w[3] = 1.0/(2.0*prj->r0);
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int aitfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double costhe, w;
if (prj->flag != PRJSET) {
if (aitset(prj)) return 1;
}
costhe = wcs_cosd(theta);
w = sqrt(prj->w[0]/(1.0 + costhe*wcs_cosd(phi/2.0)));
*x = 2.0*w*costhe*wcs_sind(phi/2.0);
*y = w*wcs_sind(theta);
return 0;
}
/*--------------------------------------------------------------------------*/
int aitrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double s, u, xp, yp, z;
const double tol = 1.0e-13;
if (prj->flag != PRJSET) {
if (aitset(prj)) return 1;
}
u = 1.0 - x*x*prj->w[2] - y*y*prj->w[1];
if (u < 0.0) {
if (u < -tol) {
return 2;
}
u = 0.0;
}
z = sqrt(u);
s = z*y/prj->r0;
if (fabs(s) > 1.0) {
if (fabs(s) > 1.0+tol) {
return 2;
}
s = wcs_copysign(1.0,s);
}
xp = 2.0*z*z - 1.0;
yp = z*x*prj->w[3];
if (xp == 0.0 && yp == 0.0) {
*phi = 0.0;
} else {
*phi = 2.0*wcs_atan2d(yp, xp);
}
*theta = wcs_asind(s);
return 0;
}
/*============================================================================
* MOL: Mollweide's projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] sqrt(2)*r0
* prj->w[1] sqrt(2)*r0/90
* prj->w[2] 1/(sqrt(2)*r0)
* prj->w[3] 90/r0
*===========================================================================*/
int molset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) prj->r0 = R2D;
prj->w[0] = SQRT2*prj->r0;
prj->w[1] = prj->w[0]/90.0;
prj->w[2] = 1.0/prj->w[0];
prj->w[3] = 90.0/prj->r0;
prj->w[4] = 2.0/PI;
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int molfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
int j;
double alpha, resid, u, v, v0, v1;
const double tol = 1.0e-13;
if (prj->flag != PRJSET) {
if (molset(prj)) return 1;
}
if (fabs(theta) == 90.0) {
*x = 0.0;
*y = wcs_copysign(prj->w[0],theta);
} else if (theta == 0.0) {
*x = prj->w[1]*phi;
*y = 0.0;
} else {
u = PI*wcs_sind(theta);
v0 = -PI;
v1 = PI;
v = u;
for (j = 0; j < 100; j++) {
resid = (v - u) + sin(v);
if (resid < 0.0) {
if (resid > -tol) break;
v0 = v;
} else {
if (resid < tol) break;
v1 = v;
}
v = (v0 + v1)/2.0;
}
alpha = v/2.0;
*x = prj->w[1]*phi*cos(alpha);
*y = prj->w[0]*sin(alpha);
}
return 0;
}
/*--------------------------------------------------------------------------*/
int molrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double s, y0, z;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (molset(prj)) return 1;
}
y0 = y/prj->r0;
s = 2.0 - y0*y0;
if (s <= tol) {
if (s < -tol) {
return 2;
}
s = 0.0;
if (fabs(x) > tol) {
return 2;
}
*phi = 0.0;
} else {
s = sqrt(s);
*phi = prj->w[3]*x/s;
}
z = y*prj->w[2];
if (fabs(z) > 1.0) {
if (fabs(z) > 1.0+tol) {
return 2;
}
z = wcs_copysign(1.0,z) + y0*s/PI;
} else {
z = asin(z)*prj->w[4] + y0*s/PI;
}
if (fabs(z) > 1.0) {
if (fabs(z) > 1.0+tol) {
return 2;
}
z = wcs_copysign(1.0,z);
}
*theta = wcs_asind(z);
return 0;
}
/*============================================================================
* CSC: COBE quadrilateralized spherical cube projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/4)
* prj->w[1] (4/pi)/r0
*===========================================================================*/
int cscset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 45.0;
prj->w[1] = 1.0/45.0;
} else {
prj->w[0] = prj->r0*PI/4.0;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int cscfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
int face;
double costhe, eta, l, m, n, rho, xi;
const float tol = 1.0e-7;
float a, a2, a2b2, a4, ab, b, b2, b4, ca2, cb2, x0, xf, y0, yf;
const float gstar = 1.37484847732;
const float mm = 0.004869491981;
const float gamma = -0.13161671474;
const float omega1 = -0.159596235474;
const float d0 = 0.0759196200467;
const float d1 = -0.0217762490699;
const float c00 = 0.141189631152;
const float c10 = 0.0809701286525;
const float c01 = -0.281528535557;
const float c11 = 0.15384112876;
const float c20 = -0.178251207466;
const float c02 = 0.106959469314;
if (prj->flag != PRJSET) {
if (cscset(prj)) return 1;
}
costhe = wcs_cosd(theta);
l = costhe*wcs_cosd(phi);
m = costhe*wcs_sind(phi);
n = wcs_sind(theta);
face = 0;
rho = n;
if (l > rho) {
face = 1;
rho = l;
}
if (m > rho) {
face = 2;
rho = m;
}
if (-l > rho) {
face = 3;
rho = -l;
}
if (-m > rho) {
face = 4;
rho = -m;
}
if (-n > rho) {
face = 5;
rho = -n;
}
if (face == 0) {
xi = m;
eta = -l;
x0 = 0.0;
y0 = 2.0;
} else if (face == 1) {
xi = m;
eta = n;
x0 = 0.0;
y0 = 0.0;
} else if (face == 2) {
xi = -l;
eta = n;
x0 = 2.0;
y0 = 0.0;
} else if (face == 3) {
xi = -m;
eta = n;
x0 = 4.0;
y0 = 0.0;
} else if (face == 4) {
xi = l;
eta = n;
x0 = 6.0;
y0 = 0.0;
} else {
xi = m;
eta = l;
x0 = 0.0;
y0 = -2.0;
}
a = xi/rho;
b = eta/rho;
a2 = a*a;
b2 = b*b;
ca2 = 1.0 - a2;
cb2 = 1.0 - b2;
/* Avoid floating underflows. */
ab = fabs(a*b);
a4 = (a2 > 1.0e-16) ? a2*a2 : 0.0;
b4 = (b2 > 1.0e-16) ? b2*b2 : 0.0;
a2b2 = (ab > 1.0e-16) ? a2*b2 : 0.0;
xf = a*(a2 + ca2*(gstar + b2*(gamma*ca2 + mm*a2 +
cb2*(c00 + c10*a2 + c01*b2 + c11*a2b2 + c20*a4 + c02*b4)) +
a2*(omega1 - ca2*(d0 + d1*a2))));
yf = b*(b2 + cb2*(gstar + a2*(gamma*cb2 + mm*b2 +
ca2*(c00 + c10*b2 + c01*a2 + c11*a2b2 + c20*b4 + c02*a4)) +
b2*(omega1 - cb2*(d0 + d1*b2))));
if (fabs(xf) > 1.0) {
if (fabs(xf) > 1.0+tol) {
return 2;
}
xf = wcs_copysign(1.0,xf);
}
if (fabs(yf) > 1.0) {
if (fabs(yf) > 1.0+tol) {
return 2;
}
yf = wcs_copysign(1.0,yf);
}
*x = prj->w[0]*(x0 + xf);
*y = prj->w[0]*(y0 + yf);
return 0;
}
/*--------------------------------------------------------------------------*/
int cscrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
int face;
double l, m, n;
float a, b, xf, xx, yf, yy, z0, z1, z2, z3, z4, z5, z6;
const float p00 = -0.27292696;
const float p10 = -0.07629969;
const float p20 = -0.22797056;
const float p30 = 0.54852384;
const float p40 = -0.62930065;
const float p50 = 0.25795794;
const float p60 = 0.02584375;
const float p01 = -0.02819452;
const float p11 = -0.01471565;
const float p21 = 0.48051509;
const float p31 = -1.74114454;
const float p41 = 1.71547508;
const float p51 = -0.53022337;
const float p02 = 0.27058160;
const float p12 = -0.56800938;
const float p22 = 0.30803317;
const float p32 = 0.98938102;
const float p42 = -0.83180469;
const float p03 = -0.60441560;
const float p13 = 1.50880086;
const float p23 = -0.93678576;
const float p33 = 0.08693841;
const float p04 = 0.93412077;
const float p14 = -1.41601920;
const float p24 = 0.33887446;
const float p05 = -0.63915306;
const float p15 = 0.52032238;
const float p06 = 0.14381585;
if (prj->flag != PRJSET) {
if (cscset(prj)) return 1;
}
xf = x*prj->w[1];
yf = y*prj->w[1];
/* Check bounds. */
if (fabs(xf) <= 1.0) {
if (fabs(yf) > 3.0) return 2;
} else {
if (fabs(xf) > 7.0) return 2;
if (fabs(yf) > 1.0) return 2;
}
/* Map negative faces to the other side. */
if (xf < -1.0) xf += 8.0;
/* Determine the face. */
if (xf > 5.0) {
face = 4;
xf = xf - 6.0;
} else if (xf > 3.0) {
face = 3;
xf = xf - 4.0;
} else if (xf > 1.0) {
face = 2;
xf = xf - 2.0;
} else if (yf > 1.0) {
face = 0;
yf = yf - 2.0;
} else if (yf < -1.0) {
face = 5;
yf = yf + 2.0;
} else {
face = 1;
}
xx = xf*xf;
yy = yf*yf;
z0 = p00 + xx*(p10 + xx*(p20 + xx*(p30 + xx*(p40 + xx*(p50 + xx*(p60))))));
z1 = p01 + xx*(p11 + xx*(p21 + xx*(p31 + xx*(p41 + xx*(p51)))));
z2 = p02 + xx*(p12 + xx*(p22 + xx*(p32 + xx*(p42))));
z3 = p03 + xx*(p13 + xx*(p23 + xx*(p33)));
z4 = p04 + xx*(p14 + xx*(p24));
z5 = p05 + xx*(p15);
z6 = p06;
a = z0 + yy*(z1 + yy*(z2 + yy*(z3 + yy*(z4 + yy*(z5 + yy*z6)))));
a = xf + xf*(1.0 - xx)*a;
z0 = p00 + yy*(p10 + yy*(p20 + yy*(p30 + yy*(p40 + yy*(p50 + yy*(p60))))));
z1 = p01 + yy*(p11 + yy*(p21 + yy*(p31 + yy*(p41 + yy*(p51)))));
z2 = p02 + yy*(p12 + yy*(p22 + yy*(p32 + yy*(p42))));
z3 = p03 + yy*(p13 + yy*(p23 + yy*(p33)));
z4 = p04 + yy*(p14 + yy*(p24));
z5 = p05 + yy*(p15);
z6 = p06;
b = z0 + xx*(z1 + xx*(z2 + xx*(z3 + xx*(z4 + xx*(z5 + xx*z6)))));
b = yf + yf*(1.0 - yy)*b;
if (face == 0) {
n = 1.0/sqrt(a*a + b*b + 1.0);
l = -b*n;
m = a*n;
} else if (face == 1) {
l = 1.0/sqrt(a*a + b*b + 1.0);
m = a*l;
n = b*l;
} else if (face == 2) {
m = 1.0/sqrt(a*a + b*b + 1.0);
l = -a*m;
n = b*m;
} else if (face == 3) {
l = -1.0/sqrt(a*a + b*b + 1.0);
m = a*l;
n = -b*l;
} else if (face == 4) {
m = -1.0/sqrt(a*a + b*b + 1.0);
l = -a*m;
n = -b*m;
} else {
n = -1.0/sqrt(a*a + b*b + 1.0);
l = -b*n;
m = -a*n;
}
if (l == 0.0 && m == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(m, l);
}
*theta = wcs_asind(n);
return 0;
}
/*============================================================================
* QSC: quadrilaterilized spherical cube projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/4)
* prj->w[1] (4/pi)/r0
*===========================================================================*/
int qscset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 45.0;
prj->w[1] = 1.0/45.0;
} else {
prj->w[0] = prj->r0*PI/4.0;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int qscfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
int face;
double chi, costhe, eta, l, m, n, p, psi, rho, rhu, t, x0, xf, xi, y0, yf;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (qscset(prj)) return 1;
}
if (fabs(theta) == 90.0) {
*x = 0.0;
*y = wcs_copysign(2.0*prj->w[0],theta);
return 0;
}
costhe = wcs_cosd(theta);
l = costhe*wcs_cosd(phi);
m = costhe*wcs_sind(phi);
n = wcs_sind(theta);
face = 0;
rho = n;
if (l > rho) {
face = 1;
rho = l;
}
if (m > rho) {
face = 2;
rho = m;
}
if (-l > rho) {
face = 3;
rho = -l;
}
if (-m > rho) {
face = 4;
rho = -m;
}
if (-n > rho) {
face = 5;
rho = -n;
}
rhu = 1.0 - rho;
if (face == 0) {
xi = m;
eta = -l;
if (rhu < 1.0e-8) {
/* Small angle formula. */
t = (90.0 - theta)*D2R;
rhu = t*t/2.0;
}
x0 = 0.0;
y0 = 2.0;
} else if (face == 1) {
xi = m;
eta = n;
if (rhu < 1.0e-8) {
/* Small angle formula. */
t = theta*D2R;
p = fmod(phi,360.0);
if (p < -180.0) p += 360.0;
if (p > 180.0) p -= 360.0;
p *= D2R;
rhu = (p*p + t*t)/2.0;
}
x0 = 0.0;
y0 = 0.0;
} else if (face == 2) {
xi = -l;
eta = n;
if (rhu < 1.0e-8) {
/* Small angle formula. */
t = theta*D2R;
p = fmod(phi,360.0);
if (p < -180.0) p += 360.0;
p = (90.0 - p)*D2R;
rhu = (p*p + t*t)/2.0;
}
x0 = 2.0;
y0 = 0.0;
} else if (face == 3) {
xi = -m;
eta = n;
if (rhu < 1.0e-8) {
/* Small angle formula. */
t = theta*D2R;
p = fmod(phi,360.0);
if (p < 0.0) p += 360.0;
p = (180.0 - p)*D2R;
rhu = (p*p + t*t)/2.0;
}
x0 = 4.0;
y0 = 0.0;
} else if (face == 4) {
xi = l;
eta = n;
if (rhu < 1.0e-8) {
/* Small angle formula. */
t = theta*D2R;
p = fmod(phi,360.0);
if (p > 180.0) p -= 360.0;
p *= (90.0 + p)*D2R;
rhu = (p*p + t*t)/2.0;
}
x0 = 6;
y0 = 0.0;
} else {
xi = m;
eta = l;
if (rhu < 1.0e-8) {
/* Small angle formula. */
t = (90.0 + theta)*D2R;
rhu = t*t/2.0;
}
x0 = 0.0;
y0 = -2;
}
if (xi == 0.0 && eta == 0.0) {
xf = 0.0;
yf = 0.0;
} else if (-xi >= fabs(eta)) {
psi = eta/xi;
chi = 1.0 + psi*psi;
xf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
yf = (xf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
} else if (xi >= fabs(eta)) {
psi = eta/xi;
chi = 1.0 + psi*psi;
xf = sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
yf = (xf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
} else if (-eta > fabs(xi)) {
psi = xi/eta;
chi = 1.0 + psi*psi;
yf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
xf = (yf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
} else {
psi = xi/eta;
chi = 1.0 + psi*psi;
yf = sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
xf = (yf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
}
if (fabs(xf) > 1.0) {
if (fabs(xf) > 1.0+tol) {
return 2;
}
xf = wcs_copysign(1.0,xf);
}
if (fabs(yf) > 1.0) {
if (fabs(yf) > 1.0+tol) {
return 2;
}
yf = wcs_copysign(1.0,yf);
}
*x = prj->w[0]*(xf + x0);
*y = prj->w[0]*(yf + y0);
return 0;
}
/*--------------------------------------------------------------------------*/
int qscrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
int direct, face;
double chi, l, m, n, psi, rho, rhu, xf, yf, w;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (qscset(prj)) return 1;
}
xf = x*prj->w[1];
yf = y*prj->w[1];
/* Check bounds. */
if (fabs(xf) <= 1.0) {
if (fabs(yf) > 3.0) return 2;
} else {
if (fabs(xf) > 7.0) return 2;
if (fabs(yf) > 1.0) return 2;
}
/* Map negative faces to the other side. */
if (xf < -1.0) xf += 8.0;
/* Determine the face. */
if (xf > 5.0) {
face = 4;
xf = xf - 6.0;
} else if (xf > 3.0) {
face = 3;
xf = xf - 4.0;
} else if (xf > 1.0) {
face = 2;
xf = xf - 2.0;
} else if (yf > 1.0) {
face = 0;
yf = yf - 2.0;
} else if (yf < -1.0) {
face = 5;
yf = yf + 2.0;
} else {
face = 1;
}
direct = (fabs(xf) > fabs(yf));
if (direct) {
if (xf == 0.0) {
psi = 0.0;
chi = 1.0;
rho = 1.0;
rhu = 0.0;
} else {
w = 15.0*yf/xf;
psi = wcs_sind(w)/(wcs_cosd(w) - SQRT2INV);
chi = 1.0 + psi*psi;
rhu = xf*xf*(1.0 - 1.0/sqrt(1.0 + chi));
rho = 1.0 - rhu;
}
} else {
if (yf == 0.0) {
psi = 0.0;
chi = 1.0;
rho = 1.0;
rhu = 0.0;
} else {
w = 15.0*xf/yf;
psi = wcs_sind(w)/(wcs_cosd(w) - SQRT2INV);
chi = 1.0 + psi*psi;
rhu = yf*yf*(1.0 - 1.0/sqrt(1.0 + chi));
rho = 1.0 - rhu;
}
}
if (rho < -1.0) {
if (rho < -1.0-tol) {
return 2;
}
rho = -1.0;
rhu = 2.0;
w = 0.0;
} else {
w = sqrt(rhu*(2.0-rhu)/chi);
}
if (face == 0) {
n = rho;
if (direct) {
m = w;
if (xf < 0.0) m = -m;
l = -m*psi;
} else {
l = w;
if (yf > 0.0) l = -l;
m = -l*psi;
}
} else if (face == 1) {
l = rho;
if (direct) {
m = w;
if (xf < 0.0) m = -m;
n = m*psi;
} else {
n = w;
if (yf < 0.0) n = -n;
m = n*psi;
}
} else if (face == 2) {
m = rho;
if (direct) {
l = w;
if (xf > 0.0) l = -l;
n = -l*psi;
} else {
n = w;
if (yf < 0.0) n = -n;
l = -n*psi;
}
} else if (face == 3) {
l = -rho;
if (direct) {
m = w;
if (xf > 0.0) m = -m;
n = -m*psi;
} else {
n = w;
if (yf < 0.0) n = -n;
m = -n*psi;
}
} else if (face == 4) {
m = -rho;
if (direct) {
l = w;
if (xf < 0.0) l = -l;
n = l*psi;
} else {
n = w;
if (yf < 0.0) n = -n;
l = n*psi;
}
} else {
n = -rho;
if (direct) {
m = w;
if (xf < 0.0) m = -m;
l = m*psi;
} else {
l = w;
if (yf < 0.0) l = -l;
m = l*psi;
}
}
if (l == 0.0 && m == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(m, l);
}
*theta = wcs_asind(n);
return 0;
}
/*============================================================================
* TSC: tangential spherical cube projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
* prj->w[0] r0*(pi/4)
* prj->w[1] (4/pi)/r0
*===========================================================================*/
int tscset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) {
prj->r0 = R2D;
prj->w[0] = 45.0;
prj->w[1] = 1.0/45.0;
} else {
prj->w[0] = prj->r0*PI/4.0;
prj->w[1] = 1.0/prj->w[0];
}
prj->flag = PRJSET;
return 0;
}
/*--------------------------------------------------------------------------*/
int tscfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
int face;
double costhe, l, m, n, rho, x0, xf, y0, yf;
const double tol = 1.0e-12;
if (prj->flag != PRJSET) {
if (tscset(prj)) return 1;
}
costhe = wcs_cosd(theta);
l = costhe*wcs_cosd(phi);
m = costhe*wcs_sind(phi);
n = wcs_sind(theta);
face = 0;
rho = n;
if (l > rho) {
face = 1;
rho = l;
}
if (m > rho) {
face = 2;
rho = m;
}
if (-l > rho) {
face = 3;
rho = -l;
}
if (-m > rho) {
face = 4;
rho = -m;
}
if (-n > rho) {
face = 5;
rho = -n;
}
if (face == 0) {
xf = m/rho;
yf = -l/rho;
x0 = 0.0;
y0 = 2.0;
} else if (face == 1) {
xf = m/rho;
yf = n/rho;
x0 = 0.0;
y0 = 0.0;
} else if (face == 2) {
xf = -l/rho;
yf = n/rho;
x0 = 2.0;
y0 = 0.0;
} else if (face == 3) {
xf = -m/rho;
yf = n/rho;
x0 = 4.0;
y0 = 0.0;
} else if (face == 4) {
xf = l/rho;
yf = n/rho;
x0 = 6.0;
y0 = 0.0;
} else {
xf = m/rho;
yf = l/rho;
x0 = 0.0;
y0 = -2.0;
}
if (fabs(xf) > 1.0) {
if (fabs(xf) > 1.0+tol) {
return 2;
}
xf = wcs_copysign(1.0,xf);
}
if (fabs(yf) > 1.0) {
if (fabs(yf) > 1.0+tol) {
return 2;
}
yf = wcs_copysign(1.0,yf);
}
*x = prj->w[0]*(xf + x0);
*y = prj->w[0]*(yf + y0);
return 0;
}
/*--------------------------------------------------------------------------*/
int tscrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double l, m, n, xf, yf;
if (prj->flag != PRJSET) {
if (tscset(prj)) return 1;
}
xf = x*prj->w[1];
yf = y*prj->w[1];
/* Check bounds. */
if (fabs(xf) <= 1.0) {
if (fabs(yf) > 3.0) return 2;
} else {
if (fabs(xf) > 7.0) return 2;
if (fabs(yf) > 1.0) return 2;
}
/* Map negative faces to the other side. */
if (xf < -1.0) xf += 8.0;
/* Determine the face. */
if (xf > 5.0) {
/* face = 4 */
xf = xf - 6.0;
m = -1.0/sqrt(1.0 + xf*xf + yf*yf);
l = -m*xf;
n = -m*yf;
} else if (xf > 3.0) {
/* face = 3 */
xf = xf - 4.0;
l = -1.0/sqrt(1.0 + xf*xf + yf*yf);
m = l*xf;
n = -l*yf;
} else if (xf > 1.0) {
/* face = 2 */
xf = xf - 2.0;
m = 1.0/sqrt(1.0 + xf*xf + yf*yf);
l = -m*xf;
n = m*yf;
} else if (yf > 1.0) {
/* face = 0 */
yf = yf - 2.0;
n = 1.0/sqrt(1.0 + xf*xf + yf*yf);
l = -n*yf;
m = n*xf;
} else if (yf < -1.0) {
/* face = 5 */
yf = yf + 2.0;
n = -1.0/sqrt(1.0 + xf*xf + yf*yf);
l = -n*yf;
m = -n*xf;
} else {
/* face = 1 */
l = 1.0/sqrt(1.0 + xf*xf + yf*yf);
m = l*xf;
n = l*yf;
}
if (l == 0.0 && m == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(m, l);
}
*theta = wcs_asind(n);
return 0;
}
/*============================================================================
* TNX: IRAF's gnomonic projection.
*
* Given and/or returned:
* prj->r0 r0; reset to 180/pi if 0.
*===========================================================================*/
int tnxset(prj)
struct prjprm *prj;
{
if (prj->r0 == 0.0) prj->r0 = R2D;
if (prj->flag == -1) {
prj->flag = -PRJSET;
} else {
prj->flag = PRJSET;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int tnxfwd(phi, theta, prj, x, y)
const double phi, theta;
struct prjprm *prj;
double *x, *y;
{
double r, s, xp[2];
if (abs(prj->flag) != PRJSET) {
if(tnxset(prj)) return 1;
}
s = wcs_sind(theta);
if (s == 0.0) return 2;
r = prj->r0*wcs_cosd(theta)/s;
xp[0] = r*wcs_sind(phi);
xp[1] = -r*wcs_cosd(phi);
*x = prj->inv_x? poly_func(prj->inv_x, xp) : xp[0];
*y = prj->inv_y? poly_func(prj->inv_y, xp) : xp[1];
if (prj->flag == PRJSET && s < 0.0) {
return 2;
}
return 0;
}
/*--------------------------------------------------------------------------*/
int tnxrev(x, y, prj, phi, theta)
const double x, y;
struct prjprm *prj;
double *phi, *theta;
{
double rp,xp,yp;
if (abs(prj->flag) != PRJSET) {
if (tanset(prj)) return 1;
}
xp = x+raw_to_tnxaxis(prj->tnx_lngcor, x, y);
yp = y+raw_to_tnxaxis(prj->tnx_latcor, x, y);
if ((rp = sqrt(xp*xp+yp*yp)) == 0.0) {
*phi = 0.0;
} else {
*phi = wcs_atan2d(xp, -yp);
}
*theta = wcs_atan2d(prj->r0, rp);
return 0;
}
/*--------------------------------------------------------------------------*/
int raw_to_pv(struct prjprm *prj, double x, double y, double *xo, double *yo)
{
int k;
double *a,*b,
r,r3,r5,r7,xy,x2,x3,x4,x5,x6,x7,y2,y3,y4,y5,y6,y7,xp,yp;
if (abs(prj->flag) != PRJSET) {
if (tanset(prj)) return 1;
}
k=prj->n;
a = prj->p; /* Longitude */
b = prj->p+100; /* Latitude */
xp = *(a++);
xp += *(a++)*x;
yp = *(b++);
yp += *(b++)*y;
if (!--k) goto poly_end;
xp += *(a++)*y;
yp += *(b++)*x;
if (!--k) goto poly_end;
r = sqrt(x*x + y*y);
xp += *(a++)*r;
yp += *(b++)*r;
if (!--k) goto poly_end;
xp += *(a++)*(x2=x*x);
yp += *(b++)*(y2=y*y);
if (!--k) goto poly_end;
xp += *(a++)*(xy=x*y);
yp += *(b++)*xy;
if (!--k) goto poly_end;
xp += *(a++)*y2;
yp += *(b++)*x2;
if (!--k) goto poly_end;
xp += *(a++)*(x3=x*x2);
yp += *(b++)*(y3=y*y2);
if (!--k) goto poly_end;
xp += *(a++)*x2*y;
yp += *(b++)*y2*x;
if (!--k) goto poly_end;
xp += *(a++)*x*y2;
yp += *(b++)*y*x2;
if (!--k) goto poly_end;
xp += *(a++)*y3;
yp += *(b++)*x3;
if (!--k) goto poly_end;
xp += *(a++)*(r3=r*r*r);
yp += *(b++)*r3;
if (!--k) goto poly_end;
xp += *(a++)*(x4=x2*x2);
yp += *(b++)*(y4=y2*y2);
if (!--k) goto poly_end;
xp += *(a++)*x3*y;
yp += *(b++)*y3*x;
if (!--k) goto poly_end;
xp += *(a++)*x2*y2;
yp += *(b++)*x2*y2;
if (!--k) goto poly_end;
xp += *(a++)*x*y3;
yp += *(b++)*y*x3;
if (!--k) goto poly_end;
xp += *(a++)*y4;
yp += *(b++)*x4;
if (!--k) goto poly_end;
xp += *(a++)*(x5=x4*x);
yp += *(b++)*(y5=y4*y);
if (!--k) goto poly_end;
xp += *(a++)*x4*y;
yp += *(b++)*y4*x;
if (!--k) goto poly_end;
xp += *(a++)*x3*y2;
yp += *(b++)*y3*x2;
if (!--k) goto poly_end;
xp += *(a++)*x2*y3;
yp += *(b++)*y2*x3;
if (!--k) goto poly_end;
xp += *(a++)*x*y4;
yp += *(b++)*y*x4;
if (!--k) goto poly_end;
xp += *(a++)*y5;
yp += *(b++)*x5;
if (!--k) goto poly_end;
xp += *(a++)*(r5=r3*r*r);
yp += *(b++)*r5;
if (!--k) goto poly_end;
xp += *(a++)*(x6=x5*x);
yp += *(b++)*(y6=y5*y);
if (!--k) goto poly_end;
xp += *(a++)*x5*y;
yp += *(b++)*y5*x;
if (!--k) goto poly_end;
xp += *(a++)*x4*y2;
yp += *(b++)*y4*x2;
if (!--k) goto poly_end;
xp += *(a++)*x3*y3;
yp += *(b++)*y3*x3;
if (!--k) goto poly_end;
xp += *(a++)*x2*y4;
yp += *(b++)*y4*x2;
if (!--k) goto poly_end;
xp += *(a++)*x*y5;
yp += *(b++)*y*x5;
if (!--k) goto poly_end;
xp += *(a++)*y6;
yp += *(b++)*x6;
if (!--k) goto poly_end;
xp += *(a++)*(x7=x6*x);
yp += *(b++)*(y7=y6*y);
if (!--k) goto poly_end;
xp += *(a++)*x6*y;
yp += *(b++)*y6*x;
if (!--k) goto poly_end;
xp += *(a++)*x5*y2;
yp += *(b++)*y5*x2;
if (!--k) goto poly_end;
xp += *(a++)*x4*y3;
yp += *(b++)*y4*x3;
if (!--k) goto poly_end;
xp += *(a++)*x3*y4;
yp += *(b++)*y3*x4;
if (!--k) goto poly_end;
xp += *(a++)*x2*y5;
yp += *(b++)*y2*x5;
if (!--k) goto poly_end;
xp += *(a++)*x*y6;
yp += *(b++)*y*x6;
if (!--k) goto poly_end;
xp += *(a++)*y7;
yp += *(b++)*x7;
if (!--k) goto poly_end;
xp += *a*(r7=r5*r*r);
yp += *b*r7;
poly_end:
*xo = xp;
*yo = yp;
return 0;
}
src/wcs/proj.h 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
* Author: Mark Calabretta, Australia Telescope National Facility
* IRAF's TNX added by E.Bertin 2000/03/28
* $Id: proj.h,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifndef WCSLIB_PROJ
#define WCSLIB_PROJ
#ifdef __cplusplus
extern "C" {
#endif
struct prjprm {
int flag;
int n;
double r0;
double p[200];
double w[10];
struct tnxaxis *tnx_latcor;
struct tnxaxis *tnx_lngcor;
struct poly *inv_x;
struct poly *inv_y;
};
#if __STDC__ || defined(__cplusplus)
int azpset(struct prjprm *);
int azpfwd(const double, const double, struct prjprm *, double *, double *);
int azprev(const double, const double, struct prjprm *, double *, double *);
int tanset(struct prjprm *);
int tanfwd(const double, const double, struct prjprm *, double *, double *);
int tanrev(const double, const double, struct prjprm *, double *, double *);
int sinset(struct prjprm *);
int sinfwd(const double, const double, struct prjprm *, double *, double *);
int sinrev(const double, const double, struct prjprm *, double *, double *);
int stgset(struct prjprm *);
int stgfwd(const double, const double, struct prjprm *, double *, double *);
int stgrev(const double, const double, struct prjprm *, double *, double *);
int arcset(struct prjprm *);
int arcfwd(const double, const double, struct prjprm *, double *, double *);
int arcrev(const double, const double, struct prjprm *, double *, double *);
int zpnset(struct prjprm *);
int zpnfwd(const double, const double, struct prjprm *, double *, double *);
int zpnrev(const double, const double, struct prjprm *, double *, double *);
int zeaset(struct prjprm *);
int zeafwd(const double, const double, struct prjprm *, double *, double *);
int zearev(const double, const double, struct prjprm *, double *, double *);
int airset(struct prjprm *);
int airfwd(const double, const double, struct prjprm *, double *, double *);
int airrev(const double, const double, struct prjprm *, double *, double *);
int cypset(struct prjprm *);
int cypfwd(const double, const double, struct prjprm *, double *, double *);
int cyprev(const double, const double, struct prjprm *, double *, double *);
int carset(struct prjprm *);
int carfwd(const double, const double, struct prjprm *, double *, double *);
int carrev(const double, const double, struct prjprm *, double *, double *);
int merset(struct prjprm *);
int merfwd(const double, const double, struct prjprm *, double *, double *);
int merrev(const double, const double, struct prjprm *, double *, double *);
int ceaset(struct prjprm *);
int ceafwd(const double, const double, struct prjprm *, double *, double *);
int cearev(const double, const double, struct prjprm *, double *, double *);
int copset(struct prjprm *);
int copfwd(const double, const double, struct prjprm *, double *, double *);
int coprev(const double, const double, struct prjprm *, double *, double *);
int codset(struct prjprm *);
int codfwd(const double, const double, struct prjprm *, double *, double *);
int codrev(const double, const double, struct prjprm *, double *, double *);
int coeset(struct prjprm *);
int coefwd(const double, const double, struct prjprm *, double *, double *);
int coerev(const double, const double, struct prjprm *, double *, double *);
int cooset(struct prjprm *);
int coofwd(const double, const double, struct prjprm *, double *, double *);
int coorev(const double, const double, struct prjprm *, double *, double *);
int bonset(struct prjprm *);
int bonfwd(const double, const double, struct prjprm *, double *, double *);
int bonrev(const double, const double, struct prjprm *, double *, double *);
int pcoset(struct prjprm *);
int pcofwd(const double, const double, struct prjprm *, double *, double *);
int pcorev(const double, const double, struct prjprm *, double *, double *);
int glsset(struct prjprm *);
int glsfwd(const double, const double, struct prjprm *, double *, double *);
int glsrev(const double, const double, struct prjprm *, double *, double *);
int parset(struct prjprm *);
int parfwd(const double, const double, struct prjprm *, double *, double *);
int parrev(const double, const double, struct prjprm *, double *, double *);
int aitset(struct prjprm *);
int aitfwd(const double, const double, struct prjprm *, double *, double *);
int aitrev(const double, const double, struct prjprm *, double *, double *);
int molset(struct prjprm *);
int molfwd(const double, const double, struct prjprm *, double *, double *);
int molrev(const double, const double, struct prjprm *, double *, double *);
int cscset(struct prjprm *);
int cscfwd(const double, const double, struct prjprm *, double *, double *);
int cscrev(const double, const double, struct prjprm *, double *, double *);
int qscset(struct prjprm *);
int qscfwd(const double, const double, struct prjprm *, double *, double *);
int qscrev(const double, const double, struct prjprm *, double *, double *);
int tscset(struct prjprm *);
int tscfwd(const double, const double, struct prjprm *, double *, double *);
int tscrev(const double, const double, struct prjprm *, double *, double *);
int tnxset(struct prjprm *);
int tnxfwd(const double, const double, struct prjprm *, double *, double *);
int tnxrev(const double, const double, struct prjprm *, double *, double *);
int raw_to_pv(struct prjprm *, double, double, double *, double *);
#else
int azpset(), azpfwd(), azprev();
int tanset(), tanfwd(), tanrev();
int sinset(), sinfwd(), sinrev();
int stgset(), stgfwd(), stgrev();
int arcset(), arcfwd(), arcrev();
int zpnset(), zpnfwd(), zpnrev();
int zeaset(), zeafwd(), zearev();
int airset(), airfwd(), airrev();
int cypset(), cypfwd(), cyprev();
int carset(), carfwd(), carrev();
int merset(), merfwd(), merrev();
int ceaset(), ceafwd(), cearev();
int copset(), copfwd(), coprev();
int codset(), codfwd(), codrev();
int coeset(), coefwd(), coerev();
int cooset(), coofwd(), coorev();
int bonset(), bonfwd(), bonrev();
int pcoset(), pcofwd(), pcorev();
int glsset(), glsfwd(), glsrev();
int parset(), parfwd(), parrev();
int aitset(), aitfwd(), aitrev();
int molset(), molfwd(), molrev();
int cscset(), cscfwd(), cscrev();
int qscset(), qscfwd(), qscrev();
int tscset(), tscfwd(), tscrev();
int tnxset(), tnxfwd(), tnxrev();
#endif
/*
extern const char *prjset_errmsg[];
extern const char *prjfwd_errmsg[];
extern const char *prjrev_errmsg[];
*/
#define PRJSET 137
#ifdef __cplusplus
};
#endif
#endif /* WCSLIB_PROJ */
src/wcs/sph.c 0 → 100644
View file @ 3b348a94
/*============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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 for the spherical coordinate transformations used by the FITS
* "World Coordinate System" (WCS) convention.
*
* Summary of routines
* -------------------
* The spherical coordinate transformations are implemented via separate
* functions for the transformation in each direction.
*
* Forward transformation; sphfwd()
* --------------------------------
* Transform celestial coordinates to the native coordinates of a projection.
*
* Given:
* lng,lat double Celestial longitude and latitude, in degrees.
* eul[5] double Euler angles for the transformation:
* 0: Celestial longitude of the native pole, in
* degrees.
* 1: Celestial colatitude of the native pole, or
* native colatitude of the celestial pole, in
* degrees.
* 2: Native longitude of the celestial pole, in
* degrees.
* 3: cos(eul[1])
* 4: sin(eul[1])
*
* Returned:
* phi, double Longitude and latitude in the native coordinate
* theta system of the projection, in degrees.
*
* Function return value:
* int Error status
* 0: Success.
*
* Reverse transformation; sphrev()
* --------------------------------
* Transform native coordinates of a projection to celestial coordinates.
*
* Given:
* phi, double Longitude and latitude in the native coordinate
* theta system of the projection, in degrees.
* eul[5] double Euler angles for the transformation:
* 0: Celestial longitude of the native pole, in
* degrees.
* 1: Celestial colatitude of the native pole, or
* native colatitude of the celestial pole, in
* degrees.
* 2: Native longitude of the celestial pole, in
* degrees.
* 3: cos(eul[1])
* 4: sin(eul[1])
*
* Returned:
* lng,lat double Celestial longitude and latitude, in degrees.
*
* Function return value:
* int Error status
* 0: Success.
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: sph.c,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include "wcstrig.h"
#include "sph.h"
#ifndef __STDC__
#ifndef const
#define const
#endif
#endif
#define wcs_copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
const double tol = 1.0e-5;
int sphfwd (lng, lat, eul, phi, theta)
const double lat, lng, eul[5];
double *phi, *theta;
{
double coslat, coslng, dlng, dphi, sinlat, sinlng, x, y, z;
coslat = wcs_cosd(lat);
sinlat = wcs_sind(lat);
dlng = lng - eul[0];
coslng = wcs_cosd(dlng);
sinlng = wcs_sind(dlng);
/* Compute the native longitude. */
x = sinlat*eul[4] - coslat*eul[3]*coslng;
if (fabs(x) < tol) {
/* Rearrange formula to reduce roundoff errors. */
x = -wcs_cosd(lat+eul[1]) + coslat*eul[3]*(1.0 - coslng);
}
y = -coslat*sinlng;
if (x != 0.0 || y != 0.0) {
dphi = wcs_atan2d(y, x);
} else {
/* Change of origin of longitude. */
dphi = dlng - 180.0;
}
*phi = eul[2] + dphi;
/* Normalize the native longitude. */
if (*phi > 180.0) {
*phi -= 360.0;
} else if (*phi < -180.0) {
*phi += 360.0;
}
/* Compute the native latitude. */
if (fmod(dlng,180.0) == 0.0) {
*theta = lat + coslng*eul[1];
if (*theta > 90.0) *theta = 180.0 - *theta;
if (*theta < -90.0) *theta = -180.0 - *theta;
} else {
z = sinlat*eul[3] + coslat*eul[4]*coslng;
if (fabs(z) > 0.99) {
/* Use an alternative formula for greater numerical accuracy. */
*theta = wcs_copysign(wcs_acosd(sqrt(x*x+y*y)), z);
} else {
*theta = wcs_asind(z);
}
}
return 0;
}
/*-----------------------------------------------------------------------*/
int sphrev (phi, theta, eul, lng, lat)
const double phi, theta, eul[5];
double *lng, *lat;
{
double cosphi, costhe, dlng, dphi, sinphi, sinthe, x, y, z;
costhe = wcs_cosd(theta);
sinthe = wcs_sind(theta);
dphi = phi - eul[2];
cosphi = wcs_cosd(dphi);
sinphi = wcs_sind(dphi);
/* Compute the celestial longitude. */
x = sinthe*eul[4] - costhe*eul[3]*cosphi;
if (fabs(x) < tol) {
/* Rearrange formula to reduce roundoff errors. */
x = -wcs_cosd(theta+eul[1]) + costhe*eul[3]*(1.0 - cosphi);
}
y = -costhe*sinphi;
if (x != 0.0 || y != 0.0) {
dlng = wcs_atan2d(y, x);
} else {
/* Change of origin of longitude. */
dlng = dphi + 180.0;
}
*lng = eul[0] + dlng;
/* Normalize the celestial longitude. */
if (eul[0] >= 0.0) {
if (*lng < 0.0) *lng += 360.0;
} else {
if (*lng > 0.0) *lng -= 360.0;
}
if (*lng > 360.0) {
*lng -= 360.0;
} else if (*lng < -360.0) {
*lng += 360.0;
}
/* Compute the celestial latitude. */
if (fmod(dphi,180.0) == 0.0) {
*lat = theta + cosphi*eul[1];
if (*lat > 90.0) *lat = 180.0 - *lat;
if (*lat < -90.0) *lat = -180.0 - *lat;
} else {
z = sinthe*eul[3] + costhe*eul[4]*cosphi;
if (fabs(z) > 0.99) {
/* Use an alternative formula for greater numerical accuracy. */
*lat = wcs_copysign(wcs_acosd(sqrt(x*x+y*y)), z);
} else {
*lat = wcs_asind(z);
}
}
return 0;
}
src/wcs/sph.h 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: sph.h,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifndef WCSLIB_SPH
#define WCSLIB_SPH
#ifdef __cplusplus
extern "C" {
#endif
#if __STDC__ || defined(__cplusplus)
int sphfwd(const double, const double,
const double [],
double *, double *);
int sphrev(const double, const double,
const double [],
double *, double *);
#else
int sphfwd(), sphrev();
#endif
#ifdef __cplusplus
}
#endif
#endif /* WCSLIB_SPH */
src/wcs/tnx.c 0 → 100644
View file @ 3b348a94
/*
tnx.c
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: WCSlib
*
* Author: E.BERTIN (IAP), based on D.Mink (SAO) WCSTools
*
* Contents: Handle TNX astrometric format (from IRAF).
*
*
* Last modify: 04/07/2006
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "tnx.h"
/******* read_tnxaxis *********************************************************
PROTO tnxaxisstruct *read_tnxaxis(char *tnxstr)
PURPOSE Read a TNX axis mapping structure.
INPUT String containing the TNX info.
OUTPUT TNXAXIS structure if OK, or NULL in case of error.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 04/07/2006
***/
tnxaxisstruct *read_tnxaxis(char *tnxstr)
{
tnxaxisstruct *tnxaxis;
char *pstr, *ptr;
double min, max;
int i, order;
if ((pstr=strpbrk(tnxstr, "1234567890-+.")))
{
if (!(tnxaxis=malloc(sizeof(tnxaxisstruct))))
return NULL;
tnxaxis->type = (int)(atof(strtok_r(pstr, " ", &ptr))+0.5);
tnxaxis->xorder = (pstr=strtok_r(NULL, " ", &ptr))?
(int)(atof(pstr)+0.5) : 0;
tnxaxis->yorder = (pstr=strtok_r(NULL, " ", &ptr))?
(int)(atof(pstr)+0.5) : 0;
tnxaxis->xterms = (pstr=strtok_r(NULL, " ", &ptr))?
(int)(atof(pstr)+0.5) : 0;
min = (pstr=strtok_r(NULL, " ", &ptr))? atof(pstr) : 0.0;
max = (pstr=strtok_r(NULL, " ", &ptr))? atof(pstr) : 0.0;
if (max <= min)
return NULL;
tnxaxis->xrange = 2.0 / (max - min);
tnxaxis->xmaxmin = - (max + min) / 2.0;
min = (pstr=strtok_r(NULL, " ", &ptr))? atof(pstr) : 0.0;
max = (pstr=strtok_r(NULL, " ", &ptr))? atof(pstr) : 0.0;
if (max <= min)
return NULL;
tnxaxis->yrange = 2.0 / (max - min);
tnxaxis->ymaxmin = - (max + min) / 2.0;
switch (tnxaxis->xterms)
{
case TNX_XNONE:
tnxaxis->ncoeff = tnxaxis->xorder + tnxaxis->yorder - 1;
break;
case TNX_XHALF:
order = tnxaxis->xorder<tnxaxis->yorder?
tnxaxis->xorder : tnxaxis->yorder;
tnxaxis->ncoeff = tnxaxis->xorder*tnxaxis->yorder - order*(order-1)/2;
break;
case TNX_XFULL:
tnxaxis->ncoeff = tnxaxis->xorder * tnxaxis->yorder;
break;
default:
return NULL;
}
/*-- Now read the mapping coefficients */
if (!(tnxaxis->coeff=malloc(tnxaxis->ncoeff*sizeof(double))))
return NULL;
for (i=0; i<tnxaxis->ncoeff && (pstr=strtok_r(NULL, " ", &ptr)); i++)
tnxaxis->coeff[i] = atof(pstr);
if (i!=tnxaxis->ncoeff)
return NULL;
if (!(tnxaxis->xbasis=malloc(tnxaxis->xorder*sizeof(double))))
return NULL;
if (!(tnxaxis->ybasis=malloc(tnxaxis->yorder*sizeof(double))))
return NULL;
return tnxaxis;
}
else
return NULL;
}
/******* copy_tnxaxis *********************************************************
PROTO tnxaxisstruct *copy_tnxaxis(tnxaxisstruct *axis)
PURPOSE Copy a TNX axis mapping structure.
INPUT TNXAXIS structure pointer.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 28/11/2003
***/
tnxaxisstruct *copy_tnxaxis(tnxaxisstruct *axis)
{
tnxaxisstruct *tnxaxis;
int i;
if (axis)
{
if (!axis->ncoeff)
return NULL;
if (!(tnxaxis=malloc(sizeof(tnxaxisstruct))))
return NULL;
*tnxaxis = *axis;
if (!(tnxaxis->coeff=malloc(tnxaxis->ncoeff*sizeof(double))))
return NULL;
for (i=0; i<tnxaxis->ncoeff; i++)
tnxaxis->coeff[i] = axis->coeff[i];
if (!(tnxaxis->xbasis=malloc(tnxaxis->xorder*sizeof(double))))
return NULL;
for (i=0; i<tnxaxis->xorder; i++)
tnxaxis->xbasis[i] = axis->xbasis[i];
if (!(tnxaxis->ybasis=malloc(tnxaxis->yorder*sizeof(double))))
return NULL;
for (i=0; i<tnxaxis->yorder; i++)
tnxaxis->ybasis[i] = axis->ybasis[i];
return tnxaxis;
}
return NULL;
}
/******* free_tnxaxis *********************************************************
PROTO void free_tnxaxis(tnxaxisstruct *axis)
PURPOSE Free a TNX axis mapping structure.
INPUT TNXAXIS structure pointer.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 09/04/2000
***/
void free_tnxaxis(tnxaxisstruct *axis)
{
if (axis)
{
free(axis->coeff);
free(axis->xbasis);
free(axis->ybasis);
free(axis);
}
return;
}
/******* raw_to_tnxaxis *******************************************************
PROTO double raw_to_tnxaxis(tnxaxisstruct *axis, double x, double y)
PURPOSE Compute the correction value on a TNX axis at current position.
INPUT TNXAXIS structure pointer,
x coordinate,
y coordinate.
OUTPUT Value on the TNXaxis.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 11/04/2000
***/
double raw_to_tnxaxis(tnxaxisstruct *axis, double x, double y)
{
double *xbasis, *ybasis,*coeff,
norm, accum, val;
int i, j, xorder,xorder0,yorder,maxorder,xterms;
xbasis = axis->xbasis;
ybasis = axis->ybasis;
xorder = axis->xorder;
yorder = axis->yorder;
xterms = axis->xterms;
switch (axis->type)
{
case TNX_CHEBYSHEV:
xbasis[0] = 1.0;
if (xorder > 1)
{
xbasis[1] = norm = (x + axis->xmaxmin)*axis->xrange;
if (xorder > 2)
for (i = 2; i < xorder; i++)
xbasis[i] = 2.0*norm*xbasis[i-1] - xbasis[i-2];
}
ybasis[0] = 1.0;
if (yorder > 1)
{
ybasis[1] = norm = (y + axis->ymaxmin)*axis->yrange;
if (yorder > 2)
for (i = 2; i < yorder; i++)
ybasis[i] = 2.0*norm*xbasis[i-1] - ybasis[i-2];
}
break;
case TNX_LEGENDRE:
xbasis[0] = 1.0;
if (xorder > 1)
{
xbasis[1] = norm = (x + axis->xmaxmin)*axis->xrange;
if (xorder > 2)
for (i = 2; (j=i) < xorder; i++)
xbasis[i] = ((2.0*j - 3.0) * norm * xbasis[i-1] -
(j - 2.0) * xbasis[i-2]) / (j - 1.0);
}
ybasis[0] = 1.0;
if (yorder > 1)
{
ybasis[1] = norm = (y + axis->ymaxmin)*axis->yrange;
if (yorder > 2)
for (i = 2; (j=i) < xorder; i++)
ybasis[i] = ((2.0*j - 3.0) * norm * ybasis[i-1] -
(j - 2.0) * ybasis[i-2]) / (j - 1.0);
}
break;
case TNX_POLYNOMIAL:
xbasis[0] = 1.0;
if (xorder > 1)
{
xbasis[1] = x;
if (xorder > 2)
for (i = 2; i < xorder; i++)
xbasis[i] = x * xbasis[i-1];
}
ybasis[0] = 1.0;
if (yorder > 1)
{
ybasis[1] = y;
if (yorder > 2)
for (i = 2; i < yorder; i++)
ybasis[i] = y * ybasis[i-1];
}
break;
default:
return 0.0;
}
/* Loop over y basis functions */
maxorder = xorder > yorder ? xorder : yorder;
xorder0 = xorder;
coeff = axis->coeff;
val = 0.0;
for (i = 0; i<yorder; i++)
{
/*-- Loop over the x basis functions */
accum = 0.0;
xbasis = axis->xbasis;
for (j = xorder; j--;)
accum += *(coeff++) * *(xbasis++);
val += accum**(ybasis++);
/*-- Elements of the coefficient vector where neither k = 1 or i = 1
are not calculated if sf->xterms = no. */
if (xterms == TNX_XNONE)
xorder = 1;
else if (xterms == TNX_XHALF && (i + 1 + xorder0) > maxorder)
xorder--;
}
return val;
}
src/wcs/tnx.h 0 → 100644
View file @ 3b348a94
/*
tnx.h
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: WCSlib
*
* Author: E.BERTIN (IAP), based on D.Mink (SAO) WCSTools
*
* Contents: Include to handle TNX astrometric format (from IRAF).
*
*
* Last modify: 28/11/2003
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifndef _TNX_H_
#define _TNX_H_
/*-------------------------------- macros -----------------------------------*/
#define TNX_MAXCHARS 2048 /* Maximum FITS "WAT" string length */
/* TNX permitted types of surfaces */
#define TNX_CHEBYSHEV 1
#define TNX_LEGENDRE 2
#define TNX_POLYNOMIAL 3
/* TNX cross-terms flags */
#define TNX_XNONE 0 /* no x-terms (old no) */
#define TNX_XFULL 1 /* full x-terms (new yes) */
#define TNX_XHALF 2 /* half x-terms (new) */
/*----------------------------- Internal constants --------------------------*/
/*------------------------------- structures --------------------------------*/
typedef struct tnxaxis
{
int type; /* Projection correction type */
int xorder,yorder; /* Polynomial orders */
int xterms; /* Well... */
int ncoeff; /* Number of polynom coefficients */
double xrange,yrange; /* Coordinate ranges */
double xmaxmin,ymaxmin; /* Well... */
double *coeff; /* Polynom coefficients */
double *xbasis,*ybasis; /* Basis function values */
} tnxaxisstruct;
/*------------------------------- functions ---------------------------------*/
tnxaxisstruct *copy_tnxaxis(tnxaxisstruct *axis),
*read_tnxaxis(char *tnxstr);
double raw_to_tnxaxis(tnxaxisstruct *axis, double x, double y);
void free_tnxaxis(tnxaxisstruct *axis);
#endif
src/wcs/wcs.c 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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.
*
* 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.
*
* 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.
* 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.
*
* 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.
*
* 4) The quadcube projections (CSC, QSC, TSC) 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 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include "stdio.h"
#include "string.h"
#include "wcsmath.h"
#include "wcstrig.h"
#include "sph.h"
#include "wcs.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 wcs_signbit(X) ((X) < 0.0 ? 1 : 0)
int wcsset (naxis, ctype, wcs)
const int naxis;
const char ctype[][9];
struct wcsprm *wcs;
{
int j, k, *ndx;
char requir[9];
strcpy(wcs->pcode, "");
strcpy(requir, "");
wcs->lng = -1;
ndx = &wcs->lng; /* to satisfy gcc -Wall */
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) {
/* Allow NCP to pass (will be converted to SIN later). */
if (strncmp(&ctype[j][5], "NCP", 3)) 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 {
/* Unrecognized celestial type. */
return 1;
}
} else {
if (strncmp(ctype[j], requir, 8) != 0) {
/* Inconsistent projection types. */
return 1;
}
*ndx = j;
strcpy(requir, "");
}
}
if (strcmp(requir, "")) {
/* Unmatched celestial axis. */
return 1;
}
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[2] == 0.0) {
return 2;
}
strcpy(wcs->pcode, "SIN");
prj->p[1] = 0.0;
prj->p[2] = wcs_cosd(cel->ref[2])/wcs_sind(cel->ref[2]);
prj->flag = 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[2] == 0.0) {
return 2;
}
strcpy(wcs->pcode, "SIN");
prj->p[1] = 0.0;
prj->p[2] = wcs_cosd(cel->ref[2])/wcs_sind(cel->ref[2]);
prj->flag = 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;
double lambda, span[2], step;
double pixmix;
double lng, lng0, lng0m, lng1, lng1m;
double lat, lat0, lat0m, lat1, lat1m;
double d, d0, d0m, d1, d1m, dx;
double dabs, dmin, lmin;
double phi0, phi1;
struct celprm cel0;
/* 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];
dx = 0.0; /* to satisfy gcc -Wall */
/* 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 (wcs_signbit(d0) != wcs_signbit(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 (wcs_signbit(d0) != wcs_signbit(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;
}
lat = lat0 + lambda*(lat1 - lat0);
world[wcs->lat] = lat;
if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
prj, imgcrd, lin, pixcrd))) {
return err;
}
d = pixcrd[mixpix] - pixmix;
/* Check for a solution. */
dabs = fabs(d);
if (dabs < tol) return 0;
/* Record the point of closest approach. */
if (dabs < dmin) {
lmin = lat;
dmin = dabs;
}
if (wcs_signbit(d0) == wcs_signbit(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 (wcs_signbit(d0) != wcs_signbit(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 (wcs_signbit(d0) != wcs_signbit(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;
}
lng = lng0 + lambda*(lng1 - lng0);
world[wcs->lng] = lng;
if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
prj, imgcrd, lin, pixcrd))) {
return err;
}
d = pixcrd[mixpix] - pixmix;
/* Check for a solution. */
dabs = fabs(d);
if (dabs < tol) return 0;
/* Record the point of closest approach. */
if (dabs < dmin) {
lmin = lng;
dmin = dabs;
}
if (wcs_signbit(d0) == wcs_signbit(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;
cel0.prjfwd = cel->prjfwd;
cel0.prjrev = cel->prjrev;
/* 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 = (float) 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 (wcs_signbit(d0) != wcs_signbit(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;
}
world[wcs->lng] = phi0 + lambda*(phi1 - phi0);
if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj,
imgcrd, lin, pixcrd))) {
return err;
}
d = pixcrd[mixpix] - pixmix;
/* Check for a solution. */
dabs = fabs(d);
if (dabs < tol) {
/* Recall saved world coordinates. */
world[wcs->lng] = lng;
world[wcs->lat] = lat;
return 0;
}
if (wcs_signbit(d0) == wcs_signbit(d)) {
phi0 = world[wcs->lng];
d0 = d;
} else {
phi1 = world[wcs->lng];
d1 = d;
}
}
}
/* No solution. */
return 5;
}
src/wcs/wcs.h 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: wcs.h,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifndef WCSLIB_WCS
#define WCSLIB_WCS
#include "cel.h"
#include "lin.h"
#ifdef __cplusplus
extern "C" {
#endif
struct wcsprm {
int flag;
char pcode[4];
char lngtyp[5], lattyp[5];
int lng, lat;
int cubeface;
};
#if __STDC__ || defined(__cplusplus)
int wcsset(const int,
const char[][9],
struct wcsprm *);
int wcsfwd(const char[][9],
struct wcsprm *,
const double[],
const double[],
struct celprm *,
double *,
double *,
struct prjprm *,
double[],
struct linprm *,
double[]);
int wcsrev(const char[][9],
struct wcsprm *,
const double[],
struct linprm *,
double[],
struct prjprm *,
double *,
double *,
const double[],
struct celprm *,
double[]);
int wcsmix(const char[][9],
struct wcsprm *,
const int,
const int,
const double[],
const double,
int,
double[],
const double[],
struct celprm *,
double *,
double *,
struct prjprm *,
double[],
struct linprm *,
double[]);
#else
int wcsset(), wcsfwd(), wcsrev(), wcsmix();
#endif
extern const char *wcsset_errmsg[];
extern const char *wcsfwd_errmsg[];
extern const char *wcsrev_errmsg[];
extern const char *wcsmix_errmsg[];
#define WCSSET 137
#ifdef __cplusplus
};
#endif
#endif /* WCSLIB_WCS */
src/wcs/wcsmath.h 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: wcsmath.h,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifndef WCSLIB_MATH
#define WCSLIB_MATH
#ifdef PI
#undef PI
#endif
#ifdef D2R
#undef D2R
#endif
#ifdef R2D
#undef R2D
#endif
#ifdef SQRT2
#undef SQRT2
#endif
#ifdef SQRT2INV
#undef SQRT2INV
#endif
#define PI 3.141592653589793238462643
#define D2R (PI/180.0)
#define R2D (180.0/PI)
#define SQRT2 1.4142135623730950488
#define SQRT2INV (1.0/SQRT2)
#endif /* WCSLIB_MATH */
src/wcs/wcstrig.c 0 → 100644
View file @ 3b348a94
/*============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
*=============================================================================
*
* The functions defined herein are trigonometric or inverse trigonometric
* functions which take or return angular arguments in decimal degrees.
*
* $Id: wcstrig.c,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*---------------------------------------------------------------------------*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include "wcstrig.h"
#ifndef PI /* EB 02/06/97 */
#define PI 3.141592653589793238462643
#endif /* EB 02/06/97 */
const double d2r = PI / 180.0;
const double r2d = 180.0 / PI;
#ifndef HAVE_MATHIMF_H
double wcs_cosd(angle)
const double angle;
{
double resid;
resid = fabs(fmod(angle,360.0));
if (resid == 0.0) {
return 1.0;
} else if (resid == 90.0) {
return 0.0;
} else if (resid == 180.0) {
return -1.0;
} else if (resid == 270.0) {
return 0.0;
}
return cos(angle*d2r);
}
/*--------------------------------------------------------------------------*/
double wcs_sind(angle)
const double angle;
{
double resid;
resid = fmod(angle-90.0,360.0);
if (resid == 0.0) {
return 1.0;
} else if (resid == 90.0) {
return 0.0;
} else if (resid == 180.0) {
return -1.0;
} else if (resid == 270.0) {
return 0.0;
}
return sin(angle*d2r);
}
/*--------------------------------------------------------------------------*/
double wcs_tand(angle)
const double angle;
{
double resid;
resid = fmod(angle,360.0);
if (resid == 0.0 || fabs(resid) == 180.0) {
return 0.0;
} else if (resid == 45.0 || resid == 225.0) {
return 1.0;
} else if (resid == -135.0 || resid == -315.0) {
return -1.0;
}
return tan(angle*d2r);
}
/*--------------------------------------------------------------------------*/
double wcs_acosd(v)
const double v;
{
if (v >= 1.0) {
if (v-1.0 < WCSTRIG_TOL) return 0.0;
} else if (v == 0.0) {
return 90.0;
} else if (v <= -1.0) {
if (v+1.0 > -WCSTRIG_TOL) return 180.0;
}
return acos(v)*r2d;
}
/*--------------------------------------------------------------------------*/
double wcs_asind(v)
const double v;
{
if (v <= -1.0) {
if (v+1.0 > -WCSTRIG_TOL) return -90.0;
} else if (v == 0.0) {
return 0.0;
} else if (v >= 1.0) {
if (v-1.0 < WCSTRIG_TOL) return 90.0;
}
return asin(v)*r2d;
}
/*--------------------------------------------------------------------------*/
double wcs_atand(v)
const double v;
{
if (v == -1.0) {
return -45.0;
} else if (v == 0.0) {
return 0.0;
} else if (v == 1.0) {
return 45.0;
}
return atan(v)*r2d;
}
/*--------------------------------------------------------------------------*/
double wcs_atan2d(y, x)
const double x, y;
{
if (y == 0.0) {
if (x >= 0.0) {
return 0.0;
} else if (x < 0.0) {
return 180.0;
}
} else if (x == 0.0) {
if (y > 0.0) {
return 90.0;
} else if (y < 0.0) {
return -90.0;
}
}
return atan2(y,x)*r2d;
}
#endif
src/wcs/wcstrig.h 0 → 100644
View file @ 3b348a94
/*=============================================================================
*
* WCSLIB - an implementation of the FITS WCS proposal.
* Copyright (C) 1995-1999, Mark Calabretta
*
* This library is free software; you can redistribute it and/or modify it
* under the terms of the GNU Library 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 Library
* General Public License for more details.
*
* You should have received a copy of the GNU Library General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 675 Mass Ave, Cambridge, MA 02139, 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
*
* Author: Mark Calabretta, Australia Telescope National Facility
* $Id: wcstrig.h,v 1.1.1.1 2002/03/15 16:33:26 bertin Exp $
*===========================================================================*/
#ifndef WCSLIB_TRIG
#define WCSLIB_TRIG
#ifdef __cplusplus
extern "C" {
#endif
#if !defined(__STDC__) && !defined(__cplusplus)
#ifndef const
#define const
#endif
#endif
#if __STDC__ || defined(__cplusplus)
#ifdef HAVE_MATHIMF_H
#define wcs_cosd cosd
#define wcs_sind sind
#define wcs_tand tand
#define wcs_acosd acosd
#define wcs_asind asind
#define wcs_atand atand
#define wcs_atan2d atan2d
#else
double wcs_cosd(const double);
double wcs_sind(const double);
double wcs_tand(const double);
double wcs_acosd(const double);
double wcs_asind(const double);
double wcs_atand(const double);
double wcs_atan2d(const double, const double);
#endif
#else
double wcs_cosd();
double wcs_sind();
double wcs_tand();
double wcs_acosd();
double wcs_asind();
double wcs_atand();
double wcs_atan2d();
#endif
/* Domain tolerance for asin and acos functions. */
#define WCSTRIG_TOL 1e-10
#endif /* TRIGD */
#ifdef __cplusplus
};
#endif /* WCSLIB_TRIG */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment