C===========================================================================
C Copyright (C) 1995,2003 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 Massachusetts 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
      PROGRAM FLIP
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program FLIP			version 1.00	881028
C  K. Banse                  	ESO - Garching
C
C.KEYWORDS
C  bulk data frame, flipping
C
C.PURPOSE
C  flip an image in x and/or y dimension
C
C.ALGORITHM
C  straight forward, data is changed in place
C
C.INPUT/OUTPUT
C  the following keywords are used:
C
C  IN_A/C/1/60				name of input frame
C  ACTION/C/1/2				action
C					= 'X ', flip in x
C					  'Y ', flip in y
C					  'XY', flip in x and y
C					  'SH', shift the image
C					  'SC', scale the image
C					  'T ', apply thinning algorithm
C					  'C ', for correlation
C
C  if ACTION = SH, then
C  OUT_A/C/1/60				name of result fram
C  INPUTI/I/1/2				xshift,yshift in [1,NPIX-1]
C
C  if ACTION = SC, then
C  OUT_A/C/1/60				name of result fram
C  INPUTI/I/1/2				x, y scaling values,
C					> 1  - enlarge image by replication
C					< -1 - decrease image size 
C					by averaging or omission
C  ACTION(3:3)				A or N for averaging or omission
C  if ACTION = T, then
C  OUT_A/C/1/60				name of result fram
C
C  if ACTION = C, then
C  IN_A/C/1/60			      	template spectrum
C  IN_B/C/1/1/60			spectrum
C  OUT_A/C/1/60				correlation table
C  INPUTI/I/1/1				max shift in pixels
C  ACTION(2:2)				Y or N for normalizing the pixels
C					of the output frame
C
C.VERSIONS
C 030731		last modif
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      IAV,NAXIS,STAT,KFLAG,FORMA,IOMOD
      INTEGER*8    PNTR,PNTRB,PNTRC,WPNTR
      INTEGER      NPIX(2),IMNO,NPIXC(2),IMNOC,SCALE(3)
      INTEGER      NAXISB,NPIXB(2),IMNOB,ISHI
      INTEGER      SIZE,XSHIFT,YSHIFT
      INTEGER      UNI(1),NULO,MADRID(1)
C
      DOUBLE PRECISION START(2),STEP(2),STARTB(2),STEPB(2)
      DOUBLE PRECISION END(2),DIFSTP,EPS
C
      REAL         CUTS(2),MINVAL,MAXVAL,LHCUTS(4)
C
      CHARACTER    CUNIT*48,IDENT*72
      CHARACTER    FRAME*84,FRAMEB*84,FRAMEC*84,ACTION*4
C
      COMMON   /VMR/  MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      DATA    CUTS   /0.,0./
      DATA    FRAME  /' '/, FRAMEB /' '/, FRAMEC /' '/
      DATA    IDENT  /' '/,  CUNIT  /' '/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  get into MIDAS
      CALL STSPRO('FLIP')
C
C  get name of input image + action flags
      CALL STKRDC('IN_A',1,1,80,IAV,FRAME,UNI,NULO,STAT)
      CALL STKRDC('ACTION',1,1,3,IAV,ACTION,UNI,NULO,STAT)
      CALL UPCAS(ACTION,ACTION)
C
C  setup action flag
      IF ((ACTION(1:1).EQ.'X') .OR. (ACTION(1:1).EQ.'Y')) THEN 
         KFLAG = 1
         FORMA = D_R4_FORMAT
         IOMOD = F_IO_MODE
      ELSE IF (ACTION(1:2).EQ.'SH') THEN
         KFLAG = 2
         FORMA = D_R4_FORMAT
         IOMOD = F_I_MODE
      ELSE IF (ACTION(1:2).EQ.'SC') THEN
         KFLAG = 3
         FORMA = D_R4_FORMAT
         IOMOD = F_I_MODE
      ELSE IF (ACTION(1:1).EQ.'T') THEN
         KFLAG = 4
         FORMA = D_I4_FORMAT
         IOMOD = F_I_MODE
      ELSE IF (ACTION(1:1).EQ.'C') THEN
         KFLAG = 5
         FORMA = D_R4_FORMAT
         IOMOD = F_I_MODE
      ELSE
         CALL STETER(77,'Invalid action ...')
      ENDIF
C 
C  map input frame 
      CALL STIGET(FRAME,FORMA,IOMOD,F_IMA_TYPE,
     +            2,NAXIS,NPIX,START,STEP,IDENT,
     +            CUNIT,PNTR,IMNO,STAT)
      CALL STDRDR(IMNO,'LHCUTS',1,4,IAV,LHCUTS,UNI,NULO,STAT) 
C
C  --------------------------
C  branch according to action
C  --------------------------
C
      IF (KFLAG.EQ.1) THEN                                    !flipping ...
         CALL STFXMP(NPIX(1),D_R4_FORMAT,WPNTR,STAT)
         CALL FLIPPI(MADRID(PNTR),MADRID(WPNTR),NPIX,ACTION)
C 
C  update START + STEP descriptors accordingly
         IF (ACTION(1:1).EQ.'Y') THEN
            END(2) = START(2) + (NPIX(2)-1)*STEP(2)
            START(2) = END(2)
            STEP(2) = -STEP(2)
         ELSE IF (ACTION(2:2).EQ.'Y') THEN
            END(1) = START(1) + (NPIX(1)-1)*STEP(1)
            END(2) = START(2) + (NPIX(2)-1)*STEP(2)
            START(1) = END(1)
            START(2) = END(2)
            STEP(1) = -STEP(1)
            STEP(2) = -STEP(2)
         ELSE
            END(1) = START(1) + (NPIX(1)-1)*STEP(1)
            START(1) = END(1)
            STEP(1) = -STEP(1)
         ENDIF
         CALL STDWRD(IMNO,'START',START,1,2,UNI,STAT) 
         CALL STDWRD(IMNO,'STEP',STEP,1,2,UNI,STAT) 
         IMNOC = IMNO                                    !needed for DSCUPT ...
         GOTO 9000
      ENDIF
C
C  for shifting, thinning, scaling and correlation get output frame
C
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT)
C
C
      IF (KFLAG.EQ.2) THEN                                    !shifting ...
         CALL STKRDI('INPUTI',1,2,IAV,NPIXC,UNI,NULO,STAT)
         XSHIFT = NPIXC(1)
         YSHIFT = NPIXC(2)
500      IF (XSHIFT.GT.NPIX(1)) THEN
            XSHIFT = XSHIFT - NPIX(1)
            GOTO 500
         ELSE IF (XSHIFT.LT.0) THEN
            XSHIFT = NPIX(1) + XSHIFT            !remember, XSHIFT is < 0
            GOTO 500
         ENDIF
600      IF (YSHIFT.GT.NPIX(2)) THEN
            YSHIFT = YSHIFT - NPIX(2)
            GOTO 600
         ELSE IF (YSHIFT.LT.0) THEN
            YSHIFT = NPIX(2) + YSHIFT            !remember, XSHIFT is < 0
            GOTO 600
         ENDIF
         IF ((XSHIFT.EQ.0).AND.(YSHIFT.EQ.0)) GOTO 9000
         CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXIS,NPIX,START,STEP,IDENT,
     +               CUNIT,PNTRC,IMNOC,STAT)
C 
C
C  if shift in both directions, get virtual memory for working buffer
         IF ((XSHIFT.NE.0).AND.(YSHIFT.NE.0)) THEN
            SIZE = NPIX(1) * NPIX(2)
            CALL STFXMP(SIZE,D_R4_FORMAT,WPNTR,STAT)
         ELSE
            WPNTR = PNTR
         ENDIF
C 
         CALL SHIFTI(MADRID(PNTR),MADRID(WPNTR),MADRID(PNTRC),
     +               NPIX,XSHIFT,YSHIFT)
         CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,1,4,UNI,STAT) 
C 
      ELSE IF (KFLAG.EQ.3) THEN                               !scaling ...
         CALL STKRDI('INPUTI',1,2,IAV,SCALE,UNI,NULO,STAT)
         IF ((SCALE(1).GE.-1).AND.(SCALE(1).LE.1)) SCALE(1) = 1
         IF ((SCALE(2).GE.-1).AND.(SCALE(2).LE.1)) SCALE(2) = 1
         IF (ACTION(3:3).EQ.'A') THEN
            SCALE(3) = 1
         ELSE
            SCALE(3) = 0
         ENDIF
C
C  calculate new nopix
         IF (SCALE(1).GT.1)  THEN
            NPIXC(1) = NPIX(1) * SCALE(1)
         ELSE IF (SCALE(1).LT.1) THEN
            NPIXC(1) = ( (NPIX(1)-1) / (-SCALE(1)) ) + 1
         ELSE
            NPIXC(1) = NPIX(1)
         ENDIF
C 
         IF (SCALE(2).GT.1)  THEN
            NPIXC(2) = NPIX(2) * SCALE(2)
            IF (NAXIS.EQ.1) THEN
               NAXIS = 2
               START(2) = START(1)
               STEP(2) = STEP(1)
               CUNIT(33:) = 'as in x-dim '
            ENDIF
         ELSE IF (SCALE(2).LT.1) THEN
            NPIXC(2) = ( (NPIX(2)-1) / (-SCALE(2)) ) + 1
         ELSE
            NPIXC(2) = NPIX(2)
         ENDIF
      write(*,123) naxis,npixc(1),npixc(2)
123   format('naxis now: ',i3,' npix = ',I5,i5)
         CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXIS,NPIXC,START,STEP,IDENT,
     +               CUNIT,PNTRC,IMNOC,STAT)
         CALL SCALI(SCALE,MADRID(PNTR),NPIX,MADRID(PNTRC))
         CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,1,4,UNI,STAT)
C
      ELSE IF (KFLAG.EQ.4) THEN                               !thinning ...
         CALL STIPUT(FRAMEC,D_I4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXIS,NPIX,START,STEP,IDENT,
     +               CUNIT,PNTRC,IMNOC,STAT)
C
C  first copy data
         SIZE = NPIX(1) * NPIX(2)
C        CALL COPYI(MADRID(PNTR),MADRID(PNTRC),SIZE)
C
C  get virtual memory for working buffer
         CALL STFXMP(SIZE,D_I4_FORMAT,WPNTR,STAT)
C
C  and do it
         CALL THINNI(MADRID(PNTRC),NPIX,MADRID(WPNTR))
C
      ELSE                                                   !correlation ...
         CALL STKRDC('IN_B',1,1,80,IAV,FRAMEB,UNI,NULO,STAT) !get spectrum
         CALL STKRDI('INPUTI',1,1,IAV,ISHI,UNI,NULO,STAT)    !get max. shift
C
         CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, !map spectrum
     +               2,NAXISB,NPIXB,STARTB,STEPB,IDENT,
     +               CUNIT,PNTRB,IMNOB,STAT)
         EPS = 0.001 * ABS(STEPB(1))                         ! 0.1 % of x-step
         DIFSTP = ABS(STEP(1)) - ABS(STEPB(1))
         IF (DIFSTP.GT.EPS)
     +      CALL STETER(1,
     +      'x-step of template and spec frame not equal...')
C
C  how many points will we have in overlap region?
         SIZE = NPIX(1) - ISHI - ISHI
         IF (SIZE.LT.3) THEN
            CALL STTPUT(' Too little overlap ',STAT)
            CALL STETER
     +      (2,'Shift too large or template frame too small...')
         ENDIF
C
         NAXIS = 1
         NPIXC(1) = 2*ISHI + 1
         START(1) = -(ISHI * STEPB(1))
         CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXIS,NPIXC,START,STEPB,IDENT,
     +               CUNIT,PNTRC,IMNOC,STAT)
C
C  and do it
         IF ((ACTION(2:2).EQ.'Y') .OR. (ACTION(2:2).EQ.'y')) THEN
            CALL CTCORR(NPIX(1),MADRID(PNTR),NPIXB(1),MADRID(PNTRB),
     +                  MADRID(PNTRC),ISHI,MINVAL,MAXVAL)
         ELSE
            CALL XSCORR(NPIX(1),MADRID(PNTR),NPIXB(1),MADRID(PNTRB),
     +                  MADRID(PNTRC),ISHI,MINVAL,MAXVAL)
         ENDIF
         LHCUTS(1) = MINVAL
         LHCUTS(2) = MAXVAL
         LHCUTS(3) = MINVAL
         LHCUTS(4) = MAXVAL
         CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,1,4,UNI,STAT) 
         
      ENDIF
C
C  that's it folks...
9000  CALL DSCUPT(IMNO,IMNOC,' ',STAT)         !copy all other descriptors
      CALL STSEPI
      END

      SUBROUTINE FLIPPI(A,R,NPIX,ACTION)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subrotine FLIPPI		version 1.0	870605
C  K. Banse                  	ESO - Garching
C
C.KEYWORDS
C  bulk data frame, flipping
C
C.PURPOSE
C  flip an image or mask in x and/or y dimension
C
C.ALGORITHM
C  straight forward, data is changed in place
C
C.INPUT/OUTPUT
C  call as FLIPPI(A,R,NPIX,ACTION)
C
C  A:		R*4 array		input image buffer
C  R:		R*4 array		buffer for single line
C  NPIX:	I*4 array		dimensions of A
C  ACTION:	char.exp		flip action
C					= 'X ', flip in x
C					  'Y ', flip in y
C					  'XY', flip in x and y
C
C.VERSIONS
C  1.00		from old FLIPPY
C	
C--------------------------------------------------
C
C	
      IMPLICIT NONE
C
      INTEGER     NPIX(2)
      INTEGER     FROM,GAP,GAPP,INDX,NTMP
      INTEGER     NX,NY,TO,XDIM,XMID,YDIM,YMID
C
      REAL        A(*),R(*)
      REAL        SAVE
C
      CHARACTER*(*)    ACTION
C
C  init
      XDIM = NPIX(1)
      YDIM = NPIX(2)
      XMID = XDIM/2
      YMID = YDIM/2
C
C  branch according to action
      IF (ACTION(1:1).EQ.'Y') GOTO 1000
      IF (ACTION(1:2).EQ.'XY') GOTO 2000
C
C  ACTION = 'X ', flip only in x
C
      FROM = 0                              !offset to first element in a row
      DO 500 NY=1,YDIM
         TO = FROM + XDIM + 1                        !point to last + 1
C
         DO 300 NX=1,XMID                            !loop till middle of row
            SAVE = A(NX+FROM)                        !flip...
            A(NX+FROM) = A(TO-NX)
            A(TO-NX) = SAVE
300      CONTINUE
         FROM = FROM + XDIM
500   CONTINUE
      RETURN
C
C  ACTION = 'Y ', flip only in y
C
1000  FROM = 0                             !offset to first element in a row
      DO 1500 NY=1,YMID                         !loop till middle of columns
         TO = YDIM + 1 - NY
         GAP = (TO-NY) * XDIM
         NTMP = FROM + GAP
C
         DO 1300 NX=1,XDIM                      !loop through row
            INDX = NX + FROM
            R(NX) = A(INDX)                     !save row elements in R
1300     CONTINUE
         DO 1330 NX=1,XDIM
            INDX = NX + FROM
            A(INDX) = A(INDX+GAP)               !flip complete row...
1330     CONTINUE
         DO 1360 NX=1,XDIM 
            INDX = NX + NTMP
            A(INDX) = R(NX)
1360     CONTINUE
         FROM = FROM + XDIM
1500  CONTINUE
      RETURN
C
C  ACTION = 'XY', flip in x and y
C
2000  FROM = 0
      DO 2500 NY=1,YMID                        !loop till middle of columns
         TO = YDIM + 1 - NY
         GAP = (TO-NY) * XDIM
         GAPP = GAP + XDIM + 1
         NTMP = FROM + GAPP
C
         DO 2300 NX=1,XDIM                      !loop through row
            INDX = NX + FROM
            R(NX) = A(INDX)                     !save row elements in R
2300     CONTINUE
         DO 2330 NX=1,XDIM
            INDX = NX + FROM
            A(INDX) = A(NTMP-NX)           !flip complete row...
2330     CONTINUE
         DO 2360 NX=1,XDIM 
            INDX = NTMP - NX
            A(INDX) = R(NX)
2360     CONTINUE
         FROM = FROM + XDIM
2500  CONTINUE
C
      IF (2*YMID.NE.YDIM) THEN                  !take care of middle row
         FROM = YMID * XDIM
         TO = FROM + XDIM + 1
         DO 2800 NX=1,XMID
            INDX = NX + FROM
            SAVE = A(INDX)
            A(INDX) = A(TO-NX)
            A(TO-NX) = SAVE
2800     CONTINUE
      ENDIF
C
      RETURN
      END

      SUBROUTINE SHIFTI(A,B,C,NPIX,XSHIFT,YSHIFT)
C
      IMPLICIT NONE
C
      INTEGER  NPIX(2),XSHIFT,YSHIFT
      INTEGER  NY,NX,NOFF,OFF1,OFF2,K1
C 
      REAL     A(*),B(*),C(*)
C 
      IF (XSHIFT.EQ.0) GOTO 3000
C 
      K1 = NPIX(1) - XSHIFT
      OFF1 = 0
C 
C  check, if we have to shift also y later on
      IF (YSHIFT.EQ.0) THEN
         DO 1800 NY=1,NPIX(2)                !shift in x for each line
            DO 1000 NX=1,K1
               NOFF = NX + OFF1
               C(NOFF+XSHIFT) = A(NOFF)
1000        CONTINUE
            DO 1500 NX=1,XSHIFT
               NOFF = NX + OFF1
               C(NOFF) = A(NOFF+K1)
1500        CONTINUE
            OFF1 = OFF1 + NPIX(1)
1800     CONTINUE
         RETURN
C 
      ELSE
         DO 2800 NY=1,NPIX(2)                !shift in x for each line
            DO 2000 NX=1,K1
               NOFF = NX + OFF1
               B(NOFF+XSHIFT) = A(NOFF)
2000        CONTINUE
            DO 2500 NX=1,XSHIFT
               NOFF = NX + OFF1
               B(NOFF) = A(NOFF+K1)
2500        CONTINUE
            OFF1 = OFF1 + NPIX(1)
2800     CONTINUE
      ENDIF
C 
3000  K1 = NPIX(2) - YSHIFT
C 
C  and now shift complete lines       (initially B = A )
      OFF1 = 0
      OFF2 = YSHIFT * NPIX(1)
      DO 4400 NY=1,K1
         DO 4000 NX=1,NPIX(1)
            C(OFF2+NX) = B(OFF1+NX)
4000     CONTINUE
         OFF1 = OFF1 + NPIX(1)               !move to next line
         OFF2 = OFF2 + NPIX(1)
4400  CONTINUE
C 
      OFF1 = 0
      OFF2 = K1 * NPIX(1)
      DO 5500 NY=1,YSHIFT
         DO 5000 NX=1,NPIX(1)
            C(OFF1+NX) = B(OFF2+NX)
5000     CONTINUE
         OFF1 = OFF1 + NPIX(1)               !move to next line
         OFF2 = OFF2 + NPIX(1)
5500  CONTINUE
C
      RETURN
      END

      SUBROUTINE SCALI(SCALE,IN,NPIX,OUT)
C	INTEGER SCALE(3) - x, yscale and method (1 = average, 0 = omit)
C	REAL IN(*) - input image
C	INTEGER NPIX(2) - dimension of input image
C	REAL OUT() - output image
C
      IMPLICIT NONE
C
      INTEGER  SCALE(3),NPIX(2),SIZE
      INTEGER  OUTNDX,INOFF,KSCAL,INDX
      INTEGER  NLINE,NY,NX,K,KK,NPX,M,JOUT
C
      REAL     IN(*),OUT(*)
      REAL     VAL,AVFACT
C
      OUTNDX = 1
      INOFF = 0
C
C  split according to y-SCALE
      IF (SCALE(2).LT.1) GOTO 10000
C
C  we keep (or increase) the no. of lines
C 
      IF (SCALE(1).GE.1) THEN                   !increase in x
C
         KSCAL = SCALE(1) - 1
         DO 2000, NLINE=1,NPIX(2)
            DO 1800, NY=1,SCALE(2)                    !replicate lines
               DO 1600, NX=1,NPIX(1)
                  INDX = INOFF + NX                   !point to input pixel
                  DO 1000, K=OUTNDX,OUTNDX+KSCAL      !replicate pixels
                     OUT(K) = IN(INDX)
1000              CONTINUE
                  OUTNDX = OUTNDX + SCALE(1)          !increment out pointer
1600           CONTINUE
1800        CONTINUE
            INOFF = INOFF + NPIX(1)                   !point to next line
2000     CONTINUE
C
      ELSE                                      !decrease in x
         SCALE(1) = -SCALE(1)
C 
         IF (SCALE(3).EQ.0) THEN                !we omit pixels
            DO 2300, NLINE=1,NPIX(2)
               DO 2200, NY=1,SCALE(2)                   !replicate lines
                  DO 2100, NX=1,NPIX(1),SCALE(1)
                     OUT(OUTNDX) = IN(INOFF+NX)
                     OUTNDX = OUTNDX + 1
2100              CONTINUE
2200           CONTINUE
               INOFF = INOFF + NPIX(1)                  !point to next line
2300        CONTINUE
C     
         ELSE                                  !averaging is harder...
            AVFACT = 1.0/SCALE(1)
            M = NPIX(1) / SCALE(1)
            NXLIM = M * SCALE(1)
            KSCAL = NPIX(1) - NXLIM            !nopts for last average on row
C 
            DO 3400,NLINE=1,NPIX(2)
               DO 3300,NY=1,SCALE(2)
                  DO 3200,NX=1,NXLIM,SCALE(1)
                     VAL = 0.
                     M = INOFF+NX-1
                     DO 3100,K=1,SCALE(1)
                        VAL = VAL + IN(M+K)
3100                 CONTINUE
                     OUT(OUTNDX) = VAL*AVFACT
                     OUTNDX = OUTNDX + 1
3200              CONTINUE
                  IF (KSCAL .GT. 0) THEN
                     VAL = 0.
                     M = INOFF+NXLIM
                     DO 3220,K=1,KSCAL               !sum up as much as we can
                        VAL = VAL + IN(M+K)
3220                 CONTINUE
                     DO 3240,K=1,SCALE(1)-KSCAL      !use last pixel for rest
                        VAL = VAL + IN(M+KSCAL)
3240                 CONTINUE
                     OUT(OUTNDX) = VAL*AVFACT
                     OUTNDX = OUTNDX + 1
                  ENDIF
3300           CONTINUE
C 
               INOFF = INOFF + NPIX(1)                  !point to next line 
3400        CONTINUE
         ENDIF
      ENDIF
      RETURN
C
C  we decrease the no. of lines
C 
10000 SCALE(2) = -SCALE(2)
      SIZE = NPIX(1) * NPIX(2)           !end of input array
      LASTROW = SIZE - NPIX(1)           !point to begin of last line
C 
      IF (SCALE(1).GE.1) THEN                   !increase in x
         KSCAL = SCALE(1) - 1
         KK = NPIX(1)*SCALE(1) - 1
C 
         IF (SCALE(3).EQ.0) THEN
            DO 12600, NLINE=1,NPIX(2),SCALE(2)
               DO 12500, NX=1,NPIX(1)
                  INDX = INOFF + NX                    !point to input pixel
                  DO 12400, K=OUTNDX,OUTNDX+KSCAL      !replicate pixels
                     OUT(K) = IN(INDX)
12400             CONTINUE
                  OUTNDX = OUTNDX + SCALE(1)           !increment out pointer
12500          CONTINUE
               INOFF = INOFF + (NPIX(1)*SCALE(2))      !point to next line
12600       CONTINUE
C
         ELSE                                  !averaging is harder...
            AVFACT = 1.0/SCALE(2)
            JOUT = 1
C 
            DO 13600, NLINE=1,NPIX(2),SCALE(2)
               DO 13100, K=OUTNDX,OUTNDX+KK            !init sums
                  OUT(K) = 0.0
13100          CONTINUE
C 
               DO 13500, NY=1,SCALE(2)
                  DO 13400, NX=1,NPIX(1)
                     INDX = INOFF + NX                 !point to input pixel
                     DO 13300, K=OUTNDX,OUTNDX+KSCAL   !sum up
                        OUT(K) = OUT(K) + IN(INDX)
13300                CONTINUE
                     OUTNDX = OUTNDX + SCALE(1)        !increment out pointer
13400             CONTINUE
                  OUTNDX = JOUT
                  INOFF = INOFF + NPIX(1)              !point to next line
                  IF (INOFF .GE. SIZE) INOFF = LASTROW  !don't fall off...
13500          CONTINUE
C 
               DO 13550, K=OUTNDX,OUTNDX+KK
                  OUT(K) = OUT(K) * AVFACT
13550          CONTINUE
               JOUT = JOUT + (NPIX(1)*SCALE(1))
               OUTNDX = JOUT
13600       CONTINUE
         ENDIF
C 
      ELSE                                      !decrease in x
         SCALE(1) = -SCALE(1)
         IF (SCALE(3).EQ.0) THEN
            DO 14400, NLINE=1,NPIX(2),SCALE(2)
               DO 14200, NX=1,NPIX(1),SCALE(1)
                  OUT(OUTNDX) = IN(INOFF+NX)
                  OUTNDX = OUTNDX + 1
14200          CONTINUE
               INOFF = INOFF + (NPIX(1)*SCALE(2))
14400       CONTINUE
C 
C average in x and y
C 
         ELSE                                  !averaging is harder...
            JOUT = 1
            AVFACT = 1.0/(SCALE(1)*SCALE(2))
            M = NPIX(1) / SCALE(1)
            NXLIM = M * SCALE(1)
            KSCAL = NPIX(1) - NXLIM            !nopts for last average on row
C 
            DO 16000, NLINE=1,NPIX(2),SCALE(2)
               DO 15100, NPX=1,NPIX(1),SCALE(1)
                  OUT(OUTNDX) = 0.0
                  OUTNDX = OUTNDX + 1
15100          CONTINUE
C
               OUTNDX = JOUT
               DO 15800,NY=1,SCALE(2)
                  DO 15700,NPX=1,NXLIM,SCALE(1)
                     INDX = INOFF + NPX - 1
                     DO 15600, NX=1,SCALE(1)
                        OUT(OUTNDX) = OUT(OUTNDX) + IN(INDX+NX)
15600                CONTINUE
                     OUTNDX = OUTNDX + 1
15700             CONTINUE
                  IF (KSCAL .GT. 0) THEN        !fill the last aver.sum
                     INDX = INOFF + NXLIM 
                     DO 15740, NX=1,KSCAL             !sum up as much as we can
                        OUT(OUTNDX) = OUT(OUTNDX) + IN(INDX+NX)
15740                CONTINUE
                     DO 15770, NX=1,SCALE(1)-KSCAL    !use last pixel for rest
                        OUT(OUTNDX) = OUT(OUTNDX) + IN(INDX+KSCAL)
15770                CONTINUE
                  ENDIF
                  OUTNDX = JOUT
                  INOFF = INOFF + NPIX(1)             !follow with in pointer
                  IF (INOFF .GE. SIZE) INOFF = LASTROW  !don't fall off...
15800          CONTINUE
C 
               DO 15900, NPX=1,NPIX(1),SCALE(1)
                  OUT(OUTNDX) = OUT(OUTNDX)*AVFACT
                  OUTNDX = OUTNDX + 1
15900          CONTINUE
               JOUT = OUTNDX                      !follow with base index
16000       CONTINUE
         ENDIF
      ENDIF
C
      RETURN
      END

      SUBROUTINE THINNI(M,NPIX,W)
C
      IMPLICIT NONE
C
      INTEGER    M(*),NPIX(2),W(*)
      INTEGER    REMAIN,INC,IOFF,IOFFX
      INTEGER    MX,MY,JN,NX,NY
      INTEGER    JOFF(4),N1,N2,N3,N4,N5,N6,N7,N8,KOFF
      INTEGER    IOFA
C
      NX = NPIX(1)
      NY = NPIX(2)
      JOFF(1) = 1
      JOFF(2) = NX
      JOFF(3) = -1
      JOFF(4) = -NX
      REMAIN = 1
C
C  main loop
200   IF (REMAIN.NE.1) RETURN                !We're done.
C
      REMAIN = 0
C
      DO 5000, JN=1,4                         !look at all direct neighbours...
         KOFF = JOFF(JN)
         INC = 0
         IOFFX = NX
         DO 4000, MY=2,NY-1                   !omit first + last line...
            DO 3000, MX=2,NX-1                !omit first + last pixel...
               IOFF = IOFFX + MX
               IF ( (M(IOFF).EQ.1) .AND.
     +              (M(IOFF+KOFF).EQ.0) ) THEN       !get all neighbours...
                  IOFA = IOFF + NX
                  N1 = M(IOFF+1)
                  N2 = M(IOFA+1)
                  N3 = M(IOFA)
                  N4 = M(IOFA-1)
                  IOFA = IOFF - NX
                  N5 = M(IOFF-1)
                  N6 = M(IOFA-1)
                  N7 = M(IOFA)
                  N8 = M(IOFA+1)
C
C  now check conditions
                  IF (N1.EQ.0) THEN
C
                     IF (N5.EQ.0) THEN                           ! 0AAA0BBB
                        IF ((N2+N3+N4.GT.0).AND.(N6+N7+N8.GT.0))
     +                     GOTO 1000
                     ELSE IF ((N7.EQ.0).AND.(N8.EQ.2)) THEN      ! 0AAAAA02
                        IF (N2+N3+N4+N5+N6.GT.0) GOTO 1000
                     ELSE IF ((N3.EQ.0).AND.(N2.EQ.2)) THEN      ! 020AAAAA
                        IF (N4+N5+N6+N7+N8.GT.0) GOTO 1000
                     ENDIF
                  ENDIF
C
                  IF (N3.EQ.0) THEN
C
                     IF (N7.EQ.0) THEN                           ! BB0AAA0B
                        IF ((N1+N2+N8.GT.0).AND.(N4+N5+N6.GT.0))
     +                     GOTO 1000
                     ELSE IF ((N5.EQ.0).AND.(N4.EQ.2)) THEN      ! AA020AAA
                        IF (N1+N2+N6+N7+N8.GT.0) GOTO 1000
                     ENDIF
                  ENDIF
C
                  IF ((N5+N7.EQ.0).AND.(N6.EQ.2)) THEN           ! AAAA020A
                     IF (N1+N2+N3+N4+N8.GT.0) GOTO 1000
                  ENDIF
C
                  M(IOFF) = 3              !no pattern match...
                  REMAIN = 1
                  INC = INC + 1
                  W(INC) = IOFF            !save index...
                  GOTO 3000
C
1000              M(IOFF) = 2              !yes, pattern match...
               ENDIF
C
3000        CONTINUE                                   !end of MX loop...
C
            IOFFX = IOFFX + NX
4000     CONTINUE                                      !end of MY loop...
C
C  change all 3's to 0's
         IF (INC.GT.0) THEN
            DO 4400, MX=1,INC
               M(W(MX)) = 0
4400        CONTINUE
         ENDIF
C
5000  CONTINUE                                         !end of JN loop...
C
C  end of main loop
      GOTO 200
C
      END


      SUBROUTINE XSCORR(NTEMP,TEMP,NSPEC,SPEC,RES,NSH,MINVAL,MAXVAL)
C
C    This program computes the cross-correlation between a template
C    and a spectrum over a bandwith of 2*NSH pixels. Identity of
C    start and step are assumed. We will have 2*NSH + 1 data.
C
C
      IMPLICIT NONE
C
      INTEGER   NTEMP,NSPEC,NSH
      INTEGER   MPS,MPE,INC,I,J
C
      REAL      TEMP(*),SPEC(*),RES(*)
      REAL      MINVAL,MAXVAL
C
      DOUBLE PRECISION  SUM
C
      MPS = 1 + NSH
      IF (NSPEC .GT. NTEMP) THEN                       !watch out for that
         MPE = NTEMP - NSH
      ELSE
         MPE = NSPEC - NSH
      ENDIF
      IF (MPE.LT.MPS)
     +   CALL STETER(2,'Shift too large or frames too small...')
      INC = 1
C
      DO 20, I=-NSH,NSH
         SUM = 0.D0
         DO 10, J=MPS,MPE
            SUM = SUM + ( SPEC(J)*TEMP(J+I) )
10       CONTINUE
C
         RES(INC) = SNGL(SUM)
         INC = INC + 1
20    CONTINUE
C
C  get min, max
      MINVAL = RES(1)
      MAXVAL = MINVAL
      DO 100, I=2,INC-1
         IF (RES(I).GT.MAXVAL) THEN
            MAXVAL = RES(I)
         ELSE IF (RES(I).LT.MINVAL) THEN
            MINVAL = RES(I)
         ENDIF
100   CONTINUE
C
      RETURN
      END


      SUBROUTINE CTCORR(NTEMP,TEMP,NSPEC,SPEC,RES,NSH,MINVAL,MAXVAL)
C
C    This program computes the cross-correlation between a template
C    and a spectrum over a bandwith of 2*NSH pixels. Identity of
C    start and step are assumed. We will have 2*NSH + 1 data.
C							M. Rosa
C
CT   Modified to compute true correlation coefficient with optimal
CT   overlap range.
CT   For two arrays x and y the correlation coefficient is
CT      r = Covar(x,y)/Sqrt[Var(x)*Var(y)], 
CT   which gives
CT      r = [n*sum(x*y) - sum(x)*sum(y)]/
CT          sqrt{[n*sum(x**2)-sum(x)**2]*[n*sum(y**2)-sum(y)**2]}
CT   Since r is normalized one can use the total overlap of two
CT   shifted arrays, varying with the shift variable.
CT   29.11.2001                                         H.-C. Thomas
CT					Max-Planck Institute for Astrophysics
C
C
      IMPLICIT NONE
C
      INTEGER   NTEMP,NSPEC,NSH
      INTEGER   MPS,MPE,INC,I,J
C
      REAL      TEMP(*),SPEC(*),RES(*)
      REAL      MINVAL,MAXVAL
C 
      DOUBLE PRECISION  SUMX,SUMY,SUMXX,SUMYY,SUMXY,NPOI,T1,T2,T3
C
CT                             compute the smallest overlap, exit if less than 1
CT
      MPS = 1
      IF (NSPEC .GT. NTEMP) THEN                       !watch out for that
         MPE = NTEMP - NSH 
      ELSE
         MPE = NSPEC - NSH
      ENDIF
      IF (MPE.LT.MPS) 
     +   CALL STETER(2,'Shift too large or frames too small...')
      INC = 1
C
CT       compute the necessary sums over a range MPS to MPE varying with shift I
CT
      DO 20, I=-NSH,NSH
         SUMX = 0.D0
         SUMXX = 0.D0
         SUMY = 0.D0
         SUMYY = 0.D0
         SUMXY = 0.D0
         MPS = 1
         IF (I .GT. 0) MPS = 1+I
         MPE = NSPEC+I
         IF (MPE .GT. NTEMP) MPE = NTEMP
         DO 10, J=MPS,MPE
            SUMX = SUMX + TEMP(J)
            SUMXX = SUMXX + TEMP(J)**2
            SUMY = SUMY + SPEC(J-I)
            SUMYY = SUMYY + SPEC(J-I)**2
            SUMXY = SUMXY + ( TEMP(J)*SPEC(J-I) )
10       CONTINUE
C
CT                                           compute the correlation coefficient
CT
         NPOI = DFLOAT(1+MPE-MPS)
         T1 = NPOI*SUMXY-SUMX*SUMY
         T2 = NPOI*SUMXX-SUMX**2
         T3 = NPOI*SUMYY-SUMY**2
         RES(INC) = SNGL(T1/DSQRT(T2*T3))
         INC = INC + 1
20    CONTINUE
C
C  get min, max
      MINVAL = RES(1)
      MAXVAL = MINVAL
      DO 100, I=2,INC-1
         IF (RES(I).GT.MAXVAL) THEN
            MAXVAL = RES(I)
         ELSE IF (RES(I).LT.MINVAL) THEN
            MINVAL = RES(I)
         ENDIF
100   CONTINUE
C
      RETURN
      END
