C @(#)gridtab.for	17.1.1.1 (ES0-DMD) 01/25/02 17:19:20
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT:    Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION: GRIDTAB
C.PURPOSE:      Produce 2-D coordinate frame with rectangular grid subjected
C               to "S-type" distortion (output to TABLE)
C.LANGUAGE:     F77+ESOext
C.AUTHOR:       D.BAADE
C.KEYWORDS:     Ingeneering frames
C.ALGORITHM:    trivial
C.PECULARITIES: The size of the coordinate frame is fixed to 510 by 510 pixels
C               as is the separation of the lines of the undistorted grid to 
C               60 pixels
C.INPUT/OUTPUT: The following key words are used:
C               ALPHA/R/1/1  angle ALPHA0 of the distortion in the following
C                            parametrizion  ALPHA = ALPHA0 * R**2, with R 
C                            being the distance from the point about which
C                            the rotation takes place. ALPHA0 to be given 
C                            in degrees/pixel.
C.VERSION:      850322  Creation                D. Baade, ST-ECF
C.VERSION:      850820  ???                     K. Banse
C.VERSION:      871123  ESO-FORTRAN Conversion, Rein H. Warmels
C --------------------------------------------------------------------
      PROGRAM GRDTAB
C
      IMPLICIT  NONE

      INTEGER   MADRID
      INTEGER   IAV,I,J
      INTEGER   NSEP
      REAL      X,Y,XD,YD,ALPHA,ALPHA0,RSQUAR,XORIG,YORIG,DATA(4)
      INTEGER   ISTAT,COLUMN(4),IROW
      INTEGER   TIDO,NCOL,NROW
      INTEGER   KUN,KNUL
      CHARACTER TBOUT*64,FORM*4
C
      INCLUDE   'MID_INCLUDE:ST_DEF.INC'
      COMMON    /VMR/MADRID(1)
      INCLUDE   'MID_INCLUDE:ST_DAT.INC'
C
      DATA      NSEP/60/
      DATA      COLUMN/1,2,3,4/    
      DATA      XORIG/100./
      DATA      YORIG/255./
      DATA      TBOUT/'GRID'/
      DATA      FORM/'F8.3'/
      DATA      NCOL/4/
      DATA      NROW/81/
C
C *** begin code
      CALL STSPRO('GRIDTAB')                        ! get into midas environment
C
      CALL STKRDR('ALPHA',1,1,IAV,ALPHA0,KUN,KNUL,ISTAT)  ! get distortion angle
      ALPHA0 = ALPHA0*3.14159/180.                        ! convert to radians
      CALL TBTINI(TBOUT,0,F_O_MODE,NCOL,NROW,TIDO,ISTAT)  ! output table (new)
C
C *** define columns
      CALL TBCINI(TIDO,D_R4_FORMAT,1,FORM,'PIXELS','X_RECTANGULAR',
     2            COLUMN(1),ISTAT)
      CALL TBCINI(TIDO,D_R4_FORMAT,1,FORM,'PIXELS','Y_RECTANGULAR',
     2            COLUMN(2),ISTAT)
      CALL TBCINI(TIDO,D_R4_FORMAT,1,FORM,'PIXELS','X_ACTUAL',
     2            COLUMN(3),ISTAT)
      CALL TBCINI(TIDO,D_R4_FORMAT,1,FORM,'PIXELS','Y_ACTUAL',
     2            COLUMN(4),ISTAT)
C
C ... The following command sequence creates 'reseau marks' in a
C ... regularly spaced rectangular grid. using the simple
C ... parametrization:
C ...                    ALPHA = ALPHA0 * R**2
C ... With R being the distance from the origin about which the
C ... rotation takes place, this first grid is distorted into a
C ... second one. Output of x,y_rectangular is to columns #1 and #2
C ... of table 'GRID.TBL'. Columns #3 and #4 receive x,y_actual.
      IROW   = 0
      DO 20 I = 1,9
         X = (I-1)*NSEP + 15. - XORIG
         DO 10 J = 1,9
            IROW    = IROW + 1
            Y       = (J-1)*NSEP + 15. - YORIG
            RSQUAR  = X**2 + Y**2
            ALPHA   = RSQUAR*ALPHA0
            XD      = X*COS(ALPHA) - Y*SIN(ALPHA) + XORIG
            YD      = X*SIN(ALPHA) + Y*COS(ALPHA) + YORIG
            DATA(1) = X + XORIG
            DATA(2) = Y + YORIG
            DATA(3) = XD
            DATA(4) = YD
            CALL TBRWRR(TIDO,IROW,4,COLUMN,DATA,ISTAT)
   10     CONTINUE
   20 CONTINUE
C
C *** release files and update keywords
      CALL TBTCLO(TIDO,ISTAT)
      CALL STSEPI
      END