C @(#)genxx2.for	17.1.1.1 (ESO-DMD) 01/25/02 17:40:58
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 Correspondence 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
      SUBROUTINE GXMOVIT(A,B,AOFF,BOFF,NOPIX,NOLINS,AXSIZE,BXSIZE,RCUTS)
C 
C  copy stuff from array A to array B
C 
C     .................
C     .................          ..............
C     ....aaaaaaa......          .bbbbbbb......
C     ....aaaaaaa......   ==>    .bbbbbbb......
C     ....aaaaaaa......          .bbbbbbb......
C     .................          ..............
C     .................          ..............
C     .................
C 
      IMPLICIT NONE
C 
      REAL    A(*),B(*),RCUTS(2)
C 
      INTEGER AOFF,BOFF,NOPIX,NOLINS,AXSIZE,BXSIZE
      INTEGER NY,NX,AAOFF,BBOFF,KOFF,LOFF,LINLIM
C 
      AAOFF = AOFF
      BBOFF = BOFF                         !we don't modify original offsets...
C 
      IF ((NOPIX+BOFF-1).GT.BXSIZE) THEN
         LINLIM = BXSIZE - BOFF + 1
      ELSE
         LINLIM = NOPIX
      ENDIF
C 
      DO 2000 NY=1,NOLINS
         KOFF = AAOFF
         LOFF = BBOFF
         DO 1000 NX=1,LINLIM
            B(LOFF) = A(KOFF)
            IF (B(LOFF).LT.RCUTS(1)) RCUTS(1) = B(LOFF)
            IF (B(LOFF).GT.RCUTS(2)) RCUTS(2) = B(LOFF)
            LOFF = LOFF + 1
            KOFF = KOFF + 1
1000     CONTINUE
         AAOFF = AAOFF + AXSIZE
         BBOFF = BBOFF + BXSIZE
2000  CONTINUE
C      
      RETURN
      END

      SUBROUTINE GXDOIT(A,B,C,THR,SIZE,TRUSIZ)
C
      IMPLICIT NONE
C
      INTEGER    SIZE,TRUSIZ
      INTEGER    N,NOUT
C
      REAL       A(SIZE),B(SIZE),C(*),THR
C
      NOUT = 0
      DO 500, N=1,SIZE
         IF (B(N).GE.THR) THEN
            NOUT = NOUT + 1
            C(NOUT) = A(N)
         ENDIF
500   CONTINUE
C
      TRUSIZ = NOUT
      RETURN
      END

      SUBROUTINE GXIZAMAP(NPIXA,NPIXB,STARTB,STEPB,A,B,C,FLAG)
C
      IMPLICIT NONE
C
      REAL         A(*),B(*),C(*)
      REAL         RVAL,ENDB,DIFF,DIFF1,DIFF2
C
      DOUBLE PRECISION STARTB,STEPB
C
      INTEGER      NPIXA(*),NPIXB(*),FLAG(*)
      INTEGER      N,KIDX,TOTAL,IOFF,NN,MCLOSE
C
      TOTAL = NPIXA(1) * NPIXA(2) * NPIXA(3)
      IF (FLAG(2).EQ.1) GOTO 3000
      IF (FLAG(2).EQ.2) GOTO 5000
C
C -------------------------------------------------------C
C
C  here for equidistant map space
C
C -------------------------------------------------------C
C
      ENDB = STARTB + ((NPIXB(1) - 1) * STEPB)
      IF (FLAG(1).EQ.1) GOTO 1000
C
      DO 400 N=1,TOTAL
         RVAL = A(N)
C
         IF (RVAL.LE.STARTB) THEN
            KIDX = 1
         ELSE IF (RVAL.GE.ENDB) THEN
            KIDX = NPIXB(1)
         ELSE
            KIDX = NINT((RVAL-STARTB)/STEPB) + 1
         ENDIF
C
         C(N) = B(KIDX)
400   CONTINUE
C
      RETURN
C
C  here we leave all pixels with intensities outside the map coordinate space
C  as they are
C
1000  DO 1400 N=1,TOTAL
         RVAL = A(N)
C
         IF ( (RVAL.LT.STARTB) .OR. (RVAL.GT.ENDB) ) THEN
            C(N) = A(N)
         ELSE
            KIDX = NINT((RVAL-STARTB)/STEPB) + 1
            C(N) = B(KIDX)
         ENDIF
1400  CONTINUE
C
      RETURN
C
C -------------------------------------------------------C
C
C  here we have the coordinates stored in the 1. line of the map frame,
C  intensities follow in the 2. line
C
C -------------------------------------------------------C
C
3000  IOFF = NPIXB(1)
      STARTB = B(1)
      ENDB = B(IOFF)
      IF (FLAG(1).EQ.1) GOTO 4000
C
      DO 3500 N=1,TOTAL
         RVAL = A(N)
C
         IF (RVAL.LE.STARTB) THEN
            KIDX = 1 + IOFF
         ELSE IF (RVAL.GE.ENDB) THEN
            KIDX = IOFF + IOFF
         ELSE
            MCLOSE = 1
            DIFF1 = ABS(RVAL-STARTB)
C
            DO 3200 NN=2,IOFF                  !find closest coordinate
               DIFF = ABS(RVAL-B(NN))
               IF (DIFF.LT.DIFF1) THEN
                  DIFF1 = DIFF
                  MCLOSE = NN
               ELSE
                  KIDX = MCLOSE + IOFF
                  GOTO 3333
               ENDIF
3200        CONTINUE
         ENDIF
C
3333     C(N) = B(KIDX)
3500  CONTINUE
C
      RETURN
C
C  here we leave all pixels with intensities outside the map coordinate space
C  as they are
C
4000  DO 4500 N=1,TOTAL
         RVAL = A(N)
C
         IF ( (RVAL.LT.STARTB) .OR. (RVAL.GT.ENDB) ) THEN
            C(N) = A(N)
C
         ELSE
            MCLOSE = 1
            DIFF1 = ABS(RVAL-STARTB)
C
            DO 4200 NN=2,IOFF                  !find closest coordinate
               DIFF = ABS(RVAL-B(NN))
               IF (DIFF.LT.DIFF1) THEN
                  DIFF1 = DIFF
                  MCLOSE = NN
               ELSE
                  KIDX = MCLOSE + IOFF
                  GOTO 4333
               ENDIF
4200        CONTINUE
4333        C(N) = B(KIDX)
         ENDIF
C
4500  CONTINUE
C
      RETURN
C
C -------------------------------------------------------C
C
C  here we have the intervals stored in the 1. line of the map frame,
C  intensities follow in the 2. line (for start,end - in between interpolation ...)
C
C -------------------------------------------------------C
C
5000  IOFF = NPIXB(1)
      STARTB = B(1)
      ENDB = B(IOFF)
      IF (FLAG(1).EQ.1) GOTO 6000
C
      DO 5555 N=1,TOTAL
         RVAL = A(N)
C
         IF (RVAL.LT.STARTB) THEN
            RVAL = STARTB
         ELSE IF (RVAL.GT.ENDB) THEN
            RVAL = ENDB
         ENDIF
C
         DO 5200 NN=1,IOFF,2                 !find enclosing interval
            IF (RVAL.LE.B(NN+1)) THEN
                KIDX = NN
                GOTO 5333
            ELSE IF (RVAL.LT.B(NN+2)) THEN
                C(N) = A(N)
                GOTO 5555                    !no interval there ...!
            ENDIF
5200     CONTINUE
C
5333     DIFF = B(KIDX+1) - B(KIDX)
         DIFF1 = (RVAL - B(KIDX)) / DIFF
         DIFF2 = (B(KIDX+1) - RVAL) / DIFF
         KIDX = KIDX + IOFF                  !point to data line
         C(N) = ( B(KIDX+1) * DIFF1 ) + ( B(KIDX) * DIFF2 )
5555  CONTINUE
C
      RETURN
C
C  here we leave all pixels with intensities outside the map coordinate space
C  as they are
C
6000  DO 6555 N=1,TOTAL
         RVAL = A(N)
C
         IF (RVAL.LT.STARTB) THEN
            C(N) = A(N)
            GOTO 6555
         ELSE IF (RVAL.GT.ENDB) THEN
            C(N) = A(N)
            GOTO 6555
         ENDIF
C
         DO 6200 NN=1,IOFF,2                  !find enclosing interval
            IF (RVAL.LE.B(NN+1)) THEN
               KIDX = NN
               GOTO 6333
            ELSE IF (RVAL.LT.B(NN+2)) THEN
               C(N) = A(N)
               GOTO 6555                      !no interval there ...!
            ENDIF
6200     CONTINUE
C
6333     DIFF = B(KIDX+1) - B(KIDX)
         DIFF1 = (RVAL - B(KIDX)) / DIFF
         DIFF2 = (B(KIDX+1) - RVAL) / DIFF
         KIDX = KIDX + IOFF                   !point to data line
         C(N) = ( B(KIDX+1) * DIFF1 ) + ( B(KIDX) * DIFF2 )
6555  CONTINUE
C
      RETURN
      END

      SUBROUTINE GXMATRX(A,NPIX,C)
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  straight forward transposition of an image (seen as a matrix)
C  this however implies, that the first element is in the upper left corner
C  since matrices are usually written like that ...
C
C -------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      NPIX(*)
      INTEGER      NX,NY
      INTEGER      XDIM,YDIM,OFFA,OFFC,WORK
C
      REAL         A(*),C(*)
C
C  init
      XDIM = NPIX(1)
      YDIM = NPIX(2)
      OFFC = 0
      WORK = XDIM * (YDIM-1)
C
C  loop through rows in output matrix
      DO 1000 NX=XDIM,1,-1
         OFFA = WORK
C
         DO 400 NY=1,YDIM                          !loop till middle of row
            C(NY+OFFC) = A(NX+OFFA)
            OFFA = OFFA - XDIM
400      CONTINUE
C
         OFFC = OFFC + YDIM
1000  CONTINUE
      RETURN
C
      END