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

merged with SExFIGI branch

parent 5ae55cd0
Hide whitespace changes
Inline Side-by-side
src/fft.h 0 → 100644
View file @ ece261cd
/*
fft.h
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: A program that uses FFTs
*
* Author: E.BERTIN (IAP)
*
* Contents: Include for fft.c.
*
* Last modify: 29/11/2006
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifndef _FITSCAT_H_
#include "fits/fitscat.h"
#endif
/*---------------------------- Internal constants ---------------------------*/
/*------------------------------- Other Macros ------------------------------*/
#define QFFTWMALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)fftw_malloc((size_t)(nel)*sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QFFTWFREE(ptr) fftw_free(ptr)
/*--------------------------- structure definitions -------------------------*/
/*---------------------------------- protos --------------------------------*/
extern void fft_conv(double *data1, double *fdata2, int *size),
fft_end(),
fft_init();
extern double *fft_rtf(double *data, int *size);
src/field.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: Handling of field structures.
*
* Last modify: 29/06/2006
* Last modify: 19/12/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -32,6 +32,7 @@
#include "back.h"
#include "field.h"
#include "filter.h"
#include "fitswcs.h"
#include "interpolate.h"
/********************************* newfield **********************************/
......@@ -44,13 +45,16 @@ picstruct *newfield(char *filename, int flags, int nok)
picstruct *field;
catstruct *cat;
tabstruct *tab;
OFF_T mefpos = 0; /* To avoid gcc -Wall warnings */
int nok2, ntab, margin;
/* Move to nok'th valid FITS image extension */
if (!(cat = read_cat(filename)))
error(EXIT_FAILURE, "*Error*: cannot open ", filename);
close_cat(cat);
/* First allocate memory for the new field (and nullify pointers) */
QCALLOC(field, picstruct, 1);
field->flags = flags;
field->cat = cat;
tab = cat->tab;
nok++; /* At least one pass through the loop */
nok2 = nok;
......@@ -60,16 +64,12 @@ picstruct *newfield(char *filename, int flags, int nok)
|| !strncmp(tab->xtension, "BINTABLE", 8)
|| !strncmp(tab->xtension, "ASCTABLE", 8))
continue;
mefpos = tab->headpos;
field->tab = tab;
nok--;
}
if (ntab<0)
error(EXIT_FAILURE, "Not enough valid FITS image extensions in ",filename);
/* First allocate memory for the new field (and nullify pointers) */
QCALLOC(field, picstruct, 1);
field->mefpos = mefpos;
field->flags = flags;
strcpy (field->filename, filename);
/* A short, "relative" version of the filename */
if (!(field->rfilename = strrchr(field->filename, '/')))
......@@ -80,9 +80,12 @@ picstruct *newfield(char *filename, int flags, int nok)
sprintf(gstr, "Looking for %s", field->rfilename);
NFPRINTF(OUTPUT, gstr);
/* Check the image exists and read important info (image size, etc...) */
field->file = cat->file;
readimagehead(field);
if (cat->ntab>1)
sprintf(gstr, "[%d/%d]", nok2, cat->ntab-1);
sprintf(gstr, "[%d/%d]", nok2, cat->tab->naxis<2? cat->ntab-1 : cat->ntab);
QPRINTF(OUTPUT, "%s \"%.20s\" %s / %d x %d / %d bits %s data\n",
flags&FLAG_FIELD? "Flagging from:" :
(flags&(RMS_FIELD|VAR_FIELD|WEIGHT_FIELD)?
......@@ -93,12 +96,8 @@ picstruct *newfield(char *filename, int flags, int nok)
cat->ntab>1? gstr : "",
field->width, field->height, field->bytepix*8,
field->bitpix>0?
(field->compress_type!=ICOMPRESS_NONE?"COMPRESSED":"INTEGER")
:"FLOATING POINT");
/* Provide a buffer for compressed data */
if (field->compress_type != ICOMPRESS_NONE)
QMALLOC(field->compress_buf, char, FBSIZE);
(field->tab->compress_type!=COMPRESS_NONE?"COMPRESSED":"INTEGER")
:"FLOATING POINT");
/* Check the astrometric system and do the setup of the astrometric stuff */
if (prefs.world_flag && (flags & (MEASURE_FIELD|DETECT_FIELD)))
......@@ -106,6 +105,17 @@ picstruct *newfield(char *filename, int flags, int nok)
else
field->pixscale=prefs.pixel_scale;
/* Gain and Saturation */
if (flags & (DETECT_FIELD|MEASURE_FIELD))
{
if (fitsread(field->tab->headbuf, prefs.gain_key, &field->gain,
H_FLOAT, T_DOUBLE) != RETURN_OK)
field->gain = prefs.gain;
if (fitsread(field->tab->headbuf, prefs.satur_key, &field->satur_level,
H_FLOAT, T_DOUBLE) !=RETURN_OK)
field->satur_level = prefs.satur_level;
}
/* Background */
if (flags & (DETECT_FIELD|MEASURE_FIELD|WEIGHT_FIELD|VAR_FIELD|RMS_FIELD))
{
......@@ -143,16 +153,12 @@ picstruct *newfield(char *filename, int flags, int nok)
if (prefs.filter_flag)
{
/*-- If filtering is on, one should consider the height of the conv. mask */
int margin;
if (field->stripheight < thefilter->convh)
field->stripheight = thefilter->convh;
if (field->stripmargin < (margin = (thefilter->convh-1)/2))
field->stripmargin = margin;
}
free_cat(&cat, 1);
return field;
}
......@@ -172,15 +178,13 @@ picstruct *inheritfield(picstruct *infield, int flags)
/* Copy what is important and reset the remaining */
*field = *infield;
field->flags = flags;
copyastrom(infield, field);
QMEMCPY(infield->fitshead, field->fitshead, char, infield->fitsheadsize);
if (infield->wcs)
field->wcs = copy_wcs(infield->wcs);
field->interp_flag = 0;
field->assoc = NULL;
field->strip = NULL;
field->fstrip = NULL;
field->reffield = infield;
field->compress_buf = NULL;
field->compress_type = ICOMPRESS_NONE;
field->file = NULL;
return field;
......@@ -194,15 +198,14 @@ Free and close everything related to a field structure.
void endfield(picstruct *field)
{
if (field->file)
fclose(field->file);
free(field->fitshead);
/* Free cat only if associated with an open file */
if (field->file)
free_cat(&field->cat, 1);
free(field->strip);
free(field->fstrip);
free(field->compress_buf);
if (field->astrom)
endastrom(field);
if (field->wcs)
end_wcs(field->wcs);
if (field->interp_flag)
end_interpolate(field);
endback(field);
......
src/fits/Makefile.in
View file @ ece261cd
# Makefile.in generated by automake 1.9.6 from Makefile.am.
# Makefile.in generated by automake 1.10.1 from Makefile.am.
# @configure_input@
# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
# 2003, 2004, 2005 Free Software Foundation, Inc.
# 2003, 2004, 2005, 2006, 2007, 2008 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.
......@@ -14,15 +14,11 @@
@SET_MAKE@
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
......@@ -34,10 +30,15 @@ POST_INSTALL = :
NORMAL_UNINSTALL = :
PRE_UNINSTALL = :
POST_UNINSTALL = :
build_triplet = @build@
host_triplet = @host@
subdir = src/fits
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/acx_prog_cc_optim.m4 \
am__aclocal_m4_deps = $(top_srcdir)/acx_atlas.m4 \
$(top_srcdir)/acx_fftw.m4 $(top_srcdir)/acx_prog_cc_optim.m4 \
$(top_srcdir)/acx_pthread.m4 \
$(top_srcdir)/acx_urbi_resolve_dir.m4 \
$(top_srcdir)/configure.ac
am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \
$(ACLOCAL_M4)
......@@ -45,7 +46,6 @@ mkinstalldirs = $(SHELL) $(top_srcdir)/autoconf/mkinstalldirs
CONFIG_HEADER = $(top_builddir)/config.h
CONFIG_CLEAN_FILES =
LIBRARIES = $(noinst_LIBRARIES)
AR = ar
ARFLAGS = cru
libfits_a_AR = $(AR) $(ARFLAGS)
libfits_a_LIBADD =
......@@ -55,22 +55,30 @@ am_libfits_a_OBJECTS = fitsbody.$(OBJEXT) fitscat.$(OBJEXT) \
fitsread.$(OBJEXT) fitstab.$(OBJEXT) fitsutil.$(OBJEXT) \
fitswrite.$(OBJEXT)
libfits_a_OBJECTS = $(am_libfits_a_OBJECTS)
DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir)
DEFAULT_INCLUDES = -I.@am__isrc@ -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)
LTCOMPILE = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
--mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) \
$(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS)
CCLD = $(CC)
LINK = $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@
LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
--mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
$(LDFLAGS) -o $@
SOURCES = $(libfits_a_SOURCES)
DIST_SOURCES = $(libfits_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@
AR = @AR@
ATLAS_CFLAGS = @ATLAS_CFLAGS@
ATLAS_ERROR = @ATLAS_ERROR@
ATLAS_LIB = @ATLAS_LIB@
ATLAS_LIBPATH = @ATLAS_LIBPATH@
AUTOCONF = @AUTOCONF@
AUTOHEADER = @AUTOHEADER@
AUTOMAKE = @AUTOMAKE@
......@@ -80,16 +88,27 @@ CCDEPMODE = @CCDEPMODE@
CFLAGS = @CFLAGS@
CPP = @CPP@
CPPFLAGS = @CPPFLAGS@
CXX = @CXX@
CXXCPP = @CXXCPP@
CXXDEPMODE = @CXXDEPMODE@
CXXFLAGS = @CXXFLAGS@
CYGPATH_W = @CYGPATH_W@
DATE2 = @DATE2@
DATE3 = @DATE3@
DEFS = @DEFS@
DEPDIR = @DEPDIR@
ECHO = @ECHO@
ECHO_C = @ECHO_C@
ECHO_N = @ECHO_N@
ECHO_T = @ECHO_T@
EGREP = @EGREP@
EXEEXT = @EXEEXT@
F77 = @F77@
FFLAGS = @FFLAGS@
FFTW_ERROR = @FFTW_ERROR@
FFTW_LIBS = @FFTW_LIBS@
GREP = @GREP@
INSTALL = @INSTALL@
INSTALL_DATA = @INSTALL_DATA@
INSTALL_PROGRAM = @INSTALL_PROGRAM@
INSTALL_SCRIPT = @INSTALL_SCRIPT@
......@@ -97,8 +116,11 @@ INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@
LDFLAGS = @LDFLAGS@
LIBOBJS = @LIBOBJS@
LIBS = @LIBS@
LIBTOOL = @LIBTOOL@
LN_S = @LN_S@
LTLIBOBJS = @LTLIBOBJS@
MAKEINFO = @MAKEINFO@
MKDIR_P = @MKDIR_P@
OBJEXT = @OBJEXT@
PACKAGE = @PACKAGE@
PACKAGER = @PACKAGER@
......@@ -108,41 +130,66 @@ PACKAGE_STRING = @PACKAGE_STRING@
PACKAGE_TARNAME = @PACKAGE_TARNAME@
PACKAGE_VERSION = @PACKAGE_VERSION@
PATH_SEPARATOR = @PATH_SEPARATOR@
PTHREAD_CC = @PTHREAD_CC@
PTHREAD_CFLAGS = @PTHREAD_CFLAGS@
PTHREAD_LIBS = @PTHREAD_LIBS@
RANLIB = @RANLIB@
SED = @SED@
SET_MAKE = @SET_MAKE@
SHELL = @SHELL@
STRIP = @STRIP@
VERSION = @VERSION@
abs_builddir = @abs_builddir@
abs_srcdir = @abs_srcdir@
abs_top_builddir = @abs_top_builddir@
abs_top_srcdir = @abs_top_srcdir@
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@
ac_ct_CXX = @ac_ct_CXX@
ac_ct_F77 = @ac_ct_F77@
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 = @build@
build_alias = @build_alias@
build_cpu = @build_cpu@
build_os = @build_os@
build_vendor = @build_vendor@
builddir = @builddir@
datadir = @datadir@
datarootdir = @datarootdir@
docdir = @docdir@
dvidir = @dvidir@
exec_prefix = @exec_prefix@
host = @host@
host_alias = @host_alias@
host_cpu = @host_cpu@
host_os = @host_os@
host_vendor = @host_vendor@
htmldir = @htmldir@
includedir = @includedir@
infodir = @infodir@
install_sh = @install_sh@
libdir = @libdir@
libexecdir = @libexecdir@
localedir = @localedir@
localstatedir = @localstatedir@
mandir = @mandir@
mkdir_p = @mkdir_p@
oldincludedir = @oldincludedir@
pdfdir = @pdfdir@
prefix = @prefix@
program_transform_name = @program_transform_name@
psdir = @psdir@
sbindir = @sbindir@
sharedstatedir = @sharedstatedir@
srcdir = @srcdir@
sysconfdir = @sysconfdir@
target_alias = @target_alias@
top_builddir = @top_builddir@
top_srcdir = @top_srcdir@
# Program Makefile for the FITScat library
# Copyright (C) 2002 Emmanuel Bertin.
......@@ -155,7 +202,7 @@ libfits_a_SOURCES = fitsbody.c fitscat.c fitscheck.c fitscleanup.c \
all: all-am
.SUFFIXES:
.SUFFIXES: .c .o .obj
.SUFFIXES: .c .lo .o .obj
$(srcdir)/Makefile.in: $(srcdir)/Makefile.am $(am__configure_deps)
@for dep in $?; do \
case '$(am__configure_deps)' in \
......@@ -165,9 +212,9 @@ $(srcdir)/Makefile.in: $(srcdir)/Makefile.am $(am__configure_deps)
exit 1;; \
esac; \
done; \
echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu src/fits/Makefile'; \
echo ' cd $(top_srcdir) && $(AUTOMAKE) --foreign src/fits/Makefile'; \
cd $(top_srcdir) && \
$(AUTOMAKE) --gnu src/fits/Makefile
$(AUTOMAKE) --foreign src/fits/Makefile
.PRECIOUS: Makefile
Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status
@case '$?' in \
......@@ -213,27 +260,39 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/fitswrite.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
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@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
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ `$(CYGPATH_W) '$<'`
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Po
@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:
.c.lo:
@am__fastdepCC_TRUE@ $(LTCOMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/$*.Tpo $(DEPDIR)/$*.Plo
@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='$<' object='$@' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(LTCOMPILE) -c -o $@ $<
mostlyclean-libtool:
-rm -f *.lo
clean-libtool:
-rm -rf .libs _libs
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; }'`; \
$(AWK) '{ files[$$0] = 1; nonemtpy = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
mkid -fID $$unique
tags: TAGS
......@@ -245,8 +304,8 @@ TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \
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; }'`; \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { 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) \
......@@ -256,13 +315,12 @@ 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; }'`; \
$(AWK) '{ files[$$0] = 1; nonempty = 1; } \
END { if (nonempty) { for (i in files) print i; }; }'`; \
test -z "$(CTAGS_ARGS)$$tags$$unique" \
|| $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \
$$tags $$unique
......@@ -276,22 +334,21 @@ 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; \
@srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
topsrcdirstrip=`echo "$(top_srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \
list='$(DISTFILES)'; \
dist_files=`for file in $$list; do echo $$file; done | \
sed -e "s|^$$srcdirstrip/||;t" \
-e "s|^$$topsrcdirstrip/|$(top_builddir)/|;t"`; \
case $$dist_files in \
*/*) $(MKDIR_P) `echo "$$dist_files" | \
sed '/\//!d;s|^|$(distdir)/|;s,/[^/]*$$,,' | \
sort -u` ;; \
esac; \
for file in $$dist_files; do \
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 \
dir=`echo "/$$file" | sed -e 's,/[^/]*$$,,'`; \
if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \
cp -pR $(srcdir)/$$file $(distdir)$$dir || exit 1; \
fi; \
......@@ -332,7 +389,8 @@ maintainer-clean-generic:
@echo "it deletes files that may require special tools to rebuild."
clean: clean-am
clean-am: clean-generic clean-noinstLIBRARIES mostlyclean-am
clean-am: clean-generic clean-libtool clean-noinstLIBRARIES \
mostlyclean-am
distclean: distclean-am
-rm -rf ./$(DEPDIR)
......@@ -352,12 +410,20 @@ info-am:
install-data-am:
install-dvi: install-dvi-am
install-exec-am:
install-html: install-html-am
install-info: install-info-am
install-man:
install-pdf: install-pdf-am
install-ps: install-ps-am
installcheck-am:
maintainer-clean: maintainer-clean-am
......@@ -367,7 +433,8 @@ maintainer-clean-am: distclean-am maintainer-clean-generic
mostlyclean: mostlyclean-am
mostlyclean-am: mostlyclean-compile mostlyclean-generic
mostlyclean-am: mostlyclean-compile mostlyclean-generic \
mostlyclean-libtool
pdf: pdf-am
......@@ -377,18 +444,22 @@ ps: ps-am
ps-am:
uninstall-am: uninstall-info-am
uninstall-am:
.MAKE: install-am install-strip
.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
clean-libtool clean-noinstLIBRARIES ctags distclean \
distclean-compile distclean-generic distclean-libtool \
distclean-tags distdir dvi dvi-am html html-am info info-am \
install install-am install-data install-data-am install-dvi \
install-dvi-am install-exec install-exec-am install-html \
install-html-am install-info install-info-am install-man \
install-pdf install-pdf-am install-ps install-ps-am \
install-strip installcheck installcheck-am installdirs \
maintainer-clean maintainer-clean-generic mostlyclean \
mostlyclean-compile mostlyclean-generic mostlyclean-libtool \
pdf pdf-am ps ps-am tags uninstall uninstall-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.
......
src/fits/fitsbody.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: Handle memory allocation for FITS bodies.
*
* Last modify: 29/06/2006
* Last modify: 10/10/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -36,7 +36,6 @@ size_t body_maxram = BODY_DEFRAM,
int body_vmnumber;
double bufdata0[DATA_BUFSIZE/sizeof(double)];
char body_swapdirname[MAXCHARS] = BODY_DEFSWAPDIR;
/******* alloc_body ***********************************************************
......@@ -185,18 +184,19 @@ void free_body(tabstruct *tab)
/******* read_body ************************************************************
PROTO read_body(tabstruct *tab, PIXTYPE *ptr, long size)
PURPOSE Read values from the body of a FITS table.
PURPOSE Read floating point values from the body of a FITS table.
INPUT A pointer to the tab structure,
a pointer to the array in memory,
the number of elements to be read.
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 06/11/2003
VERSION 10/10/2007
***/
void read_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
{
catstruct *cat;
static double bufdata0[DATA_BUFSIZE/sizeof(double)];
unsigned char cuval, cublank;
char *bufdata,
cval, cblank;
......@@ -454,6 +454,174 @@ void read_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
}
/******* read_ibody ***********************************************************
PROTO read_ibody(tabstruct *tab, FLAGTYPE *ptr, long size)
PURPOSE Read integer values from the body of a FITS table.
INPUT A pointer to the tab structure,
a pointer to the array in memory,
the number of elements to be read.
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 11/10/2007
***/
void read_ibody(tabstruct *tab, FLAGTYPE *ptr, size_t size)
{
catstruct *cat;
static int bufdata0[DATA_BUFSIZE/sizeof(int)];
char *bufdata;
short val16;
int i, bowl, spoonful, npix, curval, dval;
/* a NULL cat structure indicates that no data can be read */
if (!(cat = tab->cat))
return;
switch(tab->compress_type)
{
/*-- Uncompressed image */
case COMPRESS_NONE:
bowl = DATA_BUFSIZE/tab->bytepix;
spoonful = size<bowl?size:bowl;
for(; size>0; size -= spoonful)
{
if (spoonful>size)
spoonful = size;
bufdata = (char *)bufdata0;
QFREAD(bufdata, spoonful*tab->bytepix, cat->file, cat->filename);
switch(tab->bitpix)
{
case BP_BYTE:
for (i=spoonful; i--;)
*(ptr++) = (FLAGTYPE)*((unsigned char *)bufdata++);
break;
case BP_SHORT:
if (bswapflag)
swapbytes(bufdata, 2, spoonful);
for (i=spoonful; i--; bufdata += sizeof(unsigned short))
*(ptr++) = (FLAGTYPE)*((unsigned short *)bufdata);
break;
case BP_LONG:
if (bswapflag)
swapbytes(bufdata, 4, spoonful);
for (i=spoonful; i--; bufdata += sizeof(unsigned long))
*(ptr++) = (FLAGTYPE)*((unsigned long *)bufdata);
break;
case BP_FLOAT:
case BP_DOUBLE:
error(EXIT_FAILURE,"*Error*: I was expecting integers in ",
cat->filename);
break;
default:
error(EXIT_FAILURE,"*FATAL ERROR*: unknown BITPIX type in ",
"readdata()");
break;
}
}
break;
/*-- Compressed image */
case COMPRESS_BASEBYTE:
if (!tab->compress_buf)
QMALLOC(tab->compress_buf, char, FBSIZE);
bufdata = tab->compress_bufptr;
curval = tab->compress_curval;
npix = tab->compress_npix;
while (size--)
{
if (!(npix--))
{
if (curval != tab->compress_checkval)
error(EXIT_FAILURE, "*Error*: invalid BASEBYTE checksum in ",
cat->filename);
bufdata = tab->compress_buf;
QFREAD(bufdata, FBSIZE, cat->file, cat->filename);
curval = 0;
if (bswapflag)
swapbytes(bufdata, 4, 1);
tab->compress_checkval = *((int *)bufdata);
bufdata += 4;
if (bswapflag)
swapbytes(bufdata, 2, 1);
npix = (int)(*((short *)bufdata))-1;
bufdata+=2;
}
if ((dval=(int)*(bufdata++))==-128)
{
if (bswapflag)
swapbytes(bufdata, 2, 1);
memcpy(&val16, bufdata, 2);
dval = (int)val16;
if (dval==-32768)
{
bufdata += 2;
if (bswapflag)
swapbytes(bufdata, 4, 1);
memcpy(&dval,bufdata,4);
bufdata += 4;
}
else
bufdata += 2;
}
*(ptr++) = (FLAGTYPE)dval;
curval += dval;
}
tab->compress_curval = curval;
tab->compress_bufptr = bufdata;
tab->compress_npix = npix;
break;
case COMPRESS_PREVPIX:
if (!tab->compress_buf)
QMALLOC(tab->compress_buf, char, FBSIZE);
bufdata = tab->compress_bufptr;
curval = tab->compress_curval;
npix = tab->compress_npix;
while (size--)
{
if (!(npix--))
{
if (curval != tab->compress_checkval)
error(EXIT_FAILURE, "*Error*: invalid PREV_PIX checksum in ",
cat->filename);
bufdata = tab->compress_buf;
QFREAD(bufdata, FBSIZE, cat->file, cat->filename);
if (bswapflag)
swapbytes(bufdata, 2, 3);
curval = (int)*(short *)bufdata;
npix = (int)*(short *)(bufdata+=2)-1;
tab->compress_checkval = (int)(*(short *)(bufdata+=2));
bufdata+=4;
}
if ((dval=(int)*(bufdata++))==-128)
{
if (bswapflag)
swapbytes(bufdata, 2, 1);
memcpy(&val16, bufdata, 2);
curval = (int)val16;
bufdata += 2;
}
else
curval += dval;
*(ptr++) = (FLAGTYPE)curval;
}
tab->compress_curval = curval;
tab->compress_bufptr = bufdata;
tab->compress_npix = npix;
break;
default:
error(EXIT_FAILURE,"*Internal Error*: unknown compression mode in ",
"readdata()");
}
return;
}
/******* write_body ***********************************************************
PROTO write_body(tabstruct *tab, PIXTYPE *ptr, long size)
PURPOSE Write values to a FITS body.
......@@ -463,10 +631,11 @@ INPUT A pointer to the tab structure,
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 29/06/2006
VERSION 11/10/2007
***/
void write_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
{
static double bufdata0[DATA_BUFSIZE/sizeof(double)];
catstruct *cat;
char *cbufdata0;
size_t i, bowl, spoonful;
......
src/fits/fitscat.h
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: Simplified versin of the LDACTools: main include file
*
* Last modify: 10/07/2006
* Last modify: 10/10/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -22,7 +22,7 @@
#endif
#define MAXCHARS 256 /* max. number of characters */
#define WARNING_NMAX 100 /* max. number of recorded warnings */
#define WARNING_NMAX 1000 /* max. number of recorded warnings */
/*---------------------------- return messages ------------------------------*/
......@@ -75,6 +75,7 @@ typedef enum {SHOW_ASCII, SHOW_SKYCAT}
output_type; /* Type of output */
typedef float PIXTYPE; /* Pixel type */
typedef unsigned int FLAGTYPE; /* Flag type */
#ifdef HAVE_UNSIGNED_LONG_LONG
typedef unsigned long long KINGSIZE_T; /* for large sizes */
......@@ -212,6 +213,7 @@ extern void add_cleanupfilename(char *filename),
int nkeys, unsigned char *mask),
read_basic(tabstruct *tab),
read_body(tabstruct *tab, PIXTYPE *ptr, size_t size),
read_ibody(tabstruct *tab, FLAGTYPE *ptr, size_t size),
readbasic_head(tabstruct *tab),
remove_cleanupfilename(char *filename),
save_cat(catstruct *cat, char *filename),
......@@ -282,6 +284,7 @@ extern int about_cat(catstruct *cat, FILE *stream),
long pos),
remove_key(tabstruct *tab, char *keyname),
remove_keys(tabstruct *tab),
removekeywordfrom_head(tabstruct *tab, char *keyword),
remove_tab(catstruct *cat, char *tabname, int seg),
remove_tabs(catstruct *cat),
save_head(catstruct *cat, tabstruct *tab),
......
src/fits/fitscat_defs.h
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: Simplified version of the LDACTools: internal defs
*
* Last modify: 16/08/2004
* Last modify: 26/09/2006
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -44,6 +44,7 @@ extern int bswapflag; /* != 0 if bytes are swapped/IEEE */
#define BODY_DEFSWAPDIR "/tmp" /* OK at least for Unix systems */
#define BIG 1e+30 /* a huge number */
#define TINY (1.0/BIG) /* a tiny number */
#ifndef PI
#define PI 3.14159265359 /* never met before? */
#endif
......
src/fits/fitscleanup.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: Signal-catching routines to clean-up temporary files
*
* Last modify: 10/01/2003
* Last modify: 16/07/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -96,7 +96,7 @@ INPUT pointer to filename char string.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 21/08/2000
VERSION 16/07/2007
***/
void remove_cleanupfilename(char *filename)
{
......@@ -113,22 +113,17 @@ void remove_cleanupfilename(char *filename)
{
/* Match found: update the list and free memory is necessary*/
filename3 = filename2 - 1;
free(*filename3);
for (j=i; j--;)
*(filename3++) = *(filename2++);
if (!((--cleanup_nfiles)%CLEANUP_NFILES))
{
if (cleanup_nfiles)
{
filename2 = cleanup_filename + cleanup_nfiles;
for (i=CLEANUP_NFILES; i--;)
free(*(filename2++));
QREALLOC(cleanup_filename, char *, cleanup_nfiles);
}
else
{
free(*cleanup_filename);
free(cleanup_filename);
}
}
break;
}
......
src/fits/fitshead.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: general functions for handling FITS file headers.
*
* Last modify: 25/09/2004
* Last modify: 20/06/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -315,7 +315,7 @@ INPUT Table structure.
OUTPUT RETURN_OK if tab is a binary table, or RETURN_ERROR otherwise.
NOTES The headbuf pointer in the tabstruct might be reallocated.
AUTHOR E. Bertin (IAP & Leiden observatory)
VERSION 25/09/2004
VERSION 11/06/2007
***/
int update_head(tabstruct *tab)
......@@ -353,13 +353,13 @@ int update_head(tabstruct *tab)
}
/*First, remove all existing TTYPE, TFORM, etc...*/
fitsremove(tab->headbuf, "TTYPE???");
fitsremove(tab->headbuf, "TFORM???");
fitsremove(tab->headbuf, "TUNIT???");
fitsremove(tab->headbuf, "TZERO???");
fitsremove(tab->headbuf, "TSCAL???");
fitsremove(tab->headbuf, "TDIM???");
fitsremove(tab->headbuf, "TDISP???");
removekeywordfrom_head(tab, "TTYPE???");
removekeywordfrom_head(tab, "TFORM???");
removekeywordfrom_head(tab, "TUNIT???");
removekeywordfrom_head(tab, "TZERO???");
removekeywordfrom_head(tab, "TSCAL???");
removekeywordfrom_head(tab, "TDIM???");
removekeywordfrom_head(tab, "TDISP???");
/*Change NAXIS1 in order to take into account changes in width*/
......@@ -465,8 +465,8 @@ PURPOSE Update a FITS header to make it "primary" (not extension)
INPUT Table structure.
OUTPUT RETURN_OK if tab header was already primary, or RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP & Leiden observatory)
VERSION 08/05/2002
AUTHOR E. Bertin (IAP & Leiden observatory) C. Marmo (IAP)
VERSION 11/06/2007
***/
int prim_head(tabstruct *tab)
......@@ -477,8 +477,13 @@ int prim_head(tabstruct *tab)
{
strncpy(tab->headbuf, "SIMPLE = T "
"/ This is a FITS file ", 80);
/* fitsverify 4.13 (CFITSIO V3.002) return an error
if PCOUNT and GCOUNT are in a primary header (23/05/2007)*/
removekeywordfrom_head(tab, "PCOUNT");
removekeywordfrom_head(tab, "GCOUNT");
return RETURN_ERROR;
}
return RETURN_OK;
}
......@@ -490,8 +495,8 @@ INPUT Table structure.
OUTPUT RETURN_OK if tab header was already extension, or RETURN_ERROR
otherwise.
NOTES -.
AUTHOR E. Bertin (IAP & Leiden observatory)
VERSION 08/05/2002
AUTHOR E. Bertin (IAP & Leiden observatory) C. Marmo (IAP)
VERSION 20/06/2007
***/
int ext_head(tabstruct *tab)
......@@ -502,6 +507,15 @@ int ext_head(tabstruct *tab)
{
strncpy(tab->headbuf, "XTENSION= 'IMAGE ' "
"/ Image extension ", 80);
/* fitsverify 4.13 (CFITSIO V3.002) return an error
if EXTEND are in an extension header (20/06/2007)*/
removekeywordfrom_head(tab, "EXTEND");
/* fitsverify 4.13 (CFITSIO V3.002) return an error
if PCOUNT and GCOUNT are not in the extension header (23/05/2007) */
addkeywordto_head(tab, "PCOUNT ", "required keyword; must = 0");
addkeywordto_head(tab, "GCOUNT ", "required keyword; must = 1");
fitswrite(tab->headbuf,"PCOUNT ", &tab->pcount, H_INT, T_LONG);
fitswrite(tab->headbuf,"GCOUNT ", &tab->gcount, H_INT, T_LONG);
return RETURN_ERROR;
}
......@@ -564,6 +578,36 @@ int addkeywordto_head(tabstruct *tab, char *keyword, char *comment)
}
/****** removekeywordfrom_head ************************************************
PROTO int removekeywordfrom_head(tabstruct *tab, char *keyword)
PURPOSE Remove a keyword from a table header.
INPUT Table structure,
String containing the keyword.
OUTPUT RETURN_OK if the keyword was found, RETURN_ERROR otherwise..
NOTES The headbuf pointer in the tabstruct might be reallocated.
'?' wildcard allowed; Don't remove the ``END'' keyword with this!!!
AUTHOR E. Bertin (IAP)
VERSION 11/06/2007
***/
int removekeywordfrom_head(tabstruct *tab, char *keyword)
{
int nb;
if (fitsremove(tab->headbuf, keyword) == RETURN_OK)
{
if ((nb=fitsfind(tab->headbuf, "END ")/(FBSIZE/80)+1) < tab->headnblock)
{
tab->headnblock = nb;
QREALLOC(tab->headbuf, char, tab->headnblock*FBSIZE);
}
return RETURN_OK;
}
else
return RETURN_ERROR;
}
/****** tformof ***************************************************************
PROTO int tformof(char *str, t_type ttype, int n)
PURPOSE Return the ``TFORM'' string corresponding to a t_type
......
src/fits/fitskey.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: Functions related to the management of keys.
*
* Last modify: 15/08/2003
* Last modify: 04/06/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -575,7 +575,7 @@ NOTES This is approximately the same code as for read_keys.
A NULL keynames pointer means read ALL keys belonging to the table.
A NULL mask pointer means NO selection for reading.
AUTHOR E. Bertin (IAP & Leiden observatory)
VERSION 15/08/2003
VERSION 04/06/2007
***/
void show_keys(tabstruct *tab, char **keynames, keystruct **keys, int nkeys,
BYTE *mask, FILE *stream,
......@@ -691,10 +691,10 @@ void show_keys(tabstruct *tab, char **keynames, keystruct **keys, int nkeys,
if (banflag)
{
if (*key->unit)
fprintf(stream, "# %3d %-19.19s %-47.47s [%s]\n",
fprintf(stream, "# %3d %-15.15s %-47.47s [%s]\n",
n, key->name,key->comment, key->unit);
else
fprintf(stream, "# %3d %-19.19s %.47s\n",
fprintf(stream, "# %3d %-15.15s %.47s\n",
n, key->name,key->comment);
n += key->nbytes/t_size[key->ttype];
}
......@@ -730,10 +730,10 @@ void show_keys(tabstruct *tab, char **keynames, keystruct **keys, int nkeys,
if (banflag)
{
if (*key->unit)
fprintf(stream, "# %3d %-19.19s %-47.47s [%s]\n",
fprintf(stream, "# %3d %-15.15s %-47.47s [%s]\n",
n, key->name,key->comment, key->unit);
else
fprintf(stream, "# %3d %-19.19s %.47s\n",
fprintf(stream, "# %3d %-15.15s %.47s\n",
n, key->name,key->comment);
n += key->nbytes/t_size[key->ttype];
}
......@@ -754,10 +754,10 @@ void show_keys(tabstruct *tab, char **keynames, keystruct **keys, int nkeys,
switch (o_type) {
case SHOW_ASCII:
if (*key->unit)
fprintf(stream, "# %-19.19s %-47.47s [%s]\n",
fprintf(stream, "# %-15.15s %-47.47s [%s]\n",
key->name,key->comment, key->unit);
else
fprintf(stream, "# %-19.19s %.47s\n",
fprintf(stream, "# %-15.15s %.47s\n",
key->name,key->comment);
break;
case SHOW_SKYCAT:
......
src/fits/fitsmisc.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: miscellaneous functions.
*
* Last modify: 14/07/2006
* Last modify: 18/07/2008
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -113,8 +113,9 @@ Swap bytes for doubles, longs and shorts (for DEC machines or PC for inst.).
*/
void swapbytes(void *ptr, int nb, int n)
{
char *cp;
int j;
char *cp,
c;
int j;
cp = (char *)ptr;
......@@ -122,8 +123,12 @@ void swapbytes(void *ptr, int nb, int n)
{
for (j=n; j--; cp+=4)
{
cp[0] ^= (cp[3]^=(cp[0]^=cp[3]));
cp[1] ^= (cp[2]^=(cp[1]^=cp[2]));
c = cp[3];
cp[3] = cp[0];
cp[0] = c;
c = cp[2];
cp[2] = cp[1];
cp[1] = c;
}
return;
}
......@@ -131,7 +136,11 @@ void swapbytes(void *ptr, int nb, int n)
if (nb&2)
{
for (j=n; j--; cp+=2)
cp[0] ^= (cp[1]^=(cp[0]^=cp[1]));
{
c = cp[1];
cp[1] = cp[0];
cp[0] = c;
}
return;
}
......@@ -142,10 +151,18 @@ void swapbytes(void *ptr, int nb, int n)
{
for (j=n; j--; cp+=8)
{
cp[0] ^= (cp[7]^=(cp[0]^=cp[7]));
cp[1] ^= (cp[6]^=(cp[1]^=cp[6]));
cp[2] ^= (cp[5]^=(cp[2]^=cp[5]));
cp[3] ^= (cp[4]^=(cp[3]^=cp[4]));
c = cp[7];
cp[7] = cp[0];
cp[0] = c;
c = cp[6];
cp[6] = cp[1];
cp[1] = c;
c = cp[5];
cp[5] = cp[2];
cp[2] = c;
c = cp[4];
cp[4] = cp[3];
cp[3] = c;
}
return;
}
......
src/fits/fitsutil.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: functions for handling FITS keywords.
*
* Last modify: 21/09/2006
* Last modify: 12/06/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -38,8 +38,8 @@ NOTES For all keywords except commentary ones (like COMMENT, HISTORY or
blank), it is checked that they do not exist already.
Enough memory should be provided for the FITS header to contain one
more line of 80 char.
AUTHOR E. Bertin (IAP & Leiden observatory)
VERSION 13/06/2004
AUTHOR E. Bertin (IAP & Leiden observatory) C. Marmo (IAP)
VERSION 13/06/2007
***/
int fitsadd(char *fitsbuf, char *keyword, char *comment)
......@@ -76,6 +76,26 @@ int fitsadd(char *fitsbuf, char *keyword, char *comment)
else
return RETURN_ERROR;
}
/*-- Special case of PCOUNT/GCOUNT parameters */
if (!strncmp(keyword, "PCOUNT", 6))
{
headpos=fitsfind(fitsbuf, "NAXIS ");
sscanf(fitsbuf+80*headpos, "NAXIS = %d", &n);
if (headpos>0)
headpos+=(n+1);
else
return RETURN_ERROR;
}
if (!strncmp(keyword, "GCOUNT", 6))
{
headpos=fitsfind(fitsbuf, "NAXIS ");
sscanf(fitsbuf+80*headpos, "NAXIS = %d", &n);
if (headpos>0)
headpos+=(n+2);
else
return RETURN_ERROR;
}
key_ptr = fitsbuf+80*headpos;
memmove(key_ptr+80, key_ptr, 80*(headpos2-headpos+1));
......@@ -88,7 +108,6 @@ int fitsadd(char *fitsbuf, char *keyword, char *comment)
else
sprintf(str, "%-8.8s= %-47.47s",
keyword, " ");
memcpy(key_ptr, str, 80);
}
......@@ -162,7 +181,7 @@ OUTPUT RETURN_OK if something was found, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP),
E.R. Deul - Handling of NaN
VERSION 19/09/2006
VERSION 04/06/2007
***/
int fitspick(char *fitsline, char *keyword, void *ptr, h_type *htype,
t_type *ttype, char *comment)
......@@ -233,13 +252,7 @@ int fitspick(char *fitsline, char *keyword, void *ptr, h_type *htype,
{
for (i=j; i<80 && fitsline[i]!=(char)'/' && fitsline[i]!=(char)'.'; i++);
/*-- Handle floats*/
if (i==80)
{
*((int *)ptr) = 0;
*htype = H_INT;
*ttype = T_LONG;
}
else if (fitsline[i]==(char)'.')
if (fitsline[i]==(char)'.')
{
fixexponent(fitsline);
*((double *)ptr) = atof(fitsline+j);
......
src/fitswcs.c 0 → 100644
View file @ ece261cd
/*
fitswcs.c
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: LDACTools+
*
* Author: E.BERTIN (IAP)
*
* Contents: Read and write WCS header info.
*
* Last modify: 26/04/2008
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#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 "fits/fitscat_defs.h"
#include "fits/fitscat.h"
#include "fitswcs.h"
#include "wcscelsys.h"
#include "wcs/wcs.h"
#include "wcs/lin.h"
#include "wcs/tnx.h"
#include "wcs/poly.h"
/******* copy_wcs ************************************************************
PROTO wcsstruct *copy_wcs(wcsstruct *wcsin)
PURPOSE Copy a WCS (World Coordinate System) structure.
INPUT WCS structure to be copied.
OUTPUT pointer to a copy of the input structure.
NOTES Actually, only FITS parameters are copied. Lower-level structures
such as those created by the WCS or TNX libraries are generated.
AUTHOR E. Bertin (IAP)
VERSION 31/08/2002
***/
wcsstruct *copy_wcs(wcsstruct *wcsin)
{
wcsstruct *wcs;
/* Copy the basic stuff */
QMEMCPY(wcsin, wcs, wcsstruct, 1);
/* The PROJP WCS parameters */
QMEMCPY(wcsin->projp, wcs->projp, double, wcs->naxis*100);
/* Set other structure pointers to NULL (they'll have to be reallocated) */
wcs->wcsprm = NULL;
wcs->lin = NULL;
wcs->cel = NULL;
wcs->prj = NULL;
wcs->tnx_lngcor = copy_tnxaxis(wcsin->tnx_lngcor);
wcs->tnx_latcor = copy_tnxaxis(wcsin->tnx_latcor);
wcs->inv_x = wcs->inv_y = NULL;
QCALLOC(wcs->wcsprm, struct wcsprm, 1);
/* Test if the WCS is recognized and a celestial pair is found */
wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);
/* Initialize other WCS structures */
init_wcs(wcs);
/* Invert projection corrections */
invert_wcs(wcs);
/* Find the range of coordinates */
range_wcs(wcs);
return wcs;
}
/******* create_wcs ***********************************************************
PROTO wcsstruct *create_wcs(char **ctype, double *crval, double *crpix,
double *cdelt, int *naxisn, int naxis)
PURPOSE Generate a simple WCS (World Coordinate System) structure.
INPUT Pointer to an array of char strings with WCS projection on each axis,
pointer to an array of center coordinates (double),
pointer to an array of device coordinates (double),
pointer to an array of pixel scales (double),
pointer to an array of image dimensions (int),
number of dimensions.
OUTPUT pointer to a WCS structure.
NOTES If a pointer is set to null, the corresponding variables are set to
default values.
AUTHOR E. Bertin (IAP)
VERSION 09/08/2006
***/
wcsstruct *create_wcs(char **ctype, double *crval, double *crpix,
double *cdelt, int *naxisn, int naxis)
{
wcsstruct *wcs;
int l;
QCALLOC(wcs, wcsstruct, 1);
wcs->naxis = naxis;
QCALLOC(wcs->projp, double, naxis*100);
wcs->nprojp = 0;
wcs->longpole = wcs->latpole = 999.0;
for (l=0; l<naxis; l++)
{
wcs->naxisn[l] = naxisn? naxisn[l] : 360.0;
/*-- The default WCS projection system is an all-sky Aitoff projection */
if (ctype)
strncpy(wcs->ctype[l], ctype[l], 8);
else if (l==0)
strncpy(wcs->ctype[l], "RA---AIT", 8);
else if (l==1)
strncpy(wcs->ctype[l], "DEC--AIT", 8);
wcs->crval[l] = crval? crval[l]: 0.0;
wcs->crpix[l] = crpix? crpix[l]: 0.0;
wcs->cdelt[l] = 1.0;
wcs->cd[l*(naxis+1)] = cdelt? cdelt[l] : 1.0;
}
wcs->epoch = wcs->equinox = 2000.0;
QCALLOC(wcs->wcsprm, struct wcsprm, 1);
/* Test if the WCS is recognized and a celestial pair is found */
wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);
/* Initialize other WCS structures */
init_wcs(wcs);
/* Invert projection corrections */
invert_wcs(wcs);
/* Find the range of coordinates */
range_wcs(wcs);
return wcs;
}
/******* init_wcs ************************************************************
PROTO void init_wcs(wcsstruct *wcs)
PURPOSE Initialize astrometry and WCS (World Coordinate System) structures.
INPUT WCS structure.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 17/05/2007
***/
void init_wcs(wcsstruct *wcs)
{
int l,n,lng,lat,naxis;
naxis = wcs->naxis;
if (wcs->lin)
{
free(wcs->lin->cdelt);
free(wcs->lin->crpix);
free(wcs->lin->pc);
free(wcs->lin->piximg);
free(wcs->lin->imgpix);
free(wcs->lin);
}
QCALLOC(wcs->lin, struct linprm, 1);
QCALLOC(wcs->lin->cdelt, double, naxis);
QCALLOC(wcs->lin->crpix, double, naxis);
QCALLOC(wcs->lin->pc, double, naxis*naxis);
if (wcs->cel)
free(wcs->cel);
QCALLOC(wcs->cel, struct celprm, 1);
if (wcs->prj)
free(wcs->prj);
QCALLOC(wcs->prj, struct prjprm, 1);
if (wcs->inv_x)
{
poly_end(wcs->inv_x);
wcs->inv_x = NULL;
}
if (wcs->inv_y)
{
poly_end(wcs->inv_y);
wcs->inv_y = NULL;
}
/* Set WCS flags to 0: structures will be reinitialized by the WCS library */
wcs->lin->flag = wcs->cel->flag = wcs->prj->flag = 0;
wcs->lin->naxis = naxis;
/* wcsprm structure */
lng = wcs->lng = wcs->wcsprm->lng;
lat = wcs->lat = wcs->wcsprm->lat;
/* linprm structure */
for (l=0; l<naxis; l++)
{
wcs->lin->crpix[l] = wcs->crpix[l];
wcs->lin->cdelt[l] = 1.0;
}
for (l=0; l<naxis*naxis; l++)
wcs->lin->pc[l] = wcs->cd[l];
/* celprm structure */
if (lng>=0)
{
wcs->cel->ref[0] = wcs->crval[lng];
wcs->cel->ref[1] = wcs->crval[lat];
}
else
{
wcs->cel->ref[0] = wcs->crval[0];
wcs->cel->ref[1] = wcs->crval[1];
}
wcs->cel->ref[2] = wcs->longpole;
wcs->cel->ref[3] = wcs->latpole;
/* prjprm structure */
wcs->prj->r0 = wcs->r0;
wcs->prj->tnx_lngcor = wcs->tnx_lngcor;
wcs->prj->tnx_latcor = wcs->tnx_latcor;
if (lng>=0)
{
n = 0;
for (l=100; l--;)
{
wcs->prj->p[l] = wcs->projp[l+lat*100]; /* lat comes first for ... */
wcs->prj->p[l+100] = wcs->projp[l+lng*100];/* ... compatibility reasons */
if (!n && (wcs->prj->p[l] || wcs->prj->p[l+100]))
n = l+1;
}
wcs->nprojp = n;
}
/* Check-out chirality */
wcs->chirality = wcs_chirality(wcs);
/* Initialize Equatorial <=> Celestial coordinate system transforms */
init_wcscelsys(wcs);
return;
}
/******* init_wcscelsys *******************************************************
PROTO void init_wcscelsys(wcsstruct *wcs)
PURPOSE Initialize Equatorial <=> Celestial coordinate system transforms.
INPUT WCS structure.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 18/07/2006
***/
void init_wcscelsys(wcsstruct *wcs)
{
double *mat,
a0,d0,ap,dp,ap2,y;
int s,lng,lat;
lng = wcs->wcsprm->lng;
lat = wcs->wcsprm->lat;
/* Is it a celestial system? If not, exit! */
if (lng==lat)
{
wcs->celsysconvflag = 0;
return;
}
/* Find the celestial system */
for (s=0; *celsysname[s][0] && strncmp(wcs->ctype[lng], celsysname[s][0], 4);
s++);
/* Is it a known, non-equatorial system? If not, exit! */
if (!s || !*celsysname[s][0])
{
wcs->celsysconvflag = 0;
return;
}
wcs->celsys = (celsysenum)s;
/* Some shortcuts */
a0 = celsysorig[s][0]*DEG;
d0 = celsysorig[s][1]*DEG;
ap = celsyspole[s][0]*DEG;
dp = celsyspole[s][1]*DEG;
/* First compute in the output referential the longitude of the south pole */
y = sin(ap - a0);
/*
x = cos(d0)*(cos(d0)*sin(dp)*cos(ap-a0)-sin(d0)*cos(dp));
ap2 = atan2(y,x);
*/
ap2 = asin(cos(d0)*y) ;
/* Equatorial <=> Celestial System transformation parameters */
mat = wcs->celsysmat;
mat[0] = ap;
mat[1] = ap2;
mat[2] = cos(dp);
mat[3] = sin(dp);
wcs->celsysconvflag = 1;
return;
}
/******* read_wcs *************************************************************
PROTO wcsstruct *read_wcs(tabstruct *tab)
PURPOSE Read WCS (World Coordinate System) info in the FITS header.
INPUT tab structure.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 02/01/2008
***/
wcsstruct *read_wcs(tabstruct *tab)
{
#define FITSREADF(buf, k, val, def) \
{if (fitsread(buf,k, &val, H_FLOAT,T_DOUBLE) != RETURN_OK) \
val = def; \
}
#define FITSREADI(buf, k, val, def) \
{if (fitsread(buf,k, &val, H_INT,T_LONG) != RETURN_OK) \
val = def; \
}
#define FITSREADS(buf, k, str, def) \
{if (fitsread(buf,k,str, H_STRING,T_STRING) != RETURN_OK) \
strcpy(str, (def)); \
}
char str[MAXCHARS];
char wstr1[TNX_MAXCHARS], wstr2[TNX_MAXCHARS];
wcsstruct *wcs;
double drota;
int j, l, naxis;
char name[16],
*buf, *filename, *ptr;
buf = tab->headbuf;
filename = (tab->cat? tab->cat->filename : strcpy(name, "internal header"));
FITSREADS(buf, "OBJECT ", str, "Unnamed");
QCALLOC(wcs, wcsstruct, 1);
if (tab->naxis > NAXIS)
{
warning("Maximum number of dimensions supported by this version of the ",
"software exceeded\n");
tab->naxis = 2;
}
wcs->naxis = naxis = tab->naxis;
QCALLOC(wcs->projp, double, naxis*100);
for (l=0; l<naxis; l++)
{
wcs->naxisn[l] = tab->naxisn[l];
sprintf(str, "CTYPE%-3d", l+1);
FITSREADS(buf, str, str, "");
strncpy(wcs->ctype[l], str, 8);
sprintf(str, "CUNIT%-3d", l+1);
FITSREADS(buf, str, str, "deg");
strncpy(wcs->cunit[l], str, 32);
sprintf(str, "CRVAL%-3d", l+1);
FITSREADF(buf, str, wcs->crval[l], 0.0);
sprintf(str, "CRPIX%-3d", l+1);
FITSREADF(buf, str, wcs->crpix[l], 1.0);
sprintf(str, "CDELT%-3d", l+1);
FITSREADF(buf, str, wcs->cdelt[l], 1.0);
sprintf(str, "CRDER%-3d", l+1);
FITSREADF(buf, str, wcs->crder[l], 0.0);
sprintf(str, "CSYER%-3d", l+1);
FITSREADF(buf, str, wcs->csyer[l], 0.0);
if (fabs(wcs->cdelt[l]) < 1e-30)
error(EXIT_FAILURE, "*Error*: CDELT parameters out of range in ",
filename);
}
if (fitsfind(buf, "CD?_????")!=RETURN_ERROR)
{
/*-- If CD keywords exist, use them for the linear mapping terms... */
for (l=0; l<naxis; l++)
for (j=0; j<naxis; j++)
{
sprintf(str, "CD%d_%d", l+1, j+1);
FITSREADF(buf, str, wcs->cd[l*naxis+j], l==j?1.0:0.0)
}
}
else if (fitsfind(buf, "PC00?00?")!=RETURN_ERROR)
/*-- ...If PC keywords exist, use them for the linear mapping terms... */
for (l=0; l<naxis; l++)
for (j=0; j<naxis; j++)
{
sprintf(str, "PC%03d%03d", l+1, j+1);
FITSREADF(buf, str, wcs->cd[l*naxis+j], l==j?1.0:0.0)
wcs->cd[l*naxis+j] *= wcs->cdelt[l];
}
else
{
/*-- ...otherwise take the obsolete CROTA2 parameter */
FITSREADF(buf, "CROTA2 ", drota, 0.0)
wcs->cd[3] = wcs->cd[0] = cos(drota*DEG);
wcs->cd[1] = -(wcs->cd[2] = sin(drota*DEG));
wcs->cd[0] *= wcs->cdelt[0];
wcs->cd[2] *= wcs->cdelt[0];
wcs->cd[1] *= wcs->cdelt[1];
wcs->cd[3] *= wcs->cdelt[1];
}
QCALLOC(wcs->wcsprm, struct wcsprm, 1);
/* Test if the WCS is recognized and a celestial pair is found */
if (!wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm)
&& wcs->wcsprm->flag<999)
{
char *pstr;
double date;
int biss, dpar[3];
/*-- Coordinate reference frame */
/*-- Search for an observation date expressed in Julian days */
FITSREADF(buf, "MJD-OBS ", date, -1.0);
/*-- Precession date (defined from Ephemerides du Bureau des Longitudes) */
/*-- in Julian years from 2000.0 */
if (date>0.0)
wcs->obsdate = 2000.0 - (MJD2000 - date)/365.25;
else
{
/*---- Search for an observation date expressed in "civilian" format */
FITSREADS(buf, "DATE-OBS ", str, "");
if (*str)
{
/*------ Decode DATE-OBS format: DD/MM/YY or YYYY-MM-DD */
for (l=0; l<3 && (pstr = strtok_r(l?NULL:str,"/- ", &ptr)); l++)
dpar[l] = atoi(pstr);
if (l<3 || !dpar[0] || !dpar[1] || !dpar[2])
{
/*-------- If DATE-OBS value corrupted or incomplete, assume 2000-1-1 */
warning("Invalid DATE-OBS value in header: ", str);
dpar[0] = 2000; dpar[1] = 1; dpar[2] = 1;
}
else if (strchr(str, '/') && dpar[0]<32 && dpar[2]<100)
{
j = dpar[0];
dpar[0] = dpar[2]+1900;
dpar[2] = j;
}
biss = (dpar[0]%4)?0:1;
/*------ Convert date to MJD */
date = -678956 + (365*dpar[0]+dpar[0]/4) - biss
+ ((dpar[1]>2?((int)((dpar[1]+1)*30.6)-63+biss)
:((dpar[1]-1)*(63+biss))/2) + dpar[2]);
wcs->obsdate = 2000.0 - (MJD2000 - date)/365.25;
}
else
/*------ Well if really no date is found */
wcs->obsdate = 0.0;
}
FITSREADF(buf, "EPOCH", wcs->epoch, 2000.0);
FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->epoch);
FITSREADS(buf, "RADECSYS", str,
wcs->equinox >= 2000.0? "ICRS" : (wcs->equinox<1984.0? "FK4" : "FK5"));
if (!strcmp(str, "ICRS"))
wcs->radecsys = RDSYS_ICRS;
else if (!strcmp(str, "FK5"))
wcs->radecsys = RDSYS_FK5;
else if (!strcmp(str, "FK4"))
{
if (wcs->equinox == 2000.0)
{
FITSREADF(buf, "EPOCH ", wcs->equinox, 1950.0);
FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->equinox);
}
wcs->radecsys = RDSYS_FK4;
warning("FK4 precession formulae not yet implemented:\n",
" Astrometry may be slightly inaccurate");
}
else if (!strcmp(str, "FK4-NO-E"))
{
if (wcs->equinox == 2000.0)
{
FITSREADF(buf, "EPOCH", wcs->equinox, 1950.0);
FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->equinox);
}
wcs->radecsys = RDSYS_FK4_NO_E;
warning("FK4 precession formulae not yet implemented:\n",
" Astrometry may be slightly inaccurate");
}
else if (!strcmp(str, "GAPPT"))
{
wcs->radecsys = RDSYS_GAPPT;
warning("GAPPT reference frame not yet implemented:\n",
" Astrometry may be slightly inaccurate");
}
else
{
warning("Using ICRS instead of unknown astrometric reference frame: ",
str);
wcs->radecsys = RDSYS_ICRS;
}
/*-- Projection parameters */
if (!strcmp(wcs->wcsprm->pcode, "TNX"))
{
/*---- IRAF's TNX projection: decode these #$!?@#!! WAT parameters */
if (fitsfind(buf, "WAT?????") != RETURN_ERROR)
{
/*------ First we need to concatenate strings */
pstr = wstr1;
sprintf(str, "WAT1_001");
for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
{
sprintf(str, "WAT1_%03d", j);
pstr += strlen(pstr);
}
pstr = wstr2;
sprintf(str, "WAT2_001");
for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
{
sprintf(str, "WAT2_%03d", j);
pstr += strlen(pstr);
}
/*------ LONGPOLE defaulted to 180 deg if not found */
if ((pstr = strstr(wstr1, "longpole"))
|| (pstr = strstr(wstr2, "longpole")))
pstr = strpbrk(pstr, "1234567890-+.");
wcs->longpole = pstr? atof(pstr) : 999.0;
wcs->latpole = 999.0;
/*------ RO defaulted to 180/PI if not found */
if ((pstr = strstr(wstr1, "ro"))
|| (pstr = strstr(wstr2, "ro")))
pstr = strpbrk(pstr, "1234567890-+.");
wcs->r0 = pstr? atof(pstr) : 0.0;
/*------ Read the remaining TNX parameters */
if ((pstr = strstr(wstr1, "lngcor"))
|| (pstr = strstr(wstr2, "lngcor")))
wcs->tnx_lngcor = read_tnxaxis(pstr);
if (!wcs->tnx_lngcor)
error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
filename);
if ((pstr = strstr(wstr1, "latcor"))
|| (pstr = strstr(wstr2, "latcor")))
wcs->tnx_latcor = read_tnxaxis(pstr);
if (!wcs->tnx_latcor)
error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
filename);
}
}
else
{
FITSREADF(buf, "LONGPOLE", wcs->longpole, 999.0);
FITSREADF(buf, "LATPOLE ", wcs->latpole, 999.0);
/*---- Old convention */
if (fitsfind(buf, "PROJP???") != RETURN_ERROR)
for (j=0; j<10; j++)
{
sprintf(str, "PROJP%-3d", j);
FITSREADF(buf, str, wcs->projp[j], 0.0);
}
/*---- New convention */
if (fitsfind(buf, "PV?_????") != RETURN_ERROR)
for (l=0; l<naxis; l++)
for (j=0; j<100; j++)
{
sprintf(str, "PV%d_%d", l+1, j);
FITSREADF(buf, str, wcs->projp[j+l*100], 0.0);
}
}
}
/* Initialize other WCS structures */
init_wcs(wcs);
/* Find the range of coordinates */
range_wcs(wcs);
/* Invert projection corrections */
invert_wcs(wcs);
#undef FITSREADF
#undef FITSREADI
#undef FITSREADS
return wcs;
}
/******* write_wcs ***********************************************************
PROTO void write_wcs(tabstruct *tab, wcsstruct *wcs)
PURPOSE Write WCS (World Coordinate System) info in the FITS header.
INPUT tab structure,
WCS structure.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 17/07/2006
***/
void write_wcs(tabstruct *tab, wcsstruct *wcs)
{
char str[MAXCHARS];
int j, l, naxis;
naxis = wcs->naxis;
addkeywordto_head(tab, "BITPIX ", "Bits per pixel");
fitswrite(tab->headbuf, "BITPIX ", &tab->bitpix, H_INT, T_LONG);
addkeywordto_head(tab, "NAXIS ", "Number of axes");
fitswrite(tab->headbuf, "NAXIS ", &wcs->naxis, H_INT, T_LONG);
for (l=0; l<naxis; l++)
{
sprintf(str, "NAXIS%-3d", l+1);
addkeywordto_head(tab, str, "Number of pixels along this axis");
fitswrite(tab->headbuf, str, &wcs->naxisn[l], H_INT, T_LONG);
}
addkeywordto_head(tab, "EQUINOX ", "Mean equinox");
fitswrite(tab->headbuf, "EQUINOX ", &wcs->equinox, H_FLOAT, T_DOUBLE);
addkeywordto_head(tab, "RADECSYS", "Astrometric system");
switch(wcs->radecsys)
{
case RDSYS_ICRS:
fitswrite(tab->headbuf, "RADECSYS", "ICRS", H_STRING, T_STRING);
break;
case RDSYS_FK5:
fitswrite(tab->headbuf, "RADECSYS", "FK5", H_STRING, T_STRING);
break;
case RDSYS_FK4:
fitswrite(tab->headbuf, "RADECSYS", "FK4", H_STRING, T_STRING);
break;
case RDSYS_FK4_NO_E:
fitswrite(tab->headbuf, "RADECSYS", "FK4-NO-E", H_STRING, T_STRING);
break;
case RDSYS_GAPPT:
fitswrite(tab->headbuf, "RADECSYS", "GAPPT", H_STRING, T_STRING);
break;
default:
error(EXIT_FAILURE, "*Error*: unknown RADECSYS type in write_wcs()", "");
}
for (l=0; l<naxis; l++)
{
sprintf(str, "CTYPE%-3d", l+1);
addkeywordto_head(tab, str, "WCS projection type for this axis");
fitswrite(tab->headbuf, str, wcs->ctype[l], H_STRING, T_STRING);
sprintf(str, "CUNIT%-3d", l+1);
addkeywordto_head(tab, str, "Axis unit");
fitswrite(tab->headbuf, str, wcs->cunit[l], H_STRING, T_STRING);
sprintf(str, "CRVAL%-3d", l+1);
addkeywordto_head(tab, str, "World coordinate on this axis");
fitswrite(tab->headbuf, str, &wcs->crval[l], H_EXPO, T_DOUBLE);
sprintf(str, "CRPIX%-3d", l+1);
addkeywordto_head(tab, str, "Reference pixel on this axis");
fitswrite(tab->headbuf, str, &wcs->crpix[l], H_EXPO, T_DOUBLE);
for (j=0; j<naxis; j++)
{
sprintf(str, "CD%d_%d", l+1, j+1);
addkeywordto_head(tab, str, "Linear projection matrix");
fitswrite(tab->headbuf, str, &wcs->cd[l*naxis+j], H_EXPO, T_DOUBLE);
}
for (j=0; j<100; j++)
if (wcs->projp[j+100*l] != 0.0)
{
sprintf(str, "PV%d_%d", l+1, j);
addkeywordto_head(tab, str, "Projection distortion parameter");
fitswrite(tab->headbuf, str, &wcs->projp[j+100*l], H_EXPO, T_DOUBLE);
}
}
/* Update the tab data */
readbasic_head(tab);
return;
}
/******* end_wcs **************************************************************
PROTO void end_wcs(wcsstruct *wcs)
PURPOSE Free WCS (World Coordinate System) infos.
INPUT WCS structure.
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 24/05/2000
***/
void end_wcs(wcsstruct *wcs)
{
if (wcs)
{
if (wcs->lin)
{
free(wcs->lin->cdelt);
free(wcs->lin->crpix);
free(wcs->lin->pc);
free(wcs->lin->piximg);
free(wcs->lin->imgpix);
free(wcs->lin);
}
free(wcs->cel);
free(wcs->prj);
free(wcs->wcsprm);
free_tnxaxis(wcs->tnx_lngcor);
free_tnxaxis(wcs->tnx_latcor);
poly_end(wcs->inv_x);
poly_end(wcs->inv_y);
free(wcs->projp);
free(wcs);
}
return;
}
/******* wcs_supproj *********************************************************
PROTO int wcs_supproj(char *name)
PURPOSE Tell if a projection system is supported or not.
INPUT Proposed projection code name.
OUTPUT RETURN_OK if projection is supported, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 24/05/2000
***/
int wcs_supproj(char *name)
{
char projcode[26][5] =
{"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", "NONE"};
int i;
for (i=0; i<26; i++)
if (!strcmp(name, projcode[i]))
return RETURN_OK;
return RETURN_ERROR;
}
/******* invert_wcs ***********************************************************
PROTO void invert_wcs(wcsstruct *wcs)
PURPOSE Invert WCS projection mapping (using a polynomial).
INPUT WCS structure.
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 06/11/2003
***/
void invert_wcs(wcsstruct *wcs)
{
polystruct *poly;
double pixin[NAXIS],raw[NAXIS],rawmin[NAXIS];
double *outpos,*outpost, *lngpos,*lngpost,
*latpos,*latpost,
lngstep,latstep, rawsize, epsilon;
int group[] = {1,1};
/* Don't ask, this is needed by poly_init()! */
int i,j,lng,lat,deg, tnxflag, maxflag;
/* Check first that inversion is not straightforward */
lng = wcs->wcsprm->lng;
lat = wcs->wcsprm->lat;
if (!strcmp(wcs->wcsprm->pcode, "TNX"))
tnxflag = 1;
else if (!strcmp(wcs->wcsprm->pcode, "TAN")
&& (wcs->projp[1+lng*100] || wcs->projp[1+lat*100]))
tnxflag = 0;
else
return;
/* We define x as "longitude" and y as "latitude" projections */
/* We assume that PCxx cross-terms with additional dimensions are small */
/* Sample the whole image with a regular grid */
lngstep = wcs->naxisn[lng]/(WCS_NGRIDPOINTS-1.0);
latstep = wcs->naxisn[lat]/(WCS_NGRIDPOINTS-1.0);
QMALLOC(outpos, double, 2*WCS_NGRIDPOINTS2);
QMALLOC(lngpos, double, WCS_NGRIDPOINTS2);
QMALLOC(latpos, double, WCS_NGRIDPOINTS2);
for (i=0; i<wcs->naxis; i++)
raw[i] = rawmin[i] = 0.5;
outpost = outpos;
lngpost = lngpos;
latpost = latpos;
for (j=WCS_NGRIDPOINTS; j--; raw[lat]+=latstep)
{
raw[lng] = rawmin[lng];
for (i=WCS_NGRIDPOINTS; i--; raw[lng]+=lngstep)
{
if (linrev(raw, wcs->lin, pixin))
error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
wcs->wcsprm->pcode);
*(lngpost++) = pixin[lng];
*(latpost++) = pixin[lat];
if (tnxflag)
{
*(outpost++) = pixin[lng]
+raw_to_tnxaxis(wcs->tnx_lngcor,pixin[lng],pixin[lat]);
*(outpost++) = pixin[lat]
+raw_to_tnxaxis(wcs->tnx_latcor,pixin[lng],pixin[lat]);
}
else
{
raw_to_pv(wcs->prj,pixin[lng],pixin[lat], outpost, outpost+1);
outpost += 2;
}
}
}
/* Invert "longitude" */
/* Compute the extent of the pixel in reduced projected coordinates */
linrev(rawmin, wcs->lin, pixin);
pixin[lng] += ARCSEC/DEG;
linfwd(pixin, wcs->lin, raw);
rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
+(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*DEG/ARCSEC;
if (!rawsize)
error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
wcs->wcsprm->pcode);
epsilon = WCS_INVACCURACY/rawsize;
/* Find the lowest degree polynom */
poly = NULL; /* to avoid gcc -Wall warnings */
maxflag = 1;
for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++)
{
if (deg>1)
poly_end(poly);
poly = poly_init(group, 2, &deg, 1);
poly_fit(poly, outpos, lngpos, NULL, WCS_NGRIDPOINTS2, NULL);
maxflag = 0;
outpost = outpos;
lngpost = lngpos;
for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
if (fabs(poly_func(poly, outpost)-*(lngpost++))>epsilon)
{
maxflag = 1;
break;
}
}
if (maxflag)
warning("Significant inaccuracy likely to occur in projection","");
/* Now link the created structure */
wcs->prj->inv_x = wcs->inv_x = poly;
/* Invert "latitude" */
/* Compute the extent of the pixel in reduced projected coordinates */
linrev(rawmin, wcs->lin, pixin);
pixin[lat] += ARCSEC/DEG;
linfwd(pixin, wcs->lin, raw);
rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
+(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*DEG/ARCSEC;
if (!rawsize)
error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
wcs->wcsprm->pcode);
epsilon = WCS_INVACCURACY/rawsize;
/* Find the lowest degree polynom */
maxflag = 1;
for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++)
{
if (deg>1)
poly_end(poly);
poly = poly_init(group, 2, &deg, 1);
poly_fit(poly, outpos, latpos, NULL, WCS_NGRIDPOINTS2, NULL);
maxflag = 0;
outpost = outpos;
latpost = latpos;
for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
if (fabs(poly_func(poly, outpost)-*(latpost++))>epsilon)
{
maxflag = 1;
break;
}
}
if (maxflag)
warning("Significant inaccuracy likely to occur in projection","");
/* Now link the created structure */
wcs->prj->inv_y = wcs->inv_y = poly;
/* Free memory */
free(outpos);
free(lngpos);
free(latpos);
return;
}
/******* range_wcs ***********************************************************
PROTO void range_wcs(wcsstruct *wcs)
PURPOSE Find roughly the range of WCS coordinates on all axes,
and typical pixel scales.
INPUT WCS structure.
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 09/08/2006
***/
void range_wcs(wcsstruct *wcs)
{
double step[NAXIS], raw[NAXIS], rawmin[NAXIS],
world[NAXIS], world2[NAXIS];
double *worldmin, *worldmax, *scale, *worldc,
rad, radmax, lc;
int linecount[NAXIS];
int i,j, naxis, npoints, lng,lat;
naxis = wcs->naxis;
/* World range */
npoints = 1;
worldmin = wcs->wcsmin;
worldmax = wcs->wcsmax;
/* First, find the center and use it as a reference point for lng */
lng = wcs->lng;
lat = wcs->lat;
for (i=0; i<naxis; i++)
raw[i] = (wcs->naxisn[i]+1.0)/2.0;
if (raw_to_wcs(wcs, raw, world))
{
/*-- Oops no mapping there! So explore the image in an increasingly large */
/*-- domain to find a better "center" (now we know there must be angular */
/*-- coordinates) */
for (j=0; j<100; j++)
{
for (i=0; i<naxis; i++)
raw[i] += wcs->naxisn[i]/100.0*(0.5-(double)rand()/RAND_MAX);
if (!raw_to_wcs(wcs, raw, world))
break;
}
}
if (lng!=lat)
lc = fmod(world[lng]+180.0, 360.0);
else
{
lc = 0.0; /* to avoid gcc -Wall warnings */
lng = -1;
}
/* Pixel scales at image center */
scale = wcs->wcsscale;
for (i=0; i<naxis; i++)
{
if ((i==lng || i==lat) && lng!=lat)
wcs->pixscale = scale[i] = sqrt(wcs_scale(wcs, raw));
else
{
raw[i] += 1.0;
raw_to_wcs(wcs, raw, world2);
scale[i] = fabs(world2[i] - world[i]);
raw[i] -= 1.0;
if (lng==lat)
wcs->pixscale = scale[i];
}
wcs->wcsscalepos[i] = world[i];
}
/* Find "World limits" */
for (i=0; i<naxis; i++)
{
raw[i] = rawmin[i] = 0.5;
step[i] = wcs->naxisn[i]/(WCS_NRANGEPOINTS-1.0);
npoints *= WCS_NRANGEPOINTS;
worldmax[i] = -(worldmin[i] = 1e31);
linecount[i] = 0;
}
radmax = 0.0;
worldc = wcs->wcsscalepos;
for (j=npoints; j--;)
{
raw_to_wcs(wcs, raw, world);
for (i=0; i<naxis; i++)
{
/*---- Handle longitudes around 0 */
if (i==lng && world[i]>lc)
world[i] -= 359.9999;
if (world[i]<worldmin[i])
worldmin[i] = world[i];
if (world[i]>worldmax[i])
worldmax[i] = world[i];
}
/*-- Compute maximum distance to center */
if ((rad=wcs_dist(wcs, world, worldc)) > radmax)
radmax = rad;
for (i=0; i<naxis; i++)
{
raw[i] += step[i];
if (++linecount[i]<WCS_NRANGEPOINTS)
break;
else
{
linecount[i] = 0; /* No need to initialize it to 0! */
raw[i] = rawmin[i];
}
}
}
wcs->wcsmaxradius = radmax;
if (lng!=lat)
{
if (worldmax[lat]<-90.0)
worldmax[lat] = -90.0;
if (worldmax[lat]>90.0)
worldmax[lat] = 90.0;
}
return;
}
/******* frame_wcs ***********************************************************
PROTO void frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
PURPOSE Find the x and y limits of an input frame in an output image.
INPUT WCS structure of the input frame,
WCS structure of the output frame.
OUTPUT -.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 29/12/2004
***/
void frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
{
double rawin[NAXIS], rawout[NAXIS], world[NAXIS];
int linecount[NAXIS];
double worldc;
int *min, *max,
i,j, naxis, npoints, out, swapflag;
naxis = wcsin->naxis;
/* World range */
npoints = 1;
min = wcsin->outmin;
max = wcsin->outmax;
for (i=0; i<naxis; i++)
{
rawin[i] = 0.5; /* Lower pixel limits */
npoints *= WCS_NRANGEPOINTS;
max[i] = -(min[i] = 1<<30);
linecount[i] = 0;
}
/* Check if lng and lat are swapped between in and out wcs (vicious idea!) */
swapflag = (((wcsin->lng != wcsout->lng) || (wcsin->lat != wcsout->lat))
&& (wcsin->lng != wcsin->lat) && (wcsout->lng != wcsout->lat));
for (j=npoints; j--;)
{
if (!raw_to_wcs(wcsin, rawin, world))
{
if (swapflag)
{
worldc = world[wcsout->lat];
world[wcsout->lat] = world[wcsin->lat];
world[wcsin->lat] = worldc;
}
if (!wcs_to_raw(wcsout, world, rawout))
for (i=0; i<naxis; i++)
{
if ((out=(int)(rawout[i]+0.499))<min[i])
min[i] = out;
if (out>max[i])
max[i] = out;
}
}
for (i=0; i<naxis; i++)
{
rawin[i] = 0.5 + 0.5*wcsin->naxisn[i]
*(1-cos(PI*(linecount[i]+1.0)/(WCS_NRANGEPOINTS-1)));
if (++linecount[i]<WCS_NRANGEPOINTS)
break;
else
{
linecount[i] = 0; /* No need to initialize it to 0! */
rawin[i] = 0.5;
}
}
}
/* Just add a little margin, in case of... */
for (i=0; i<naxis; i++)
{
if (min[i]>-2147483647)
min[i] -= 2;
if (max[i]>2147483647)
max[i] += 2;
}
return;
}
/******* reaxe_wcs ***********************************************************
PROTO int reaxe_wcs(wcsstruct *wcs, int lng, int lat)
PURPOSE Reformulate a wcs structure to match lng and lat axis indices
INPUT WCS structure,
longitude index,
latitude index.
OUTPUT RETURN_OK if successful, RETURN_ERROR otherwise.
NOTES .
AUTHOR E. Bertin (IAP)
VERSION 20/11/2003
***/
int reaxe_wcs(wcsstruct *wcs, int lng, int lat)
{
char strlng[80], strlat[80];
double dlng,dlat,dlng2,dlat2;
int l, ilng,ilat,olng,olat, naxis;
olng = wcs->lng;
olat = wcs->lat;
if (lng<0 || lat<0 || olng<0 || olat<0)
return RETURN_ERROR;
ilng = wcs->naxisn[olng];
ilat = wcs->naxisn[olat];
wcs->naxisn[lng] = ilng;
wcs->naxisn[lat] = ilat;
strcpy(strlng, wcs->ctype[olng]);
strcpy(strlat, wcs->ctype[olat]);
strcpy(wcs->ctype[lng], strlng);
strcpy(wcs->ctype[lat], strlat);
dlng = wcs->crval[olng];
dlat = wcs->crval[olat];
wcs->crval[lng] = dlng;
wcs->crval[lat] = dlat;
naxis = wcs->naxis;
dlng = wcs->cd[olng+olng*naxis];
dlng2 = wcs->cd[olng+olat*naxis];
dlat = wcs->cd[olat+olat*naxis];
dlat2 = wcs->cd[olat+olng*naxis];
wcs->cd[lng+lng*naxis] = dlng2;
wcs->cd[lng+lat*naxis] = dlng;
wcs->cd[lat+lat*naxis] = dlat2;
wcs->cd[lat+lng*naxis] = dlat;
for (l=0; l<100; l++)
{
dlng = wcs->projp[l+olng*100];
dlat = wcs->projp[l+olat*100];
wcs->projp[l+lng*100] = dlng;
wcs->projp[l+lat*100] = dlat;
}
/*-- Reinitialize wcs */
wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);
/*-- Initialize other WCS structures */
init_wcs(wcs);
/*-- Find the range of coordinates */
range_wcs(wcs);
return RETURN_OK;
}
/******* celsys_to_eq *********************************************************
PROTO int celsys_to_eq(wcsstruct *wcs, double *wcspos)
PURPOSE Convert arbitrary celestial coordinates to equatorial.
INPUT WCS structure,
Coordinate vector.
OUTPUT RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 08/02/2007
***/
int celsys_to_eq(wcsstruct *wcs, double *wcspos)
{
double *mat,
a2,d2,sd2,cd2cp,sd,x,y;
int lng, lat;
mat = wcs->celsysmat;
a2 = wcspos[lng = wcs->wcsprm->lng]*DEG - mat[1];
d2 = wcspos[lat = wcs->wcsprm->lat]*DEG;
/* A bit of spherical trigonometry... */
/* Compute the latitude... */
sd2 = sin(d2);
cd2cp = cos(d2)*mat[2];
sd = sd2*mat[3]-cd2cp*cos(a2);
/* ...and the longitude */
y = cd2cp*sin(a2);
x = sd2 - sd*mat[3];
wcspos[lng] = fmod((atan2(y,x) + mat[0])/DEG+360.0, 360.0);
wcspos[lat] = asin(sd)/DEG;
return RETURN_OK;
}
/******* eq_to_celsys *********************************************************
PROTO int eq_to_celsys(wcsstruct *wcs, double *wcspos)
PURPOSE Convert equatorial to arbitrary celestial coordinates.
INPUT WCS structure,
Coordinate vector.
OUTPUT RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 08/02/2007
***/
int eq_to_celsys(wcsstruct *wcs, double *wcspos)
{
double *mat,
a,d,sd2,cdcp,sd,x,y;
int lng, lat;
mat = wcs->celsysmat;
a = wcspos[lng = wcs->wcsprm->lng]*DEG - mat[0];
d = wcspos[lat = wcs->wcsprm->lat]*DEG;
/* A bit of spherical trigonometry... */
/* Compute the latitude... */
sd = sin(d);
cdcp = cos(d)*mat[2];
sd2 = sd*mat[3]+cdcp*cos(a);
/* ...and the longitude */
y = cdcp*sin(a);
x = sd2*mat[3]-sd;
wcspos[lng] = fmod((atan2(y,x) + mat[1])/DEG+360.0, 360.0);
wcspos[lat] = asin(sd2)/DEG;
return RETURN_OK;
}
/******* raw_to_wcs ***********************************************************
PROTO int raw_to_wcs(wcsstruct *, double *, double *)
PURPOSE Convert raw (pixel) coordinates to WCS (World Coordinate System).
INPUT WCS structure,
Pointer to the array of input coordinates,
Pointer to the array of output coordinates.
OUTPUT RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 08/02/2007
***/
int raw_to_wcs(wcsstruct *wcs, double *pixpos, double *wcspos)
{
double imgcrd[NAXIS],
phi,theta;
int i;
if (wcsrev((const char(*)[9])wcs->ctype, wcs->wcsprm, pixpos,
wcs->lin,imgcrd, wcs->prj, &phi, &theta, wcs->crval, wcs->cel, wcspos))
{
for (i=0; i<wcs->naxis; i++)
wcspos[i] = WCS_NOCOORD;
return RETURN_ERROR;
}
/* If needed, convert from a different coordinate system to equatorial */
if (wcs->celsysconvflag)
celsys_to_eq(wcs, wcspos);
return RETURN_OK;
}
/******* wcs_to_raw ***********************************************************
PROTO int wcs_to_raw(wcsstruct *, double *, double *)
PURPOSE Convert WCS (World Coordinate System) coords to raw (pixel) coords.
INPUT WCS structure,
Pointer to the array of input coordinates,
Pointer to the array of output coordinates.
OUTPUT RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 08/02/2007
***/
int wcs_to_raw(wcsstruct *wcs, double *wcspos, double *pixpos)
{
double imgcrd[NAXIS],
phi,theta;
int i;
/* If needed, convert to a coordinate system different from equatorial */
if (wcs->celsysconvflag)
eq_to_celsys(wcs, wcspos);
if (wcsfwd((const char(*)[9])wcs->ctype,wcs->wcsprm,wcspos,
wcs->crval, wcs->cel,&phi,&theta,wcs->prj, imgcrd,wcs->lin,pixpos))
{
for (i=0; i<wcs->naxis; i++)
pixpos[i] = WCS_NOCOORD;
return RETURN_ERROR;
}
return RETURN_OK;
}
/******* red_to_raw **********************************************************
PROTO int red_to_raw(wcsstruct *, double *, double *)
PURPOSE Convert reduced (World Coordinate System) coords to raw (pixel)
coords.
INPUT WCS structure,
Pointer to the array of input (reduced) coordinates,
Pointer to the array of output (pixel) coordinates.
OUTPUT RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 23/10/2003
***/
int red_to_raw(wcsstruct *wcs, double *redpos, double *pixpos)
{
struct wcsprm *wcsprm;
double offset;
wcsprm = wcs->wcsprm;
/* Initialize if required */
if (wcsprm && wcsprm->flag != WCSSET)
{
if (wcsset(wcs->naxis, (const char(*)[9])wcs->ctype, wcsprm))
return RETURN_ERROR;
}
if (wcsprm && wcsprm->flag != 999 && wcsprm->cubeface != -1)
{
/*-- Separation between faces */
offset = (wcs->prj->r0 == 0.0 ? 90.0 : wcs->prj->r0*PI/2.0);
/*-- Stack faces in a cube */
if (redpos[wcs->lat] < -0.5*offset)
{
redpos[wcs->lat] += offset;
redpos[wcsprm->cubeface] = 5.0;
}
else if (redpos[wcs->lat] > 0.5*offset)
{
redpos[wcs->lat] -= offset;
redpos[wcsprm->cubeface] = 0.0;
}
else if (redpos[wcs->lng] > 2.5*offset)
{
redpos[wcs->lng] -= 3.0*offset;
redpos[wcsprm->cubeface] = 4.0;
}
else if (redpos[wcs->lng] > 1.5*offset)
{
redpos[wcs->lng] -= 2.0*offset;
redpos[wcsprm->cubeface] = 3.0;
}
else if (redpos[wcs->lng] > 0.5*offset)
{
redpos[wcs->lng] -= offset;
redpos[wcsprm->cubeface] = 2.0;
}
else
redpos[wcsprm->cubeface] = 1.0;
}
/* Apply forward linear transformation */
if (linfwd(redpos, wcs->lin, pixpos))
return RETURN_ERROR;
return RETURN_OK;
}
/******* raw_to_red **********************************************************
PROTO int raw_to_red(wcsstruct *, double *, double *)
PURPOSE Convert raw (pixel) coordinates to reduced WCS coordinates.
INPUT WCS structure,
Pointer to the array of input (pixel) coordinates,
Pointer to the array of output (reduced) coordinates.
OUTPUT RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 23/10/2003
***/
int raw_to_red(wcsstruct *wcs, double *pixpos, double *redpos)
{
struct wcsprm *wcsprm;
double offset;
int face;
wcsprm = wcs->wcsprm;
/* Initialize if required */
if (wcsprm && wcsprm->flag != WCSSET)
{
if (wcsset(wcs->naxis, (const char(*)[9])wcs->ctype, wcsprm))
return RETURN_ERROR;
}
/* Apply reverse linear transformation */
if (linrev(pixpos, wcs->lin, redpos))
return RETURN_ERROR;
if (wcsprm && wcsprm->flag != 999 && wcsprm->cubeface != -1)
{
/*-- Do we have a CUBEFACE axis? */
face = (int)(redpos[wcsprm->cubeface] + 0.5);
if (fabs(redpos[wcsprm->cubeface]-face) > 1e-10)
return RETURN_ERROR;
/*-- Separation between faces. */
offset = (wcs->prj->r0 == 0.0 ? 90.0 : wcs->prj->r0*PI/2.0);
/*-- Lay out faces in a plane. */
switch (face)
{
case 0:
redpos[wcs->lat] += offset;
break;
case 1:
break;
case 2:
redpos[wcs->lng] += offset;
break;
case 3:
redpos[wcs->lng] += offset*2;
break;
case 4:
redpos[wcs->lng] += offset*3;
break;
case 5:
redpos[wcs->lat] -= offset;
break;
default:
return RETURN_ERROR;
}
}
return RETURN_OK;
}
/******* wcs_dist ***********************************************************
PROTO double wcs_dist(wcsstruct *wcs, double *wcspos1, double *wcspos2)
PURPOSE Compute the angular distance between 2 points on the sky.
INPUT WCS structure,
Pointer to the first array of world coordinates,
Pointer to the second array of world coordinates.
OUTPUT Angular distance (in degrees) between points.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 24/07/2002
***/
double wcs_dist(wcsstruct *wcs, double *wcspos1, double *wcspos2)
{
double d, dp;
int i, lng, lat;
lng = wcs->lng;
lat = wcs->lat;
if (lat!=lng)
{
/*-- We are operating in angular coordinates */
d = sin(wcspos1[lat]*DEG)*sin(wcspos2[lat]*DEG)
+ cos(wcspos1[lat]*DEG)*cos(wcspos2[lat]*DEG)
*cos((wcspos1[lng]-wcspos2[lng])*DEG);
return d>-1.0? (d<1.0 ? acos(d)/DEG : 0.0) : 180.0;
}
else
{
d = 0.0;
for (i=0; i<wcs->naxis; i++)
{
dp = wcspos1[i] - wcspos2[i];
d += dp*dp;
}
return sqrt(d);
}
}
/******* wcs_scale ***********************************************************
PROTO double wcs_scale(wcsstruct *wcs, double *pixpos)
PURPOSE Compute the sky area equivalent to a local pixel.
INPUT WCS structure,
Pointer to the array of local raw coordinates,
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 03/01/2008
***/
double wcs_scale(wcsstruct *wcs, double *pixpos)
{
double wcspos[NAXIS], wcspos1[NAXIS], wcspos2[NAXIS], pixpos2[NAXIS];
double dpos1,dpos2;
int lng, lat;
if (raw_to_wcs(wcs, pixpos, wcspos))
return 0.0;
lng = wcs->lng;
lat = wcs->lat;
if (lng == lat)
{
lng = 0;
lat = 1;
}
/* Compute pixel scale */
pixpos2[lng] = pixpos[lng] + 1.0;
pixpos2[lat] = pixpos[lat];
if (raw_to_wcs(wcs, pixpos2, wcspos1))
return 0.0;
pixpos2[lng] -= 1.0;
pixpos2[lat] += 1.0;
if (raw_to_wcs(wcs, pixpos2, wcspos2))
return 0.0;
dpos1 = wcspos1[lng]-wcspos[lng];
dpos2 = wcspos2[lng]-wcspos[lng];
if (wcs->lng!=wcs->lat)
{
if (dpos1>180.0)
dpos1 -= 360.0;
else if (dpos1<-180.0)
dpos1 += 360.0;
if (dpos2>180.0)
dpos2 -= 360.0;
else if (dpos2<-180.0)
dpos2 += 360.0;
return fabs((dpos1*(wcspos2[lat]-wcspos[lat])
-(wcspos1[lat]-wcspos[lat])*dpos2)*cos(wcspos[lat]*DEG));
}
else
return fabs((dpos1*(wcspos2[lat]-wcspos[lat])
-(wcspos1[lat]-wcspos[lat])*dpos2));
}
/****** wcs jacobian *********************************************************
PROTO double wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
PURPOSE Compute the local Jacobian matrice of the astrometric deprojection.
INPUT WCS structure,
Pointer to the array of local raw coordinates,
Pointer to the jacobian array (output).
OUTPUT Determinant over spatial coordinates (=pixel area), or -1.0 if mapping
was unsuccesful.
NOTES Memory must have been allocated (naxis*naxis*sizeof(double)) for the
Jacobian array.
AUTHOR E. Bertin (IAP)
VERSION 11/10/2007
***/
double wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
{
double pixpos0[NAXIS], wcspos0[NAXIS], wcspos[NAXIS],
dpos;
int i,j, lng,lat,naxis;
lng = wcs->lng;
lat = wcs->lat;
naxis = wcs->naxis;
for (i=0; i<naxis; i++)
pixpos0[i] = pixpos[i];
if (raw_to_wcs(wcs, pixpos0, wcspos0) == RETURN_ERROR)
return -1.0;
for (i=0; i<naxis; i++)
{
pixpos0[i] += 1.0;
if (raw_to_wcs(wcs, pixpos0, wcspos) == RETURN_ERROR)
return -1.0;
pixpos0[i] -= 1.0;
for (j=0; j<naxis; j++)
{
dpos = wcspos[j]-wcspos0[j];
if (lng!=lat && j==lng)
{
if (dpos>180.0)
dpos -= 360.0;
else if (dpos<-180.0)
dpos += 360.0;
dpos *= cos(wcspos0[lat]*DEG);
}
jacob[j*naxis+i] = dpos;
}
}
if (lng==lat)
{
lng = 0;
lat = 1;
}
return fabs(jacob[lng+naxis*lng]*jacob[lat+naxis*lat]
- jacob[lat+naxis*lng]*jacob[lng+naxis*lat]);
}
/******* wcs_chirality *******************************************************
PROTO int wcs_chirality(wcsstruct *wcs)
PURPOSE Compute the chirality of a WCS projection.
INPUT WCS structure.
OUTPUT +1 if determinant of matrix is positive, -1 if negative, 0 if null.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 26/09/2006
***/
int wcs_chirality(wcsstruct *wcs)
{
double a;
int lng,lat, naxis;
lng = wcs->lng;
lat = wcs->lat;
naxis = wcs->naxis;
if (lng==lat && naxis>=2)
{
lng = 0;
lat = 1;
}
a = wcs->cd[lng*naxis+lng]*wcs->cd[lat*naxis+lat]
- wcs->cd[lng*naxis+lat]*wcs->cd[lat*naxis+lng];
return a>TINY? 1 : (a<-TINY? -1 : 0);
}
/****** precess_wcs **********************************************************
PROTO void precess_wcs(wcsstruct *wcs, double yearin, double yearout)
PURPOSE Precess the content of a WCS structure according to the equinox.
INPUT WCS structure,
Input year,
Output year.
OUTPUT -.
NOTES Epoch for coordinates should be J2000 (FK5 system).
AUTHOR E. Bertin (IAP)
VERSION 04/01/2008
***/
void precess_wcs(wcsstruct *wcs, double yearin, double yearout)
{
double crval[NAXIS],a[NAXIS*NAXIS],b[NAXIS*NAXIS],
*c,*at,
val, cas, sas, angle, dalpha;
int i,j,k, lng,lat, naxis;
lng = wcs->lng;
lat = wcs->lat;
if (lat==lng || yearin==yearout)
return;
naxis = wcs->naxis;
/* Precess to year out */
precess(yearin, wcs->crval[lng], wcs->crval[lat], yearout,
&crval[lng], &crval[lat]);
dalpha = (crval[lng] - wcs->crval[lng])*DEG;
/* Compute difference angle with the north axis between start and end */
angle = (dalpha!=0.0 && (crval[lat] - wcs->crval[lat])*DEG != 0.0) ?
180.0 - (atan2(sin(dalpha),
cos(crval[lat]*DEG)*tan(wcs->crval[lat]*DEG)
- sin(crval[lat]*DEG)*cos(dalpha))
+ atan2(sin(dalpha),
cos(wcs->crval[lat]*DEG)*tan(crval[lat]*DEG)
- sin(wcs->crval[lat]*DEG)*cos(dalpha)))/DEG
: 0.0;
/* A = C*B */
c = wcs->cd;
/* The B matrix is made of 2 numbers */
cas = cos(angle*DEG);
sas = sin(-angle*DEG);
for (i=0; i<naxis; i++)
b[i+i*naxis] = 1.0;
b[lng+lng*naxis] = cas;
b[lat+lng*naxis] = -sas;
b[lng+lat*naxis] = sas;
b[lat+lat*naxis] = cas;
at = a;
for (j=0; j<naxis; j++)
for (i=0; i<naxis; i++)
{
val = 0.0;
for (k=0; k<naxis; k++)
val += c[k+j*naxis]*b[i+k*naxis];
*(at++) = val;
}
at = a;
for (i=0; i<naxis*naxis; i++)
*(c++) = *(at++);
wcs->crval[lng] = crval[lng];
wcs->crval[lat] = crval[lat];
wcs->equinox = yearout;
/* Initialize other WCS structures */
init_wcs(wcs);
/* Find the range of coordinates */
range_wcs(wcs);
/* Invert projection corrections */
invert_wcs(wcs);
return;
}
/********************************* precess ***********************************/
/*
precess equatorial coordinates according to the equinox (from Ephemerides du
Bureau des Longitudes 1992). Epoch for coordinates should be J2000
(FK5 system).
*/
void precess(double yearin, double alphain, double deltain,
double yearout, double *alphaout, double *deltaout)
{
double dzeta,theta,z, t1,t1t1, t2,t2t2,t2t2t2,
cddsadz, cddcadz, cdd, sdd, adz, cdin,sdin,ct,st,caindz;
alphain *= DEG;
deltain *= DEG;
t1 = (yearin - 2000.0)/1000.0;
t2 = (yearout - yearin)/1000.0;
t1t1 = t1*t1;
t2t2t2 = (t2t2 = t2*t2)*t2;
theta = (97171.735e-06 - 413.691e-06*t1 - 1.052e-06 * t1t1) * t2
+ (-206.846e-06 - 1.052e-06*t1) * t2t2 - 202.812e-06 * t2t2t2;
dzeta = (111808.609e-06 + 677.071e-06*t1 - 0.674e-06 * t1t1) * t2
+ (146.356e-06 - 1.673e-06*t1) * t2t2 + 87.257e-06 * t2t2t2;
z = (111808.609e-06 +677.071e-06*t1 - 0.674e-06 * t1t1) * t2
+ (530.716e-06 + 0.320e-06*t1) * t2t2 + 88.251e-06 * t2t2t2;
cddsadz = (cdin=cos(deltain)) * sin(alphain+dzeta);
cddcadz = -(sdin=sin(deltain))*(st=sin(theta))
+cdin*(ct=cos(theta))*(caindz=cos(alphain+dzeta));
sdd = sdin*ct + cdin*st*caindz;
cdd = cos(*deltaout = asin(sdd));
adz = asin(cddsadz/cdd);
if (cddcadz<0.0)
adz = PI - adz;
if (adz<0.0)
adz += 2.0*PI;
adz += z;
*alphaout = adz/DEG;
*deltaout /= DEG;
return;
}
/********************************* b2j ***********************************/
/*
conver equatorial coordinates from equinox and epoch B1950 to equinox and
epoch J2000 for extragalactic sources (from Aoki et al. 1983).
*/
void b2j(double yearobs, double alphain, double deltain,
double *alphaout, double *deltaout)
{
int i,j;
double a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6},
ap[3] = {1.245e-3, -1.580e-3, -0.659e-3},
m[6][6] = {
{ 0.9999256782, -0.0111820611, -0.0048579477,
0.00000242395018, -0.00000002710663, -0.00000001177656},
{ 0.0111820610, 0.9999374784, -0.0000271765,
0.00000002710663, 0.00000242397878, -0.00000000006587},
{ 0.0048579479, -0.0000271474, 0.9999881997,
0.00000001177656, -0.00000000006582, 0.00000242410173},
{-0.000551, -0.238565, 0.435739,
0.99994704, -0.01118251, -0.00485767},
{ 0.238514, -0.002662, -0.008541,
0.01118251, 0.99995883, -0.00002718},
{-0.435623, 0.012254, 0.002117,
0.00485767, -0.00002714, 1.00000956}},
a1[3], r[3], ro[3], r1[3], r2[3], v1[3], v[3];
double cai, sai, cdi, sdi, dotp, rmod, alpha, delta,
t1 = (yearobs - 1950.0)/100.0;
alphain *= PI/180.0;
deltain *= PI/180.0;
cai = cos(alphain);
sai = sin(alphain);
cdi = cos(deltain);
sdi = sin(deltain);
ro[0] = cdi*cai;
ro[1] = cdi*sai;
ro[2] = sdi;
dotp = 0.0;
for (i=0; i<3; i++)
{
a1[i] = a[i]+ap[i]*ARCSEC*t1;
dotp += a1[i]*ro[i];
}
for (i=0; i<3; i++)
{
r1[i] = ro[i] - a1[i] + dotp*ro[i];
r[i] = v[i] = v1[i] = 0.0;
}
for (j=0; j<6; j++)
for (i=0; i<6; i++)
{
if (j<3)
r[j] += m[j][i]*(i<3?r1[i]:v1[i-3]);
else
v[j-3] += m[j][i]*(i<3?r1[i]:v1[i-3]);
}
rmod = 0.0;
for (i=0; i<3; i++)
{
r2[i] = r[i]+v[i]*ARCSEC*(t1-0.5);
rmod += r2[i]*r2[i];
}
rmod = sqrt(rmod);
delta = asin(r2[2]/rmod);
alpha = acos(r2[0]/cos(delta)/rmod);
if (r2[1]<0)
alpha = 2*PI - alpha;
*alphaout = alpha*180.0/PI;
*deltaout = delta*180.0/PI;
return;
}
/*********************************** j2b *************************************/
/*
conver equatorial coordinates from equinox and epoch J2000 to equinox and
epoch B1950 for extragalactic sources (from Aoki et al. 1983, after
inversion of their matrix and some custom arrangements).
*/
void j2b(double yearobs, double alphain, double deltain,
double *alphaout, double *deltaout)
{
int i,j;
double a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6},
ap[3] = {1.245e-3, -1.580e-3, -0.659e-3},
m[6][6] = {
{ 0.9999256794678425, 0.01118148281196562, 0.004859003848996022,
-2.423898417033081e-06,-2.710547600126671e-08,-1.177738063266745e-08},
{-0.01118148272969232, 0.9999374849247641, -2.717708936468247e-05,
2.710547578707874e-08,-2.423927042585208e-06, 6.588254898401055e-11},
{-0.00485900399622881, -2.715579322970546e-05, 0.999988194643078,
1.177738102358923e-08, 6.582788892816657e-11,-2.424049920613325e-06},
{-0.0005508458576414713, 0.2384844384742432, -0.4356144527773499,
0.9999043171308133, 0.01118145410120206, 0.004858518651645554},
{-0.2385354433560954, -0.002664266996872802, 0.01225282765749546,
-0.01118145417187502, 0.9999161290795875, -2.717034576263522e-05},
{ 0.4357269351676567, -0.008536768476441086, 0.002113420799663768,
-0.004858518477064975, -2.715994547222661e-05, 0.9999668385070383}},
a1[3], r[3], ro[3], r1[3], r2[3], v1[3], v[3];
double cai, sai, cdi, sdi, dotp, rmod, alpha, delta, t1;
/* Convert Julian years from J2000.0 to tropic centuries from B1950.0 */
t1 = ((yearobs - 2000.0) + (MJD2000 - MJD1950)/365.25)*JU2TROP/100.0;
alphain *= DEG;
deltain *= DEG;
cai = cos(alphain);
sai = sin(alphain);
cdi = cos(deltain);
sdi = sin(deltain);
r[0] = cdi*cai;
r[1] = cdi*sai;
r[2] = sdi;
for (i=0; i<3; i++)
v[i] = r2[i] = v1[i] = 0.0;
for (j=0; j<6; j++)
for (i=0; i<6; i++)
if (j<3)
r2[j] += m[j][i]*(i<3?r[i]:v[i-3]);
else
v1[j-3] += m[j][i]*(i<3?r[i]:v[i-3]);
for (i=0; i<3; i++)
r1[i] = r2[i]+v1[i]*ARCSEC*t1;
dotp = 0.0;
for (i=0; i<3; i++)
{
a1[i] = a[i]+ap[i]*ARCSEC*t1;
dotp += a1[i]*(r1[i]+a1[i]);
}
dotp = 2.0/(sqrt(1+4.0*dotp)+1.0);
rmod = 0.0;
for (i=0; i<3; i++)
{
ro[i] = dotp*(r1[i]+a1[i]);
rmod += ro[i]*ro[i];
}
rmod = sqrt(rmod);
delta = asin(ro[2]/rmod);
alpha = acos(ro[0]/cos(delta)/rmod);
if (ro[1]<0)
alpha = 2.0*PI - alpha;
*alphaout = alpha/DEG;
*deltaout = delta/DEG;
return;
}
/******************************** degtosexal *********************************/
/*
Convert degrees to hh mm ss.xx alpha coordinates.
*/
char *degtosexal(double alpha, char *str)
{
int hh, mm;
double ss;
if (alpha>=0.0 && alpha <360.0)
{
hh = (int)(alpha/15.0);
mm = (int)(60.0*(alpha/15.0 - hh));
ss = 60.0*(60.0*(alpha/15.0 - hh) - mm);
}
else
hh = mm = ss = 0.0;
sprintf(str,"%02d:%02d:%05.2f", hh, mm, ss);
return str;
}
/******************************** degtosexde *********************************/
/*
Convert degrees to dd dm ds.x delta coordinates.
*/
char *degtosexde(double delta, char *str)
{
char sign;
double ds;
int dd, dm;
sign = delta<0.0?'-':'+';
delta = fabs(delta);
if (delta>=-90.0 && delta <=90.0)
{
dd = (int)delta;
dm = (int)(60.0*(delta - dd));
ds = 60.0*fabs(60.0*(delta - dd) - dm);
}
else
dd = dm = ds = 0.0;
sprintf(str,"%c%02d:%02d:%04.1f", sign, dd, dm, ds);
return str;
}
/******************************** sextodegal *********************************/
/*
Convert hh mm ss.xxx alpha coordinates to degrees.
*/
double sextodegal(char *hms)
{
double val;
char *ptr;
val = atof(strtok_r(hms, ": \t", &ptr))*15.0; /* Hours */
val += atof(strtok_r(NULL, ": \t", &ptr))/4.0; /* Minutes */
val += atof(strtok_r(NULL, ": \t", &ptr))/240.0; /* Seconds */
return val;
}
/******************************** sextodegde *********************************/
/*
Convert dd dm ds.xxx delta coordinates to degrees.
*/
double sextodegde(char *dms)
{
double val, sgn;
char *str, *ptr;
str = strtok_r(dms, ": \t", &ptr);
sgn = (strchr(str, '-') ? -1.0:1.0);
val = atof(dms); /* Degrees */
val += atof(strtok_r(NULL, ": \t", &ptr))*sgn/60.0; /* Minutes */
val += atof(strtok_r(NULL, ": \t", &ptr))*sgn/3600.0; /* Seconds */
return val;
}
/******************************** fmod_0_p360 *******************************/
/*
Fold input angle in the [0,+360[ domain.
*/
double fmod_0_p360(double angle)
{
return angle>0.0? fmod(angle,360.0) : fmod(angle,360.0)+360.0;
}
/******************************** fmod_m90_p90 *******************************/
/*
Fold input angle in the [-90,+90[ domain.
*/
double fmod_m90_p90(double angle)
{
return angle>0.0? fmod(angle+90.0,180.0)-90.0 : fmod(angle-90.0,180.0)+90.0;
}
src/fitswcs.h 0 → 100644
View file @ ece261cd
/*
fitswcs.h
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: LDACTools+
*
* Author: E.BERTIN (IAP)
*
* Contents: Include file for fitswcs.c
*
* Last modify: 26/04/2008
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifndef _FITSWCS_H_
#define _FITSWCS_H_
/*-------------------------------- macros -----------------------------------*/
/*----------------------------- Internal constants --------------------------*/
#define NAXIS 2 /* Max number of FITS axes */
#define DEG (PI/180.0) /* 1 deg in radians */
#define ARCMIN (DEG/60.0) /* 1 arcsec in radians */
#define ARCSEC (DEG/3600.0) /* 1 arcsec in radians */
#define MAS (ARCSEC/1000.0) /* 1 mas in radians */
#define MJD2000 51544.50000 /* Modified Julian date for J2000.0 */
#define MJD1950 33281.92346 /* Modified Julian date for B1950.0 */
#define JU2TROP 1.0000214 /* 1 Julian century in tropical units*/
#define WCS_NOCOORD 1e31 /* Code for non-existing coordinates */
#define WCS_NGRIDPOINTS 12 /* Number of WCS grid points / axis */
#define WCS_NGRIDPOINTS2 (WCS_NGRIDPOINTS*WCS_NGRIDPOINTS)
#define WCS_INVMAXDEG 9 /* Maximum inversion polynom degree */
#define WCS_INVACCURACY 0.04 /* Maximum inversion error (pixels) */
#define WCS_NRANGEPOINTS 32 /* Number of WCS range points / axis */
/*-------------------------------- typedefs ---------------------------------*/
typedef enum {CELSYS_NATIVE, CELSYS_PIXEL, CELSYS_EQUATORIAL, CELSYS_GALACTIC,
CELSYS_ECLIPTIC, CELSYS_SUPERGALACTIC} celsysenum;
/*------------------------------- structures --------------------------------*/
typedef struct wcs
{
int naxis; /* Number of image axes */
int naxisn[NAXIS]; /* FITS NAXISx parameters */
char ctype[NAXIS][9]; /* FITS CTYPE strings */
char cunit[NAXIS][32]; /* FITS CUNIT strings */
double crval[NAXIS]; /* FITS CRVAL parameters */
double cdelt[NAXIS]; /* FITS CDELT parameters */
double crpix[NAXIS]; /* FITS CRPIX parameters */
double crder[NAXIS]; /* FITS CRDER parameters */
double csyer[NAXIS]; /* FITS CSYER parameters */
double cd[NAXIS*NAXIS]; /* FITS CD matrix */
double *projp; /* FITS PV/PROJP mapping parameters */
int nprojp; /* number of useful projp parameters */
double longpole,latpole; /* FITS LONGPOLE and LATPOLE */
double wcsmin[NAXIS]; /* minimum values of WCS coords */
double wcsmax[NAXIS]; /* maximum values of WCS coords */
double wcsscale[NAXIS]; /* typical pixel scale at center */
double wcsscalepos[NAXIS]; /* WCS coordinates of scaling point */
double wcsmaxradius; /* Maximum distance to wcsscalepos */
int outmin[NAXIS]; /* minimum output pixel coordinate */
int outmax[NAXIS]; /* maximum output pixel coordinate */
int lat,lng; /* longitude and latitude axes # */
double r0; /* projection "radius" */
double lindet; /* Determinant of the local matrix */
int chirality; /* Chirality of the CD matrix */
double pixscale; /* (Local) pixel scale */
double ap2000,dp2000; /* J2000 coordinates of pole */
double ap1950,dp1950; /* B1950 coordinates of pole */
double obsdate; /* Date of observations */
double equinox; /* Equinox of observations */
double epoch; /* Epoch of observations (deprec.) */
enum {RDSYS_ICRS, RDSYS_FK5, RDSYS_FK4, RDSYS_FK4_NO_E, RDSYS_GAPPT}
radecsys; /* FITS RADECSYS reference frame */
celsysenum celsys; /* Celestial coordinate system */
double celsysmat[4]; /* Equ. <=> Cel. system parameters */
int celsysconvflag; /* Equ. <=> Cel. conversion needed? */
struct wcsprm *wcsprm; /* WCSLIB's wcsprm structure */
struct linprm *lin; /* WCSLIB's linprm structure */
struct celprm *cel; /* WCSLIB's celprm structure */
struct prjprm *prj; /* WCSLIB's prjprm structure */
struct tnxaxis *tnx_latcor; /* IRAF's TNX latitude corrections */
struct tnxaxis *tnx_lngcor; /* IRAF's TNX longitude corrections */
struct poly *inv_x; /* Proj. correction polynom in x */
struct poly *inv_y; /* Proj. correction polynom in y */
} wcsstruct;
/*------------------------------- functions ---------------------------------*/
extern wcsstruct *create_wcs(char **ctype, double *crval, double *crpix,
double *cdelt, int *naxisn, int naxis),
*copy_wcs(wcsstruct *wcsin),
*read_wcs(tabstruct *tab);
extern double fmod_0_p360(double angle),
fmod_m90_p90(double angle),
sextodegal(char *hms),
sextodegde(char *dms),
wcs_dist(wcsstruct *wcs,
double *wcspos1, double *wcspos2),
wcs_jacobian(wcsstruct *wcs, double *pixpos,
double *jacob),
wcs_scale(wcsstruct *wcs, double *pixpos);
extern int celsys_to_eq(wcsstruct *wcs, double *wcspos),
eq_to_celsys(wcsstruct *wcs, double *wcspos),
raw_to_red(wcsstruct *wcs,
double *pixpos, double *redpos),
raw_to_wcs(wcsstruct *wcs,
double *pixpos, double *wcspos),
reaxe_wcs(wcsstruct *wcs, int lng, int lat),
red_to_raw(wcsstruct *wcs,
double *redpos, double *pixpos),
wcs_chirality(wcsstruct *wcs),
wcs_supproj(char *name),
wcs_to_raw(wcsstruct *wcs,
double *wcspos, double *pixpos);
extern char *degtosexal(double alpha, char *str),
*degtosexde(double delta, char *str);
extern void b2j(double yearobs, double alphain, double deltain,
double *alphaout, double *deltaout),
end_wcs(wcsstruct *wcs),
init_wcs(wcsstruct *wcs),
init_wcscelsys(wcsstruct *wcs),
invert_wcs(wcsstruct *wcs),
frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout),
j2b(double yearobs, double alphain, double deltain,
double *alphaout, double *deltaout),
precess(double yearin, double alphain, double deltain,
double yearout,
double *alphaout, double *deltaout),
precess_wcs(wcsstruct *wcs, double yearin,
double yearout),
range_wcs(wcsstruct *wcs),
write_wcs(tabstruct *tab, wcsstruct *wcs);
#endif
src/globals.h
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: global declarations.
*
* Last modify: 14/07/2006
* Last modify: 11/05/2008
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -27,13 +27,16 @@ float ctg[37], stg[37];
char gstr[MAXCHAR];
/*------------------------------- functions ---------------------------------*/
extern void allocparcelout(void),
extern void alloccatparams(void),
allocparcelout(void),
analyse(picstruct *, picstruct *, int, objliststruct *),
blankit(char *, int),
endcat(char *error),
reendcat(void),
changecatparamarrays(char *keyword, int *axisn, int naxis),
closecheck(void),
copydata(picstruct *, int, int),
dumpparams(void),
endfield(picstruct *),
endobject(picstruct *, picstruct *, picstruct *, picstruct *,
int, objliststruct *),
......
src/growth.c
View file @ ece261cd
......@@ -9,7 +9,7 @@
*
* Contents: Make growth curves.
*
* Last modify: 15/02/2005
* Last modify: 19/12/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
......@@ -103,7 +103,7 @@ void makeavergrowth(picstruct *field, picstruct *wfield, objstruct *obj)
pflag = (prefs.detect_type==PHOTO)? 1:0;
corrflag = (prefs.mask_type==MASK_CORRECT);
var = backnoise2 = field->backsig*field->backsig;
gain = prefs.gain;
gain = field->gain;
/* Integration radius */
rlim = GROWTH_NSIG*obj->a;
......
src/ldactoasc.c 0 → 100644
View file @ ece261cd
/*
ldactoasc.c
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: LDACtoASC
*
* Author: E.BERTIN (IAP)
*
* Contents: Convert LDAC binary format to ASCII.
*
* Last modify: 19/12/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "ldactoasc.h"
#include "fits/fitscat.h"
#define SYNTAX "ldactoasc catalog\n"
extern const char notokstr[];
/********************************** main ************************************/
int main(int argc, char *argv[])
{
catstruct *cat;
tabstruct *tab;
unsigned short ashort=1;
char catname[MAXCHAR];
int a, t, opt,opt2, flag;
if (argc<2)
{
fprintf(OUTPUT, "\n %s Version %s (%s)\n",
BANNER, MYVERSION, DATE);
fprintf(OUTPUT, "\nFor information, please contact: %s\n", COPYRIGHT);
error(EXIT_SUCCESS, "SYNTAX: ", SYNTAX);
}
/* Test if byteswapping will be needed */
bswapflag = *((char *)&ashort);
/* Default parameters */
for (a=1; a<argc; a++)
{
if (*(argv[a]) == '-')
{
opt = (int)argv[a][1];
{
opt2 = (int)tolower((int)argv[a][2]);
if (opt == '-')
{
opt = opt2;
opt2 = (int)tolower((int)argv[a][3]);
}
switch(opt)
{
case 'v':
printf("%s version %s (%s)\n", BANNER,MYVERSION,DATE);
exit(EXIT_SUCCESS);
break;
case 'h':
fprintf(OUTPUT, "\nSYNTAX: %s", SYNTAX);
exit(EXIT_SUCCESS);
break;
default:
error(EXIT_SUCCESS, "SYNTAX: ", SYNTAX);
}
}
}
else
strcpy(catname, argv[a]);
}
flag = 1; /* display banner of first extension */
if ((cat = read_cat(catname)))
{
tab = cat->tab;
for (t=cat->ntab; t--; tab=tab->nexttab)
if (!strcmp("LDAC_OBJECTS", tab->extname)
|| !strcmp("OBJECTS", tab->extname))
{
show_keys(tab, NULL, NULL, 0, NULL, stdout, 1, flag, 0, SHOW_ASCII);
flag = 0;
}
free_cat(&cat, 1);
}
else
error(EXIT_FAILURE,"Cannot open ",catname);
return EXIT_SUCCESS;
}
src/ldactoasc.h 0 → 100644
View file @ ece261cd
/*
ldactoasc.h
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
* Part of: LDACtoASC
*
* Author: E.BERTIN (IAP)
*
* Contents: global definitions.
*
* Last modify: 04/06/2007
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
/* Check if we are using a configure script here */
#ifndef HAVE_CONFIG_H
#define VERSION "1.x"
#define DATE "2007-06-04"
#define THREADS_NMAX 16 /* max. number of threads */
#endif
/*------------------------ what, who, when and where ------------------------*/
#define BANNER "LDACtoASC"
#ifdef USE_THREADS
#define MYVERSION VERSION "-MP"
#else
#define MYVERSION VERSION
#endif
#define COPYRIGHT "Emmanuel BERTIN <bertin@iap.fr>"
#define WEBSITE "http://terapix.iap.fr/soft/sextractor/"
#define INSTITUTE "TERAPIX team at IAP http://terapix.iap.fr"
/*----------------------------- Physical constants --------------------------*/
#ifndef PI
#define PI 3.1415926535898
#endif
/*----------------------------- Internal constants --------------------------*/
#define BIG 1e+30 /* a huge number */
#define TINY (1.0/BIG) /* a small number */
#define OUTPUT stdout /* where all msgs are sent */
#define MAXCHAR 512 /* max. number of characters */
#define MAXFILE 32768 /* max number of input files */
/*------------ Set defines according to machine's specificities -------------*/
#if 0
#define NO_ENVVAR
#endif
/*--------------------- in case of missing constants ------------------------*/
#ifndef SEEK_SET
#define SEEK_SET 0
#endif
#ifndef SEEK_CUR
#define SEEK_CUR 1
#endif
#ifndef EXIT_SUCCESS
#define EXIT_SUCCESS 0
#endif
#ifndef EXIT_FAILURE
#define EXIT_FAILURE -1
#endif
/*---------------------------- return messages ------------------------------*/
#define RETURN_OK 0
#define RETURN_ERROR (-1)
#define RETURN_FATAL_ERROR (-2)
/*------------------------------- Other Macros ------------------------------*/
#define DEXP(x) exp(2.30258509299*(x)) /* 10^x */
#define QFREAD(ptr, size, afile, fname) \
if (fread(ptr, (size_t)(size), (size_t)1, afile)!=1) \
error(EXIT_FAILURE, "*Error* while reading ", fname)
#define QFWRITE(ptr, size, afile, fname) \
if (fwrite(ptr, (size_t)(size), (size_t)1, afile)!=1) \
error(EXIT_FAILURE, "*Error* while writing ", fname)
#define QFSEEK(afile, offset, pos, fname) \
if (fseek(afile, (offset), pos)) \
error(EXIT_FAILURE,"*Error*: file positioning failed in ", \
fname)
#define QFTELL(pos, afile, fname) \
if ((pos=ftell(afile))==-1) \
error(EXIT_FAILURE,"*Error*: file position unknown in ", \
fname)
#define QCALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)calloc((size_t)(nel),sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QMALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)malloc((size_t)(nel)*sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QREALLOC(ptr, typ, nel) \
{if (!(ptr = (typ *)realloc(ptr, (size_t)(nel)*sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptr " (" #nel " elements) !");;}
#define QMEMCPY(ptrin, ptrout, typ, nel) \
{if (ptrin) \
{if (!(ptrout = (typ *)malloc((size_t)(nel)*sizeof(typ)))) \
error(EXIT_FAILURE, "Not enough memory for ", \
#ptrout " (" #nel " elements) !"); \
memcpy(ptrout, ptrin, (size_t)(nel)*sizeof(typ));};;}
#define QPOPEN(file, cmdline, flag) \
{if (!(file=popen(cmdline, flag))) \
error(EXIT_FAILURE, "*Error*: cannot execute ", cmdline);;}
#define RINT(x) (int)(floor(x+0.5))
#define NPRINTF if (prefs.verbose_type == NORM) fprintf
#define NFPRINTF(w,x) {if (prefs.verbose_type == NORM) \
fprintf(w, "\33[1M> %s\n\33[1A",x);}
#define FPRINTF if (prefs.verbose_type == FULL) fprintf
#define QPRINTF if (prefs.verbose_type != QUIET) fprintf
#define QIPRINTF(w,x) {if (prefs.verbose_type == NORM) \
fprintf(w, "\33[7m%s\33[0m\n", x); \
else if (prefs.verbose_type == LOG) \
fprintf(w, "%s\n", x);}
#define QBPRINTF(w,x) {if (prefs.verbose_type == NORM) \
fprintf(w, "\33[01;31m%s\33[0m\n", x); \
else if (prefs.verbose_type == LOG) \
fprintf(w, "%s\n", x);}
src/levmar/Axb.c 0 → 100644
View file @ ece261cd
/////////////////////////////////////////////////////////////////////////////////
//
// Solution of linear systems involved in the Levenberg - Marquardt
// minimization algorithm
// Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr)
// Institute of Computer Science, Foundation for Research & Technology - Hellas
// Heraklion, Crete, Greece.
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
/////////////////////////////////////////////////////////////////////////////////
/********************************************************************************
* LAPACK-based implementations for various linear system solvers. The same core
* code is used with appropriate #defines to derive single and double precision
* solver versions, see also Axb_core.c
********************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lm.h"
#include "misc.h"
#if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC)
#error At least one of LM_DBL_PREC, LM_SNGL_PREC should be defined!
#endif
#ifdef LM_DBL_PREC
/* double precision definitions */
#define LM_REAL double
#define LM_PREFIX d
#define LM_CNST(x) (x)
#ifndef HAVE_LAPACK
#include <float.h>
#define LM_REAL_EPSILON DBL_EPSILON
#endif
#include "Axb_core.c"
#undef LM_REAL
#undef LM_PREFIX
#undef LM_CNST
#undef LM_REAL_EPSILON
#endif /* LM_DBL_PREC */
#ifdef LM_SNGL_PREC
/* single precision (float) definitions */
#define LM_REAL float
#define LM_PREFIX s
#define __SUBCNST(x) x##F
#define LM_CNST(x) __SUBCNST(x) // force substitution
#ifndef HAVE_LAPACK
#define LM_REAL_EPSILON FLT_EPSILON
#endif
#include "Axb_core.c"
#undef LM_REAL
#undef LM_PREFIX
#undef __SUBCNST
#undef LM_CNST
#undef LM_REAL_EPSILON
#endif /* LM_SNGL_PREC */
src/levmar/Axb_core.c 0 → 100644
View file @ ece261cd
/////////////////////////////////////////////////////////////////////////////////
//
// Solution of linear systems involved in the Levenberg - Marquardt
// minimization algorithm
// Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr)
// Institute of Computer Science, Foundation for Research & Technology - Hellas
// Heraklion, Crete, Greece.
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program 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 General Public License for more details.
//
/////////////////////////////////////////////////////////////////////////////////
/* Solvers for the linear systems Ax=b. Solvers should NOT modify their A & B arguments! */
#ifndef LM_REAL // not included by Axb.c
#error This file should not be compiled directly!
#endif
#ifdef LINSOLVERS_RETAIN_MEMORY
#define __STATIC__ static
#else
#define __STATIC__ // empty
#endif /* LINSOLVERS_RETAIN_MEMORY */
#ifdef HAVE_LAPACK
/* prototypes of LAPACK routines */
#define GEQRF LM_ADD_PREFIX(geqrf_)
#define ORGQR LM_ADD_PREFIX(orgqr_)
#define TRTRS LM_ADD_PREFIX(trtrs_)
#define POTF2 LM_ADD_PREFIX(potf2_)
#define POTRF LM_ADD_PREFIX(potrf_)
#define GETRF LM_ADD_PREFIX(getrf_)
#define GETRS LM_ADD_PREFIX(getrs_)
#define GESVD LM_ADD_PREFIX(gesvd_)
#define GESDD LM_ADD_PREFIX(gesdd_)
/* QR decomposition */
extern int GEQRF(int *m, int *n, LM_REAL *a, int *lda, LM_REAL *tau, LM_REAL *work, int *lwork, int *info);
extern int ORGQR(int *m, int *n, int *k, LM_REAL *a, int *lda, LM_REAL *tau, LM_REAL *work, int *lwork, int *info);
/* solution of triangular systems */
extern int TRTRS(char *uplo, char *trans, char *diag, int *n, int *nrhs, LM_REAL *a, int *lda, LM_REAL *b, int *ldb, int *info);
/* cholesky decomposition */
extern int POTF2(char *uplo, int *n, LM_REAL *a, int *lda, int *info);
extern int POTRF(char *uplo, int *n, LM_REAL *a, int *lda, int *info); /* block version of dpotf2 */
/* LU decomposition and systems solution */
extern int GETRF(int *m, int *n, LM_REAL *a, int *lda, int *ipiv, int *info);
extern int GETRS(char *trans, int *n, int *nrhs, LM_REAL *a, int *lda, int *ipiv, LM_REAL *b, int *ldb, int *info);
/* Singular Value Decomposition (SVD) */
extern int GESVD(char *jobu, char *jobvt, int *m, int *n, LM_REAL *a, int *lda, LM_REAL *s, LM_REAL *u, int *ldu,
LM_REAL *vt, int *ldvt, LM_REAL *work, int *lwork, int *info);
/* lapack 3.0 new SVD routine, faster than xgesvd().
* In case that your version of LAPACK does not include them, use the above two older routines
*/
extern int GESDD(char *jobz, int *m, int *n, LM_REAL *a, int *lda, LM_REAL *s, LM_REAL *u, int *ldu, LM_REAL *vt, int *ldvt,
LM_REAL *work, int *lwork, int *iwork, int *info);
/* precision-specific definitions */
#define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR)
#define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS)
#define AX_EQ_B_CHOL LM_ADD_PREFIX(Ax_eq_b_Chol)
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU)
#define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD)
/*
* This function returns the solution of Ax = b
*
* The function is based on QR decomposition with explicit computation of Q:
* If A=Q R with Q orthogonal and R upper triangular, the linear system becomes
* Q R x = b or R x = Q^T b.
* The last equation can be solved directly.
*
* A is mxm, b is mx1
*
* The function returns 0 in case of error, 1 if successfull
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_QR(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
LM_REAL *a, *qtb, *tau, *r, *work;
int a_sz, qtb_sz, tau_sz, r_sz, tot_sz;
register int i, j;
int info, worksz, nrhs=1;
register LM_REAL sum;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf_sz=0;
buf=NULL;
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
/* calculate required memory size */
a_sz=m*m;
qtb_sz=m;
tau_sz=m;
r_sz=m*m; /* only the upper triangular part really needed */
worksz=3*m; /* this is probably too much */
tot_sz=a_sz + qtb_sz + tau_sz + r_sz + worksz;
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_QR) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_QR) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
a=buf;
qtb=a+a_sz;
tau=qtb+qtb_sz;
r=tau+tau_sz;
work=r+r_sz;
/* store A (column major!) into a */
for(i=0; i<m; i++)
for(j=0; j<m; j++)
a[i+j*m]=A[i*m+j];
/* QR decomposition of A */
GEQRF((int *)&m, (int *)&m, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GEQRF) " in ", AX_EQ_B_QR) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("Unknown LAPACK error %d for ", GEQRF) " in ", AX_EQ_B_QR) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* R is stored in the upper triangular part of a; copy it in r so that ORGQR() below won't destroy it */
for(i=0; i<r_sz; i++)
r[i]=a[i];
/* compute Q using the elementary reflectors computed by the above decomposition */
ORGQR((int *)&m, (int *)&m, (int *)&m, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", ORGQR) " in ", AX_EQ_B_QR) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT("Unknown LAPACK error (%d) in ", AX_EQ_B_QR) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* Q is now in a; compute Q^T b in qtb */
for(i=0; i<m; i++){
for(j=0, sum=0.0; j<m; j++)
sum+=a[i*m+j]*B[j];
qtb[i]=sum;
}
/* solve the linear system R x = Q^t b */
TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, r, (int *)&m, qtb, (int *)&m, &info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_QR) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_QR) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* copy the result in x */
for(i=0; i<m; i++)
x[i]=qtb[i];
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 1;
}
/*
* This function returns the solution of min_x ||Ax - b||
*
* || . || is the second order (i.e. L2) norm. This is a least squares technique that
* is based on QR decomposition:
* If A=Q R with Q orthogonal and R upper triangular, the normal equations become
* (A^T A) x = A^T b or (R^T Q^T Q R) x = A^T b or (R^T R) x = A^T b.
* This amounts to solving R^T y = A^T b for y and then R x = y for x
* Note that Q does not need to be explicitly computed
*
* A is mxn, b is mx1
*
* The function returns 0 in case of error, 1 if successfull
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_QRLS(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m, int n)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
LM_REAL *a, *atb, *tau, *r, *work;
int a_sz, atb_sz, tau_sz, r_sz, tot_sz;
register int i, j;
int info, worksz, nrhs=1;
register LM_REAL sum;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
if(m<n){
fprintf(stderr, RCAT("Normal equations require that the number of rows is greater than number of columns in ", AX_EQ_B_QRLS) "() [%d x %d]! -- try transposing\n", m, n);
exit(1);
}
/* calculate required memory size */
a_sz=m*n;
atb_sz=n;
tau_sz=n;
r_sz=n*n;
worksz=3*n; /* this is probably too much */
tot_sz=a_sz + atb_sz + tau_sz + r_sz + worksz;
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_QRLS) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_QRLS) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
a=buf;
atb=a+a_sz;
tau=atb+atb_sz;
r=tau+tau_sz;
work=r+r_sz;
/* store A (column major!) into a */
for(i=0; i<m; i++)
for(j=0; j<n; j++)
a[i+j*m]=A[i*n+j];
/* compute A^T b in atb */
for(i=0; i<n; i++){
for(j=0, sum=0.0; j<m; j++)
sum+=A[j*n+i]*B[j];
atb[i]=sum;
}
/* QR decomposition of A */
GEQRF((int *)&m, (int *)&n, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GEQRF) " in ", AX_EQ_B_QRLS) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("Unknown LAPACK error %d for ", GEQRF) " in ", AX_EQ_B_QRLS) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* R is stored in the upper triangular part of a. Note that a is mxn while r nxn */
for(j=0; j<n; j++){
for(i=0; i<=j; i++)
r[i+j*n]=a[i+j*m];
/* lower part is zero */
for(i=j+1; i<n; i++)
r[i+j*n]=0.0;
}
/* solve the linear system R^T y = A^t b */
TRTRS("U", "T", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, atb, (int *)&n, &info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_QRLS) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_QRLS) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* solve the linear system R x = y */
TRTRS("U", "N", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, atb, (int *)&n, &info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_QRLS) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_QRLS) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* copy the result in x */
for(i=0; i<n; i++)
x[i]=atb[i];
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 1;
}
/*
* This function returns the solution of Ax=b
*
* The function assumes that A is symmetric & postive definite and employs
* the Cholesky decomposition:
* If A=U^T U with U upper triangular, the system to be solved becomes
* (U^T U) x = b
* This amount to solving U^T y = b for y and then U x = y for x
*
* A is mxm, b is mx1
*
* The function returns 0 in case of error, 1 if successfull
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_CHOL(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
LM_REAL *a, *b;
int a_sz, b_sz, tot_sz;
register int i, j;
int info, nrhs=1;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
/* calculate required memory size */
a_sz=m*m;
b_sz=m;
tot_sz=a_sz + b_sz;
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_CHOL) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_CHOL) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
a=buf;
b=a+a_sz;
/* store A (column major!) into a anb B into b */
for(i=0; i<m; i++){
for(j=0; j<m; j++)
a[i+j*m]=A[i*m+j];
b[i]=B[i];
}
/* Cholesky decomposition of A */
POTF2("U", (int *)&m, a, (int *)&m, (int *)&info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", POTF2) " in ", AX_EQ_B_CHOL) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("LAPACK error: the leading minor of order %d is not positive definite,\nthe factorization could not be completed for ", POTF2) " in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* solve the linear system U^T y = b */
TRTRS("U", "T", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (int *)&m, &info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_CHOL) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* solve the linear system U x = y */
TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, b, (int *)&m, &info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) "in ", AX_EQ_B_CHOL) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* copy the result in x */
for(i=0; i<m; i++)
x[i]=b[i];
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 1;
}
/*
* This function returns the solution of Ax = b
*
* The function employs LU decomposition:
* If A=L U with L lower and U upper triangular, then the original system
* amounts to solving
* L y = b, U x = y
*
* A is mxm, b is mx1
*
* The function returns 0 in case of error,
* 1 if successfull
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_LU(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
int a_sz, ipiv_sz, b_sz, work_sz, tot_sz;
register int i, j;
int info, *ipiv, nrhs=1;
LM_REAL *a, *b, *work;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
/* calculate required memory size */
ipiv_sz=m;
a_sz=m*m;
b_sz=m;
work_sz=100*m; /* this is probably too much */
tot_sz=ipiv_sz + a_sz + b_sz + work_sz; // ipiv_sz counted as LM_REAL here, no harm is done though
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
ipiv=(int *)buf;
a=(LM_REAL *)(ipiv + ipiv_sz);
b=a+a_sz;
work=b+b_sz;
/* store A (column major!) into a and B into b */
for(i=0; i<m; i++){
for(j=0; j<m; j++)
a[i+j*m]=A[i*m+j];
b[i]=B[i];
}
/* LU decomposition for A */
GETRF((int *)&m, (int *)&m, a, (int *)&m, ipiv, (int *)&info);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("argument %d of ", GETRF) " illegal in ", AX_EQ_B_LU) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* solve the system with the computed LU */
GETRS("N", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, b, (int *)&m, (int *)&info);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("argument %d of ", GETRS) " illegal in ", AX_EQ_B_LU) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("unknown error for ", GETRS) " in ", AX_EQ_B_LU) "()\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* copy the result in x */
for(i=0; i<m; i++){
x[i]=b[i];
}
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 1;
}
/*
* This function returns the solution of Ax = b
*
* The function is based on SVD decomposition:
* If A=U D V^T with U, V orthogonal and D diagonal, the linear system becomes
* (U D V^T) x = b or x=V D^{-1} U^T b
* Note that V D^{-1} U^T is the pseudoinverse A^+
*
* A is mxm, b is mx1.
*
* The function returns 0 in case of error, 1 if successfull
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_SVD(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
static LM_REAL eps=LM_CNST(-1.0);
register int i, j;
LM_REAL *a, *u, *s, *vt, *work;
int a_sz, u_sz, s_sz, vt_sz, tot_sz;
LM_REAL thresh, one_over_denom;
register LM_REAL sum;
int info, rank, worksz, *iwork, iworksz;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
/* calculate required memory size */
worksz=16*m; /* more than needed */
iworksz=8*m;
a_sz=m*m;
u_sz=m*m; s_sz=m; vt_sz=m*m;
tot_sz=iworksz*sizeof(int) + (a_sz + u_sz + s_sz + vt_sz + worksz)*sizeof(LM_REAL);
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_SVD) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_SVD) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
iwork=(int *)buf;
a=(LM_REAL *)(iwork+iworksz);
/* store A (column major!) into a */
for(i=0; i<m; i++)
for(j=0; j<m; j++)
a[i+j*m]=A[i*m+j];
u=a + a_sz;
s=u+u_sz;
vt=s+s_sz;
work=vt+vt_sz;
/* SVD decomposition of A */
GESVD("A", "A", (int *)&m, (int *)&m, a, (int *)&m, s, u, (int *)&m, vt, (int *)&m, work, (int *)&worksz, &info);
//GESDD("A", (int *)&m, (int *)&m, a, (int *)&m, s, u, (int *)&m, vt, (int *)&m, work, (int *)&worksz, iwork, &info);
/* error treatment */
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GESVD), "/" GESDD) " in ", AX_EQ_B_SVD) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT("LAPACK error: dgesdd (dbdsdc)/dgesvd (dbdsqr) failed to converge in ", AX_EQ_B_SVD) "() [info=%d]\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
if(eps<0.0){
LM_REAL aux;
/* compute machine epsilon */
for(eps=LM_CNST(1.0); aux=eps+LM_CNST(1.0), aux-LM_CNST(1.0)>0.0; eps*=LM_CNST(0.5))
;
eps*=LM_CNST(2.0);
}
/* compute the pseudoinverse in a */
for(i=0; i<a_sz; i++) a[i]=0.0; /* initialize to zero */
for(rank=0, thresh=eps*s[0]; rank<m && s[rank]>thresh; rank++){
one_over_denom=LM_CNST(1.0)/s[rank];
for(j=0; j<m; j++)
for(i=0; i<m; i++)
a[i*m+j]+=vt[rank+i*m]*u[j+rank*m]*one_over_denom;
}
/* compute A^+ b in x */
for(i=0; i<m; i++){
for(j=0, sum=0.0; j<m; j++)
sum+=a[i*m+j]*B[j];
x[i]=sum;
}
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 1;
}
/* undefine all. IT MUST REMAIN IN THIS POSITION IN FILE */
#undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS
#undef AX_EQ_B_CHOL
#undef AX_EQ_B_LU
#undef AX_EQ_B_SVD
#undef GEQRF
#undef ORGQR
#undef TRTRS
#undef POTF2
#undef POTRF
#undef GETRF
#undef GETRS
#undef GESVD
#undef GESDD
#else // no LAPACK
/* precision-specific definitions */
/* Added by EB */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include ATLAS_LAPACK_H
#define ATLAS_POTRF LM_CAT_(clapack_, LM_ADD_PREFIX(potrf))
#define ATLAS_POTRS LM_CAT_(clapack_, LM_ADD_PREFIX(potrs))
#define GETRF LM_ADD_PREFIX(getrf_)
#define GETRS LM_ADD_PREFIX(getrs_)
extern int GETRF(int *m, int *n, LM_REAL *a, int *lda, int *ipiv, int *info);
extern int GETRS(char *trans, int *n, int *nrhs, LM_REAL *a, int *lda, int *ipiv, LM_REAL *b, int *ldb, int *info);
/* End added by EB */
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
/*
* This function returns the solution of Ax = b
*
* The function employs LU decomposition followed by forward/back substitution (see
* also the LAPACK-based LU solver above)
*
* A is mxm, b is mx1
*
* The function returns 0 in case of error,
* 1 if successfull
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_LU(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
int a_sz, ipiv_sz, b_sz, work_sz, tot_sz;
register int i, j;
int info, *ipiv;
LM_REAL *a, *b, *work;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
/* calculate required memory size */
ipiv_sz=m;
a_sz=m*m;
b_sz=m;
work_sz=100*m; /* this is probably too much */
tot_sz=ipiv_sz + a_sz + b_sz + work_sz; // ipiv_sz counted as LM_REAL here, no harm is done though
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
ipiv=(int *)buf;
a=(LM_REAL *)(ipiv + ipiv_sz);
b=a+a_sz;
work=b+b_sz;
/* store A (column major!) into a and B into b */
for(i=0; i<m; i++){
for(j=0; j<m; j++)
a[i+j*m]=A[i*m+j];
b[i]=B[i];
}
/* LU decomposition for A */
info = ATLAS_POTRF(CblasRowMajor, CblasUpper, m, a, m);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("argument %d of ", GETRF) " illegal in ", AX_EQ_B_LU) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* solve the system with the computed LU */
info = ATLAS_POTRS(CblasRowMajor, CblasUpper, m, 1, a, m, b, m);
if(info!=0){
if(info<0){
fprintf(stderr, RCAT(RCAT("argument %d of ", GETRS) " illegal in ", AX_EQ_B_LU) "()\n", -info);
exit(1);
}
else{
fprintf(stderr, RCAT(RCAT("unknown error for ", GETRS) " in ", AX_EQ_B_LU) "()\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
}
/* copy the result in x */
for(i=0; i<m; i++){
x[i]=b[i];
}
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 1;
}
/* undefine all. IT MUST REMAIN IN THIS POSITION IN FILE */
#undef AX_EQ_B_LU
/* Added by EB */
#undef GETRF
#undef GETRS
/* End Added by EB */
#endif /* HAVE_LAPACK */
src/levmar/CMakeLists.txt 0 → 100644
View file @ ece261cd
# levmar CMake file; see http://www.cmake.org and
# http://www.insightsoftwareconsortium.org/wiki/index.php/CMake_Tutorial
PROJECT(LEVMAR)
#CMAKE_MINIMUM_REQUIRED(VERSION 1.4)
# compiler flags
ADD_DEFINITIONS(-DLINSOLVERS_RETAIN_MEMORY) # do not free memory between linear solvers calls
#REMOVE_DEFINITIONS(-DLINSOLVERS_RETAIN_MEMORY)
# f2c is sometimes equivalent to libF77 & libI77; in that case, set HAVE_F2C to 0
SET(HAVE_F2C 1 CACHE BOOL "Do we have f2c or F77/I77?" )
# the directory where the lapack/blas/f2c libraries reside
SET(LAPACKBLAS_DIR /usr/lib CACHE PATH "Path to lapack/blas libraries")
# actual names for the lapack/blas/f2c libraries
SET(LAPACK_LIB lapack CACHE STRING "The name of the lapack library")
SET(BLAS_LIB blas CACHE STRING "The name of the blas library")
IF(HAVE_F2C)
SET(F2C_LIB f2c CACHE STRING "The name of the f2c library")
ELSE(HAVE_F2C)
SET(F77_LIB libF77 CACHE STRING "The name of the F77 library")
SET(I77_LIB libI77 CACHE STRING "The name of the I77 library")
ENDIF(HAVE_F2C)
########################## NO CHANGES BEYOND THIS POINT ##########################
#INCLUDE_DIRECTORIES(/usr/include)
LINK_DIRECTORIES(${LAPACKBLAS_DIR})
# levmar library source files
ADD_LIBRARY(levmar STATIC
lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c
lm.h misc.h compiler.h
)
# demo program
ADD_EXECUTABLE(lmdemo lmdemo.c lm.h)
# libraries the demo depends on
IF(HAVE_F2C)
TARGET_LINK_LIBRARIES(lmdemo levmar ${LAPACK_LIB} ${BLAS_LIB} ${F2C_LIB})
ELSE(HAVE_F2C)
TARGET_LINK_LIBRARIES(lmdemo levmar ${LAPACK_LIB} ${BLAS_LIB} ${F77_LIB} ${I77_LIB})
ENDIF(HAVE_F2C)
# make sure that the library is built before the demo
ADD_DEPENDENCIES(lmdemo levmar)
#SUBDIRS(matlab)
#ADD_TEST(levmar_tst lmdemo)
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