*=======================================================================
*
*   WCSLIB 4.3 - an implementation of the FITS WCS standard.
*   Copyright (C) 1995-2007, Mark Calabretta
*
*   This file is part of WCSLIB.
*
*   WCSLIB is free software: you can redistribute it and/or modify
*   it under the terms of the GNU Lesser General Public License as
*   published by the Free Software Foundation, either version 3 of
*   the License, or (at your option) any later version.
*
*   WCSLIB is distributed in the hope that it will be useful, but
*   WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*   GNU Lesser General Public License for more details.
*
*   You should have received a copy of the GNU Lesser General Public
*   License along with WCSLIB.  If not, see http://www.gnu.org/licenses.
*
*   Correspondence concerning WCSLIB may be directed to:
*      Internet email: mcalabre@atnf.csiro.au
*      Postal address: Dr. Mark Calabretta
*                      Australia Telescope National Facility, CSIRO
*                      PO Box 76
*                      Epping NSW 1710
*                      AUSTRALIA
*
*   Author: Mark Calabretta, Australia Telescope National Facility
*   http://www.atnf.csiro.au/~mcalabre/index.html
*   $Id: twcstab.f,v 4.3 2007/12/27 05:42:52 cal103 Exp $
*=======================================================================

      PROGRAM TWCSTAB
*-----------------------------------------------------------------------
*
*   TWCSTAB tests WCSTAB and also provides sample code for using it in
*   conjunction with WCSPIH and FTWCST.  Although this example and
*   FTWCST are based on the CFITSIO library, WCSTAB itself is completely
*   independent of it.
*
*   We assume that the input file, ../C/wcstab.fits, has already been
*   generated by running the C version of twcstab.
*
*=======================================================================

      LOGICAL   GOTEND
      INTEGER   ARRAY, BLOKSZ, DIMLEN, EXTLEV, EXTNAM, EXTVER, I, J, K,
     :          IERR, IUNIT, NKEYRC, NDIM, NREJECT, NWCS, NWTB, STATUS,
     :          TTYPE, WCSP, WTB
      CHARACTER KEYREC*80, HEADER*28801, INFILE*16

      INCLUDE 'wcs.inc'
      INCLUDE 'wcshdr.inc'
      INCLUDE 'wcsfix.inc'
      INCLUDE 'getwcstab.inc'
      INTEGER STAT(WCSFIX_NWCS)
      INTEGER WCS(WCSLEN)

      DATA INFILE /'../C/wcstab.fits'/
*-----------------------------------------------------------------------
      WRITE (*, 10)
 10   FORMAT ('Testing WCSTAB and associated routines (twcstab.f)',/,
     :        '--------------------------------------------------',/)

*     Open the FITS test file.
      IUNIT = 1
      STATUS = 0
      CALL FTOPEN (IUNIT, INFILE, 0, BLOKSZ, STATUS)
      IF (STATUS.NE.0) THEN
         CALL FTRPRT ('STDERR', STATUS)
         GO TO 999
      END IF

*     Read the primary header; unfortunately the is no FITSIO equivalent
*     of fits_hdr2str() in CFITSIO so do it the long way.
      OPEN (UNIT=1, FILE=INFILE, FORM='FORMATTED', ACCESS='DIRECT',
     :      RECL=80, IOSTAT=IERR)
      IF (IERR.NE.0) THEN
         WRITE (*, 20) IERR, INFILE
 20      FORMAT ('ERROR',I3,' opening ',A)
         GO TO 999
      END IF

*     Read in the FITS header, excluding COMMENT and HISTORY keyrecords.
      K = 1
      NKEYRC = 0
      GOTEND = .FALSE.
      DO 50 J = 0, 100
         DO 40 I = 1, 36
            READ (1, '(A80)', REC=36*J+I, IOSTAT=IERR) KEYREC
            IF (IERR.NE.0) THEN
               WRITE (*, 30) IERR
 30            FORMAT ('ERROR',I3,' reading header.')
               GO TO 999
            END IF

            IF (KEYREC(:8).EQ.'        ') GO TO 40
            IF (KEYREC(:8).EQ.'COMMENT ') GO TO 40
            IF (KEYREC(:8).EQ.'HISTORY ') GO TO 40

            HEADER(K:) = KEYREC
            K = K + 80
            NKEYRC = NKEYRC + 1

            IF (KEYREC(:8).EQ.'END     ') THEN
*              An END keyrecord was read, read the rest of the block.
               GOTEND = .TRUE.
            END IF
 40      CONTINUE

         IF (GOTEND) GO TO 60
 50   CONTINUE

 60   CLOSE (UNIT=1)

*-----------------------------------------------------------------------
* Basic steps required to interpret a FITS WCS header, including -TAB.
*-----------------------------------------------------------------------

*     Parse the primary header of the FITS file.
      STATUS = WCSPIH (HEADER, NKEYRC, WCSHDR_all, 2, NREJECT, NWCS,
     :                 WCSP)
      IF (STATUS.NE.0) THEN
         WRITE (*, 70) STATUS, WCSHDR_ERRMSG(STATUS)
 70      FORMAT ('WCSPIH ERROR',I2,A)
         GO TO 999
      END IF

*     Copy into our WCSPRM struct.
      IERR = WCSVCOPY (WCSP, 0, WCS)

*     Read coordinate arrays from the binary table extension.
      STATUS = WCSGET (WCS, WCS_NWTB, NWTB)
      STATUS = WCSGET (WCS, WCS_WTB,  WTB)
      STATUS = FTWCST (IUNIT, NWTB, WTB, STATUS)
      IF (STATUS.NE.0) THEN
         CALL FTRPRT ('STDERR', STATUS)
         GO TO 999
      END IF

*     Fix non-standard WCS keyvalues.
      STATUS = WCSFIX (7, 0, WCS, STAT)
      IF (STATUS.NE.0) THEN
         WRITE (*, 80) (STAT(I), I=1,WCSFIX_NWCS)
 80      FORMAT ('WCSFIX ERROR, status returns: ',10(I2,:,','))
         RETURN
      END IF

*-----------------------------------------------------------------------
* The wcsprm struct is now ready for use.
*-----------------------------------------------------------------------

*     Finished with the FITS file.
      CALL FTCLOS (IUNIT, STATUS)

*     Initialize the wcsprm struct, taking memory allocated by FTWCST.
      IF (STATUS.EQ.0) STATUS = WCSSET (WCS)
      IF (STATUS.NE.0) THEN
         WRITE (*, 80) STATUS, WCS_ERRMSG(STATUS)
 90      FORMAT ('WCSSET ERROR',I2,A)
         GO TO 999
      END IF

*     Do something with it.
      STATUS = WCSPRT (WCS)

*     Clean up.
      STATUS = WCSFREE (WCS)
      STATUS = WCSVFREE (NWCS, WCSP)

999   RETURN
      END