C  @(#)genxx1.for	19.1 (ESO-DMD) 02/25/03 14:01:43
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 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 GENXX1
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program  GENXX1			 version 1.00	890704
C  K. Banse                 		 ESO - Garching
C  1.10	891124		
C
C.KEYWORDS
C  image processing functions
C
C.PURPOSE
C  merger of several functions, see description there
C
C.ALGORITHM
C  see local descriptions
C
C.INPUT/OUTPUT
C  the following keys are used:
C
C  ACTION/C/1/2		action flag
C			= ED, for edge detection
C			= EX, for general extract routines
C			= EZ, for reference frame extraction
C			= IN, for inserting an image
C			= MP, for mapping images
C			= MA, for matrix operations
C
C.VERSION
C  1.00		merge former individual main programmes
C
C 021212	last modif
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER    IAV,STAT
      INTEGER    UNI(1),NULO
C
      CHARACTER  ACTION*2
C
C  set up MIDAS environment + enable automatic error abort
      CALL STSPRO('GENXX1')
C
C  get option
      CALL STKRDC('ACTION',1,1,2,IAV,ACTION,UNI,NULO,STAT)
C 
      IF (ACTION.EQ.'ED') THEN
         CALL SUBEDG
      ELSE IF (ACTION.EQ.'EX') THEN
         CALL SUBEXT
      ELSE IF (ACTION.EQ.'EZ') THEN
         CALL SUBEZT
      ELSE IF (ACTION.EQ.'IN') THEN
         CALL SUBINS
      ELSE IF (ACTION.EQ.'MP') THEN
         CALL SUBMAP
      ELSE IF (ACTION.EQ.'MA') THEN
         CALL SUBMAT
      ENDIF
C
C  get away 
      CALL STSEPI
      END

      SUBROUTINE SUBEDG
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine SUBEDG         	 version 1.00	881028
C  K. Banse                 	 ESO - Garching
C
C.KEYWORDS
C  edge detection
C
C.PURPOSE
C  detect edges in a given data frame
C
C.ALGORITHM
C  get the names of input + output frames from IN_A + OUT_A
C  get the threshold for algorithm from INPUTR
C  and calculate edge pixels via subroutine FNDEDGE (see description
C  of algorithm there)
C
C.INPUT/OUTPUT
C  the following keys are used:
C
C  IN_A/C/1/80			input data frame
C  P4/C/1/1			edge detection method
C				= B, use band - method
C				= S, use Sobel operator
C  OUT_A/C/1/80			output data frame
C  DEFAULT/C/1/1		default flag, if = Y, use half of dynamic range of input frame
C				as threshold, else read it from INPUTR
C  INPUTR/R/1/1			threshold for edge detection
C
C.VERSION
C  1.00		from version 2.00 as of 841120
C		move to FORTRAN 77 + new ST interfaces
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER    NAXISA,IAV,STAT
      INTEGER*8  PNTRA,PNTRC
      INTEGER    IMNOA,IMNOC
      INTEGER    NPIXA(3)
      INTEGER    UNI(1),NULO,MADRID(1)
C
      CHARACTER*80   FRAMEA,FRAMEC
      CHARACTER      CUNITA*64,IDENTA*72
      CHARACTER      DEFAUL*1,ACTION*1
C
      REAL       THRESH,LHCUTS(2)
      REAL       OLDCUT(2)
C
      DOUBLE PRECISION STEPA(3),STARTA(3)
C
      COMMON /VMR/ MADRID
C
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C
      DATA       LHCUTS    /0.,1./
      DATA       IDENTA    /' '/,   CUNITA  /' '/
C
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C
C  get input frame + map it
      CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT)
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            3,NAXISA,NPIXA,STARTA,STEPA,
     +            IDENTA,CUNITA,PNTRA,IMNOA,STAT)
C
C  get output frame + map it (=edge file) + copy standard descriptors
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT)
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXISA,NPIXA,STARTA,STEPA,
     +            IDENTA,CUNITA,PNTRC,IMNOC,STAT)
C
C  get threshold + method
      CALL STKRDC('DEFAULT',1,1,1,IAV,DEFAUL,UNI,NULO,STAT)
      IF ( (DEFAUL.EQ.'Y') .OR. (DEFAUL.EQ.'y') ) THEN
         CALL STDRDR(IMNOA,'LHCUTS',3,2,IAV,OLDCUT,UNI,NULO,STAT)
         THRESH = (OLDCUT(2) - OLDCUT(1)) * .5
      ELSE
         CALL STKRDR('INPUTR',1,1,IAV,THRESH,UNI,NULO,STAT)
      ENDIF
      CALL STKRDC('P4',1,1,1,IAV,ACTION,UNI,NULO,STAT)
      CALL UPCAS(ACTION,ACTION)
C
C  now do the edge detection
      CALL FNDEDG(ACTION,MADRID(PNTRA),NPIXA(1),NPIXA(2),
     +            THRESH,MADRID(PNTRC))
      CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,3,2,UNI,STAT)
C
      RETURN
      END

      SUBROUTINE SUBEXT
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine SUBEXT			version 1.00	881028
C  K. Banse            		 	ESO - Garching
C  981016
C
C.KEYWORDS
C  bulk data frames
C
C.PURPOSE
C  extract a subframe from a given frame
C
C.ALGORITHM
C  get output frame from OUT_A
C  get relevant coordinates as character data from INPUTC or via cursor
C
C.INPUT/OUTPUT
C  the following keys are used:
C  
C  IN_A/C/1/80		input frame
C  OUT_A/C/1/80		result frame
C  ACTION/C/3/2		action flag
C  action flag = IM: extract subframe from coord interval
C              = IC: extract subframe around center coord
C              = LI: extract a line in any direction from a higher
C                    dimensioned frame
C  if action flag = IM
C  INPUTC/C/1/80	character string of coordinates
C			in the form: [x1,y1,z1:x2,y2,z2] 
C			x,y,z according to the MIDAS standard for 
C			specifying coordinates
C  or if centered extract
C  INPUTC/C/1/80	[xc,yc,zc]
C  INQXI/I/I/1/3	xl,yl,zl        (0,...,NPIX)
C  INQXI/I/I/4/3	xr,yr,zr        (0,...,NPIX)

C  if action flag = LI
C  INPUTC/C/1/80	character string of coordinates
C			in the form: [x1,y1,z1:x2,y2,z2] 
C  INPUTD/D/1/1		stepsize for line
C	
C.VERSIONS
C 
C 021212	last modif	25 years ...
C 
C--------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER      IAV,IAX,N,NINDX,NN,REFSET
      INTEGER      NAXISA,NAXISC,STAT
      INTEGER*8    PNTRA,PNTRC,PNTRX,PNTRY
      INTEGER      NPIXA(3),NPIXC(3)
      INTEGER      RETBUF(3),SUBLO(3),SUBHI(3)
      INTEGER      MADRID(1),UNI(1),NULO
      INTEGER      IMNOA,IMNOC
      INTEGER      NPXY,K2,RMAIND,CHUCHO,CHUCHI
      INTEGER      NOLIN,TRUSIZ,CHUNKI,CHUNKO,K1,KOFF,MOFF,LOOPLM
      INTEGER      PXOFF(6),PCENT(3),CLONY
C
      CHARACTER*80   FRAMEA,FRAMEC,EQUAL
      CHARACTER      CUNITA*64,ACTION*2
      CHARACTER      INPUTC*80
      CHARACTER      ERROR1*40,ERROR2*40,ERROR3*40
      CHARACTER*72   IDENTA
C      
      DOUBLE PRECISION STARTA(3),STARTC(3),STEPA(3),STEPC(3),STEPLN
      DOUBLE PRECISION DIN(4),DOUT(4)
C 
      REAL           CUTS(4),RM,DIST,WW(2)
      REAL           XA,XB,YA,YB
C 
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C 
      INCLUDE  'MID_INCLUDE:IDIDEV.INC'
      INCLUDE  'MID_INCLUDE:IDIMEM.INC'
C 
      COMMON   /VMR/  MADRID
C
      DATA      ERROR1   /'invalid coordinate string ... '/
      DATA      ERROR2   /'problems with dynamic range... '/
      DATA      ERROR3   /'max. 3-dim images supported... '/
      DATA      STEPA    /3*1.D0/, STARTA /3*0.D0/
      DATA      STEPC    /3*1.D0/, STARTC /3*0.D0/
      DATA      NPIXA    /3*1/,    NPIXC /3*1/
      DATA      SUBLO    /1,1,1/,  SUBHI /1,1,1/
      DATA      CUTS     /0.,0.,0.,0./
      DATA      IDENTA   /' '/
      DATA      CUNITA   /' '/
      DATA      PXOFF    /6*0/
C 
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C
C  get name of input frame, result frame, descr_copy_flag + action flag
      CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT)
      CALL CLNFRA(FRAMEA,FRAMEA,0)
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT)
      CALL STKRDC('ACTION',1,3,2,IAV,ACTION,UNI,NULO,STAT)
      CLONY = 0
C      
C  get string with coordinates
      REFSET = 0                     !set to 1, if we update REFPIX descr.
      CALL STKRDC('INPUTC',1,1,80,IAV,INPUTC,UNI,NULO,STAT)
C
C  check if we shrink same frame
      CALL GENEQF(FRAMEA,FRAMEC,STAT)
      IF (STAT.EQ.1) THEN                !if same name, use temporary one
         EQUAL(1:) = FRAMEC(1:)
         FRAMEC(1:) = 'middummeq.bdf '
      ELSE
         EQUAL(1:) = '  '
      ENDIF
C      
C  and open input frame or open + already map it
      IF (ACTION(1:1).EQ.'I') THEN
         CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT)
         CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXISA,UNI,NULO,STAT)
         IF (NAXISA.GT.3) CALL STETER(44,ERROR3)
         CALL STDRDI(IMNOA,'NPIX',1,NAXISA,IAV,NPIXA,UNI,NULO,STAT)
         CALL STDRDD(IMNOA,'START',1,NAXISA,IAV,STARTA,UNI,NULO,STAT)
         CALL STDRDD(IMNOA,'STEP',1,NAXISA,IAV,STEPA,UNI,NULO,STAT)
         CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENTA,UNI,NULO,STAT)
         N = (NAXISA+1) * 16
         CALL STDRDC(IMNOA,'CUNIT',1,1,N,IAV,CUNITA,UNI,NULO,STAT)
C  
C  check if centered extract
         IF (ACTION(2:2).EQ.'C') THEN             !Yes, centered extraction
            CALL STKRDI('INQXI',1,6,IAV,PXOFF,UNI,NULO,STAT)
            CALL EXTCO1(IMNOA,INPUTC,NAXISA,NAXISC,PCENT,STAT)
            IF (STAT.NE.0) CALL STETER(1,ERROR1)
            SUBLO(1) = PCENT(1) - PXOFF(1)
            SUBHI(1) = PCENT(1) + PXOFF(4)
            IF (NAXISC.GT.1) THEN
               SUBLO(2) = PCENT(2) - PXOFF(2)
               SUBHI(2) = PCENT(2) + PXOFF(5)
               IF (NAXISC.GT.2) THEN
                  SUBLO(3) = PCENT(3) - PXOFF(3)
                  SUBHI(3) = PCENT(3) + PXOFF(6)
               ENDIF
            ENDIF
            DO 440, N=1,NAXISC
               IF ( (SUBLO(N).LT.1) .OR.
     +              (SUBLO(N).GT.SUBHI(N)) .OR.
     +              (SUBHI(N).GT.NPIXA(N)) )
     +            CALL STETER(4,'invalid pixel offsets ... ')
440         CONTINUE
            GOTO 480
         ENDIF
      ELSE
         CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,3,
     +            NAXISA,NPIXA,STARTA,STEPA,
     +            IDENTA,CUNITA,PNTRA,IMNOA,STAT)
         IF (NAXISA.GT.3) CALL STETER(44,ERROR3)
      ENDIF
      CALL FPXWCO(0,IMNOA,DIN,DOUT,STAT)               !init wcs conversion
C 
C  copy subframe or extract a line
      CALL EXTCOO(IMNOA,INPUTC,NAXISA,NAXISC,SUBLO,SUBHI,STAT)
      IF (STAT.NE.0) CALL STETER(1,ERROR1)
C      
C  get start + end pixels of subframe within father frame 
480   DO 500 N=1,NAXISC
         DIN(N) = SUBLO(N)                             !in fp's
         NPIXC(N) = SUBHI(N) - SUBLO(N) + 1
500   CONTINUE
      CALL FPXWCO(1,0,DIN,DOUT,STAT)                   !fp -> wc
      DO 510 N=1,NAXISC
         STARTC(N) = STARTA(N) + (SUBLO(N)-1)*STEPA(N)
510   CONTINUE
      CALL STDFND(IMNOA,'REFPIX',INPUTC,N,NN,STAT)
      IF (INPUTC(1:1).EQ.'D') THEN
         CALL STDRDD(IMNOA,'REFPIX',1,NAXISA,IAV,DOUT,UNI,NULO,STAT)
         DO 8080, N=1,NAXISA
            DOUT(N) =  DOUT(N) + 1.0D0 - DIN(N)
8080     CONTINUE
         REFSET = 1
      ENDIF
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT)
C 
C  branch according to line or subframe extraction...
      IF (ACTION(1:1).EQ.'L') GOTO 5000
C  
C -------------------------------------------------------
C  
C  here to extract a subframe
C      
C -------------------------------------------------------
C  
      TRUSIZ = 1                           !get size of input frame
      DO 600 N=1,NAXISA
         TRUSIZ = TRUSIZ * NPIXA(N)
600   CONTINUE
C 
C  get buffersize of at least one line but no more than total size
      CALL HACKUP(NPIXA,D_R4_FORMAT,RETBUF)
      NOLIN = RETBUF(1)                      !buffer = NOLIN * NPIXA(1)
      IF (NOLIN.GT.NPIXC(2)) NOLIN = NPIXC(2)
      CHUNKI = NOLIN * NPIXA(1)
      CHUNKO = NOLIN * NPIXC(1)
      LOOPLM = NPIXC(2) / NOLIN              !get no. of passes through loop
      RMAIND = NPIXC(2) - (LOOPLM*NOLIN)     !and remainder
C 
      CALL STFXMP(CHUNKI,D_R4_FORMAT,PNTRX,STAT)
      CALL STFXMP(CHUNKO,D_R4_FORMAT,PNTRY,STAT)
C 
      TRUSIZ = 1
      DO 800 N=1,NAXISC
         TRUSIZ = TRUSIZ * NPIXC(N)
800   CONTINUE
      CALL GENBLD(IMNOA,FRAMEC,D_R4_FORMAT,TRUSIZ,IMNOC,CLONY,STAT)
C  
C  now we loop through input frame
      NPXY = NPIXA(1) * NPIXA(2)
      K2 = (SUBLO(3)-1) * NPXY                   !skip unused planes
      K1 = (SUBLO(2)-1) * NPIXA(1) + 1           !point to 1. used y-lines
      MOFF = 1                                   !pointer for output frame
C 
      CUTS(3) = 999999.
      CUTS(4) = -999999.
C 
      DO 4000 IAX=1,NPIXC(3)                     !loop over all planes
         KOFF = K2 + K1
         IF (LOOPLM.LE.0) GOTO 3000
C 
         DO 2000 N=1,LOOPLM
            CALL STFGET(IMNOA,KOFF,CHUNKI,TRUSIZ,MADRID(PNTRX),STAT)  !data in
            CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),SUBLO(1),1,
     +               NPIXA(1),NOLIN,NPIXA(1),NPIXC(1),CUTS(3))   !rearrange data
            CALL STFPUT(IMNOC,MOFF,CHUNKO,MADRID(PNTRY),STAT)         !data out
            KOFF = KOFF + CHUNKI                     !follow with indices...
            MOFF = MOFF + CHUNKO
2000     CONTINUE
C 
3000     IF (RMAIND.GT.0) THEN
            CHUCHI = RMAIND * NPIXA(1)
            CHUCHO = RMAIND * NPIXC(1)
            CALL STFGET(IMNOA,KOFF,CHUCHI,TRUSIZ,MADRID(PNTRX),STAT) 
            CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),SUBLO(1),1,
     +               NPIXA(1),RMAIND,NPIXA(1),NPIXC(1),CUTS(3))      
            CALL STFPUT(IMNOC,MOFF,CHUCHO,MADRID(PNTRY),STAT)  
            MOFF = MOFF + CHUCHO
         ENDIF
C 
         K2 = K2 + NPXY                                !move to next plane
4000  CONTINUE
C 
C  compress NAXISC + write descriptors
4200  IF (NAXISC.GT.1) THEN
         IF (NPIXC(1).EQ.1) THEN
            NPIXC(1) = NPIXC(2)
            NPIXC(2) = NPIXC(3)
            NPIXC(3) = 1
            STARTC(1) = STARTC(2)
            STARTC(2) = STARTC(3)
            STARTC(3) = 0.D0
            STEPA(1) = STEPA(2)
            STEPA(2) = STEPA(3)
            STEPA(3) = 1.D0
            CUNITA(17:) = CUNITA(33:)//' '
            NAXISC = NAXISC - 1
            GOTO 4200
         ENDIF
         IF (NPIXC(2).EQ.1) THEN
            NPIXC(2) = NPIXC(3)
            NPIXC(3) = 1
            STARTC(2) = STARTC(3)
            STARTC(3) = 0.D0
            STEPA(2) = STEPA(3)
            STEPA(3) = 1.D0
            CUNITA(33:) = CUNITA(49:)//' '
            NAXISC = NAXISC - 1
            GOTO 4200
         ENDIF
         IF (NPIXC(NAXISC).EQ.1) NAXISC = NAXISC - 1
      ENDIF
      IF (NAXISC.LT.3) NPIXC(3) = 1
      IF (NAXISC.LT.2) NPIXC(2) = 1
C 
      CALL STDWRI(IMNOC,'NAXIS',NAXISC,1,1,UNI,STAT)
      CALL STDWRI(IMNOC,'NPIX',NPIXC,1,NAXISA,UNI,STAT)      !still use NAXISA
      CALL STDWRD(IMNOC,'START',STARTC,1,NAXISA,UNI,STAT)
      CALL STDWRD(IMNOC,'STEP',STEPA,1,NAXISA,UNI,STAT)
C 
      GOTO 8000
C  
C -------------------------------------------------------
C  
C  here to extract a line
C      
C -------------------------------------------------------
C  
5000  IF ( (NAXISA.GT.2) .AND.         !test that input frame is a 2-d frame...
     +     (NPIXA(3).GT.1) ) THEN
         NN = INDEX(FRAMEA,' ') - 1
         IF (NN.LE.0) NN = LEN(FRAMEA)
         INPUTC(1:) = 'Input frame '//FRAMEA(1:NN)//
     +                ' has to be 2-dim frame... '
         CALL STETER(3,INPUTC)
      ENDIF
C 
C  get stepsize for line
      CALL STKRDD('INPUTD',1,1,IAV,STEPLN,UNI,NULO,STAT)
      XA = SUBLO(1)
      YA = SUBLO(2)
      XB = SUBHI(1)
      YB = SUBHI(2)
C      
C  determine main axis (depends upon orientation of line in real space)
      IF (ABS(XA-XB).LT.10.E-5) THEN
         IAX = 2                                 !parallel to y-axis...
      ELSE
         RM = (YB-YA) / (XB-XA)
         IF (ABS(RM).GT.1.0001) THEN
            IAX = 2                              !steep ramp, so we use y..
         ELSE
            IAX = 1                              !low ramp, so we use x...
         ENDIF
      ENDIF
      IF (ABS(STEPLN).LE.10.E-20) THEN        !use orig_step if = 0.0
         STEPLN = STEPA(IAX)
         RM = 1.                              !RM holds scaled stepsieze
      ELSE
         RM = STEPLN/STEPA(IAX)
      ENDIF
C  
C  get no. of points and map index arrays as dummy frame
      WW(1) = XB - XA     
      WW(2) = YB - YA
      DIST = SQRT( (WW(1)*WW(1)) + (WW(2)*WW(2)) )
      WW(1) = DIST/RM
      IF (WW(1).LT.0.) WW(1) = -WW(1)
      NN = IFIX( WW(1) )
      WW(1) = 0.01 * WW(1)
      IF (WW(1).GT.3.0) THEN
         NN = NN + IFIX( WW(1) )                  !add 1 %  for safety check
      ELSE
         NN = NN + 3                                !add 3 for safety check
      ENDIF
C      
      CALL STFXMP(NN,D_R4_FORMAT,PNTRX,STAT)
      CALL STFXMP(NN,D_R4_FORMAT,PNTRY,STAT)
C      
C  extract the indices of the complete line from input frame
      CALL PIXLIN(XA,YA,XB,YB,RM,MADRID(PNTRX),MADRID(PNTRY),NN,NINDX)
      IF (NINDX.GE.NN) 
     +   CALL STETER(3,'too many elements in index array...')
C  
C  create 1-dim output file + get z-values for pixels recorded in XINDX + YINDX
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            1,NINDX,STARTC,STEPLN,
     +            IDENTA,CUNITA,PNTRC,IMNOC,STAT)
      CALL ZIMA(MADRID(PNTRA),NPIXA,MADRID(PNTRX),MADRID(PNTRY),NINDX,
     +          MADRID(PNTRC),CUTS(3),CUTS(4))
C 
C  also save the Y-start coordinate
      STARTC(2) = STARTA(2) + (SUBLO(2)-1)*STEPA(2)
      CALL STDWRD(IMNOC,'START',STARTC,1,2,UNI,STAT)  !STIPUT only fills 1 value
C 
C  common section
C  
C  copy non-standard descriptors, update LHCUTS + append HISTORY
8000  IF (CLONY.EQ.0) THEN
         CALL DSCUPT(IMNOA,IMNOC,' ',STAT)
      ELSE
         CALL DSCUPT(IMNOC,IMNOC,' ',STAT)
      ENDIF
      CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT)
C
      IF (REFSET.EQ.1) THEN          !don't forget to update REFPIX
         CALL STDWRD(IMNOC,'REFPIX',DOUT,1,NAXISA,UNI,STAT)
      ENDIF
C
      IF (EQUAL(1:).NE.' ') THEN
         CALL STFCLO(IMNOA,STAT)
         CALL STFCLO(IMNOC,STAT)
         NN = INDEX(EQUAL,' ') - 1
         IF (NN.LT.1) NN = LEN(EQUAL)
         IDENTA(1:) = ' '                   !IDENTA is at least 4 chars. longer
         IDENTA(1:) = EQUAL(1:)
         DO 8400 N=NN,1,-1
            IF (EQUAL(N:N).EQ.'.') GOTO 8800
            IF ((EQUAL(N:N).EQ.'/') .OR.
     +          (EQUAL(N:N).EQ.']')) GOTO 8440
8400     CONTINUE
8440     IDENTA(NN+1:) = '.bdf'
8800     CALL STFRNM(FRAMEC,IDENTA,STAT)
      ENDIF
C      
      RETURN
      END

      SUBROUTINE SUBEZT
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine SUBEZT		version 1.00	881028
C  K. Banse			ESO - Garching
C
C.KEYWORDS
C  bulk data frames
C
C.PURPOSE
C  extract only "referenced" pixels from a frame
C
C.ALGORITHM
C  Pull out only those pixels from inframe which are indicated in
C  a reference frame. Result frame will be a 1-dim frame.
C
C.INPUT/OUTPUT
C  the following keys are used:
C
C  IN_A/C/1/80			input frame
C  IN_B/C/1/80			reference frame
C  OUT_A/C/1/80			result frame
C  THR/R/1/1			threshold value in reference frame
C
C.VERSIONS
C  1.0		from version 1.00 as of 860110
C               move to FORTRAN 77 + new ST-Interfaces
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER    IAV,N,NAXISA,NAXISB
      INTEGER*8  PNTRA,PNTRB,PNTRC
      INTEGER    IMNOA,IMNOB,STAT
      INTEGER    SIZE,TRUESZ
      INTEGER    NPIXA(3),NPIXB(3)
      INTEGER    UNI(1),NULO,MADRID(1)
C
      CHARACTER*80 FRAMEA,FRAMEB,FRAMEC
      CHARACTER    IDENTA*72,IDENTB*72,CUNITA*64,CUNITB*64
      CHARACTER    CBUF*100
C
      DOUBLE PRECISION STARTA(3),STARTB(3),STEPA(3),STEPB(3)
C
      REAL       THR
C
      COMMON  /VMR/  MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      DATA    NPIXA   /1,1,1/,     NPIXB  /1,1,1/
      DATA    STARTA  /3*0.D0/,    STEPA  /3*1.D0/
      DATA    STARTB  /3*0.D0/,    STEPB  /3*1.D0/
      DATA    IDENTA  /' '/,       CUNITA /' '/
      DATA    IDENTB  /' '/,       CUNITB /' '/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT)
      CALL STKRDC('IN_B',1,1,80,IAV,FRAMEB,UNI,NULO,STAT)
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT)
      CALL STKRDR('THR',1,1,IAV,THR,UNI,NULO,STAT)
C
C  map input + reference
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            3,NAXISA,NPIXA,STARTA,STEPA,
     +            IDENTA,CUNITA,PNTRA,IMNOA,STAT)
      IF (FRAMEA.NE.FRAMEB) THEN
         CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +               3,NAXISB,NPIXB,STARTB,STEPB,
     +               IDENTB,CUNITB,PNTRB,IMNOB,STAT)
      ELSE
         IMNOB = IMNOA
         PNTRB = PNTRA
         NAXISB = NAXISA
         DO 20 N=1,NAXISA
            NPIXB(N) = NPIXA(N)
20       CONTINUE
      ENDIF
C
C  get total size + map intermediate frame of that size
      SIZE = 1
      DO 200 N=1,NAXISA
         SIZE = SIZE * NPIXA(N)
200   CONTINUE
      TRUESZ = 1
      DO 300 N=1,NAXISB
         TRUESZ = TRUESZ * NPIXB(N)
300   CONTINUE
      IF (TRUESZ.LT.SIZE) CALL STETER
     +    (8,'refframe smaller than input frame - we abort ...')
C 
      CALL STFXMP(SIZE,D_R4_FORMAT,PNTRC,STAT)
C
C  now do the actual operation
      CALL GXDOIT(MADRID(PNTRA),MADRID(PNTRB),MADRID(PNTRC),
     +          THR,SIZE,TRUESZ)
      IF (TRUESZ.LE.0)
     +   CALL STETER(1,'result frame empty...')
C
      IF (IMNOA.NE.IMNOB)
     +   CALL STFCLO(IMNOB,STAT)            !we don't need the refframe anymore
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            1,TRUESZ,STARTA,STEPA,
     +            IDENTA,CUNITA,PNTRB,IMNOB,STAT)
C
      CALL COPYF(MADRID(PNTRC),MADRID(PNTRB),TRUESZ)
C
C  update descriptors
      CALL DSCUPT(IMNOA,IMNOB,' ',STAT)
      N = INDEX(FRAMEC//' ',' ') - 1
      WRITE(CBUF,10000) FRAMEC(1:N),TRUESZ
      CALL STTPUT(CBUF,STAT)
C
C  that's it folks
      RETURN
C
10000 FORMAT('frame ',A,' created with',I5,' pixels ...')
      END

      SUBROUTINE SUBINS
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine SUBINS		version 1.00	881028
C  K. Banse			ESO - Garching
C
C.KEYWORDS
C  bulk data frames
C
C.PURPOSE
C  insert a frame into another one
C
C.ALGORITHM
C  straight forward
C
C.INPUT/OUTPUT
C  the following keys are used:
C
C  IN_A/C/1/80			input frame
C  OUT_A/C/1/80			result frame
C  P3/C/1/80			string with start pixels
C
C.VERSIONS
C  1.00		from version 1.20 as of 860109
C  		move to FORTRAN 77 + new ST Interfaces
C	
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      IAV,LSUB,N,NAXISA,NAXISC
      INTEGER*8    PNTRX,PNTRY
      INTEGER      NPIXA(3),NPIXC(3),IMNOA,IMNOC,STAT
      INTEGER      SUBLO(3),IAX,NN
      INTEGER      UNI(1),NULO,MADRID(1)
      INTEGER      NPXY,K2,RMAIND,CHUCHO,CHUCHI,IDIM(3)
      INTEGER      NOLIN,TRUSIZ,CHUNKI,CHUNKO,K1
      INTEGER      LOFF,KOFF,MOFF,ZOFF,LOOPLM,NPA(3),PXOFA(3)
C
      CHARACTER*80 FRAMEA,FRAMEC
      CHARACTER    PARM*80,OUTPUT*80
      CHARACTER    SUBSTR*30,ERROR3*40
C
      DOUBLE PRECISION STARTA(3),STARTC(3),STEPA(3),STEPC(3)
      DOUBLE PRECISION EPS,STA(3)
C 
      REAL         CUTS(4),CUTVLS(2),RBUF
C
      COMMON   /VMR/ MADRID
C
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C
      DATA      ERROR3   /'max. 3-dim images supported... '/
      DATA      SUBLO   /1,1,1/,   IDIM   /1,1,1/
      DATA      PARM    /' '/     
      DATA      NPIXA   /1,1,1/,   NPIXC  /1,1,1/
      DATA      NPA     /1,1,1/,   PXOFA  /0,0,0/
      DATA      CUTVLS  /+99999.,-99999./
      DATA      STARTA  /3*0.0D0/, STEPA  /3*1.0D0/
C
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C
C  open son + father frame          (or daughter + mother if you prefer...)
      CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT)
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT)
C
      CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT)     !son
      CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXISA,UNI,NULO,STAT)
      IF (NAXISA.GT.3) CALL STETER(44,ERROR3)
      CALL STDRDI(IMNOA,'NPIX',1,NAXISA,IAV,NPIXA,UNI,NULO,STAT)
      CALL STDRDD(IMNOA,'START',1,NAXISA,IAV,STARTA,UNI,NULO,STAT)
      CALL STDRDD(IMNOA,'STEP',1,NAXISA,IAV,STEPA,UNI,NULO,STAT)
C 
      CALL STFOPN(FRAMEC,D_R4_FORMAT,0,F_IMA_TYPE,IMNOC,STAT)     !father
      CALL STDRDI(IMNOC,'NAXIS',1,1,IAV,NAXISC,UNI,NULO,STAT)
      IF (NAXISC.GT.3) CALL STETER(44,ERROR3)
      CALL STDRDI(IMNOC,'NPIX',1,NAXISC,IAV,NPIXC,UNI,NULO,STAT)
      CALL STDRDD(IMNOC,'START',1,NAXISC,IAV,STARTC,UNI,NULO,STAT)
      CALL STDRDD(IMNOC,'STEP',1,NAXISC,IAV,STEPC,UNI,NULO,STAT)
C 
C  if son has less dimensions, copy stuff from father
      IF (NAXISA.LT.NAXISC) THEN
         DO 300, N=NAXISA+1,NAXISC
            NPIXA(N) = 1
            STARTA(N) = STARTC(N)
            STEPA(N) = STEPC(N)
300      CONTINUE
         NAXISA = NAXISC
      ENDIF
C
C  stepsizes of son should not differ from father
      DO 400, N=1,NAXISA
         EPS = ABS(0.001 * STEPA(N))
         IF (ABS(STEPA(N)-STEPC(N)).GT.EPS) THEN
            CALL STTPUT
     +      ('Warning: stepsizes do not match...',STAT)
            GOTO 1000
         ENDIF
400   CONTINUE
C
C  get start of insertion in FRAMEC
1000  DO 1020 N=1,NAXISC
         STA(N) = STARTC(N)
1020  CONTINUE
C 
      CALL STKRDC('P3',1,1,80,IAV,PARM,UNI,NULO,STAT)
      IF (PARM(1:1).NE.'?') THEN
         IAX = 1
         DO 1100 N=1,NAXISC
            CALL EXTRSS(PARM,',',IAX,SUBSTR,LSUB)
            IF (LSUB.LE.0) GOTO 1200
            IF (SUBSTR(1:1).EQ.'@') THEN
               CALL GENCNV(SUBSTR(2:),1,1,SUBLO(N),RBUF,STA(N),IAV)
               IF (IAV.LE.0)
     +            CALL STETER(46,'Invalid syntax in start coords. ')
               STA(N) = STARTC(N) + ((SUBLO(N)-1)*STEPC(N))
            ELSE IF (SUBSTR(1:1).EQ.'<') THEN
               STA(N) = STARTC(N)
            ELSE IF (SUBSTR(1:1).EQ.'>') THEN
               STA(N) = STARTC(N) + ((NPIXC(N)-1)*STEPC(N))
            ELSE
               CALL GENCNV(SUBSTR(1:),4,1,IAV,RBUF,STA(N),IAV)
               IF (IAV.LE.0)
     +            CALL STETER(46,'Invalid syntax in start coords. ')
            ENDIF
1100     CONTINUE
C
C  if no start values given, use descriptor START of FRAMEA
      ELSE
         DO 1150 N=1,NAXISA                   !note, NAXISA = NAXISC
            STA(N) = STARTA(N)
1150     CONTINUE
      ENDIF
C 
C  find offsets inside FRAMEA
1200  DO 1300 N=1,NAXISA
         IF (ABS(STEPC(N)).LT.1.E-20) THEN
            WRITE(OUTPUT,10123) N,STEPC(N)
            CALL STTPUT(OUTPUT,STAT)
            SUBLO(N) = 0
         ELSE
            SUBLO(N) = NINT( (STA(N) - STARTC(N))/STEPC(N) ) 
         ENDIF
         IF (SUBLO(N).GE.0) THEN
            SUBLO(N) = SUBLO(N) + 1
         ELSE
            PXOFA(N) = - SUBLO(N)
            SUBLO(N) = 1
         ENDIF
         NPA(N) = NPIXA(N) - PXOFA(N)
1300  CONTINUE
C 
C  check, if subframe overlaps with result frame
      DO 2500 N=1,NAXISA
         IF (SUBLO(N).GT.NPIXC(N))
     +      CALL STETER(3,
     +      'invalid coords.: sub- and modframe do not overlap...')
         NN = SUBLO(N) + NPA(N) - 1          !check for overflow
         IF ( NN.GT.NPIXC(N) ) THEN
            IDIM(N) = NPIXC(N) - SUBLO(N) + 1
            CALL STTPUT
     +      ('Warning: son does not fit completely in father...',STAT)
         ELSE
            IDIM(N) = NPA(N)
         ENDIF
2500  CONTINUE
C 
C  get buffersize of at least one line but no more than total size in y
      TRUSIZ = NPIXA(1)*IDIM(2)  
      IAX = 20000                              !first guess ...
      NN = IAX   
      IF (NN.LT.NPIXA(1)) NN = NPIXA(1)
      IF (NN.GT.TRUSIZ) NN = TRUSIZ
      NOLIN = NN / NPIXA(1)                    !buffer = NOLIN * NPIXA(1)
      IF (NOLIN.GT.IDIM(2)) NOLIN = IDIM(2)    !avoid overflow in y ...
C 
3500  CHUNKI = NOLIN * NPIXA(1)
      CHUNKO = NOLIN * NPIXC(1)
      IF ((NOLIN.GT.1).AND.(CHUNKO.GT.IAX*2)) THEN
         NOLIN = NOLIN - 1
         GOTO 3500
      ENDIF
      LOOPLM = IDIM(2) / NOLIN                 !get no. of passes through loop
      RMAIND = IDIM(2) - (LOOPLM * NOLIN)
C 
      CALL STFXMP(CHUNKI,D_R4_FORMAT,PNTRX,STAT)
      CALL STFXMP(CHUNKO,D_R4_FORMAT,PNTRY,STAT)
C 
C  now we loop through input (son) frame
      NPXY = NPIXC(1) * NPIXC(2)
      K2 = (SUBLO(3)-1) * NPXY                   !skip unused planes
      K1 = (SUBLO(2)-1) * NPIXC(1) + 1           !point to 1. used y-lines
      ZOFF = PXOFA(3)*NPIXA(2)*NPIXA(1)
      LOFF = PXOFA(1) + 1                        !first pixel in x-line
      TRUSIZ = NPIXA(1) - LOFF + 1               !no. of pixels per line to move
C 
      DO 5000 IAX=1,IDIM(3)                      !loop over all planes
         MOFF = ZOFF + (PXOFA(2)*NPIXA(1)) + 1   !pointer for input frame
         KOFF = K2 + K1
         IF (LOOPLM.LE.0) GOTO 4400
C 
         DO 4000 N=1,LOOPLM
            CALL STFGET(IMNOA,MOFF,CHUNKI,NN,MADRID(PNTRX),STAT)     !son in
            CALL STFGET(IMNOC,KOFF,CHUNKO,IAV,MADRID(PNTRY),STAT)    !father in
            CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),LOFF,SUBLO(1),
     +              TRUSIZ,NOLIN,NPIXA(1),NPIXC(1),CUTVLS)      !rearrange data
            CALL STFPUT(IMNOC,KOFF,IAV,MADRID(PNTRY),STAT)           !father out
            MOFF = MOFF + CHUNKI
            KOFF = KOFF + CHUNKO                     !follow with indices...
4000     CONTINUE
C 
4400     IF (RMAIND.GT.0) THEN
            CHUCHI = RMAIND * NPIXA(1)
            CHUCHO = RMAIND * NPIXC(1)
            CALL STFGET(IMNOA,MOFF,CHUCHI,NN,MADRID(PNTRX),STAT)
            CALL STFGET(IMNOC,KOFF,CHUCHO,IAV,MADRID(PNTRY),STAT)  
            CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),LOFF,SUBLO(1),
     +               TRUSIZ,RMAIND,NPIXA(1),NPIXC(1),CUTVLS)                
            CALL STFPUT(IMNOC,KOFF,IAV,MADRID(PNTRY),STAT)     
            MOFF = MOFF + CHUCHI
         ENDIF
C 
         K2 = K2 + NPXY                                !move to next plane
         ZOFF = ZOFF + (NPIXA(1)*NPIXA(2))
5000  CONTINUE
C
C  do not copy all descriptors, just update descr HISTORY
      CALL DSCUPT(IMNOC,IMNOC,' ',STAT)
C 
C  and LHCUTS also
      CALL STDRDR(IMNOC,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT)
      IAX = 0
      IF (CUTVLS(1).LT.CUTS(3)) THEN
         CUTS(3) = CUTVLS(1)
         IAX = 1
      ENDIF
      IF (CUTVLS(2).GT.CUTS(4)) THEN
         CUTS(4) = CUTVLS(2)
         IAX = 1
      ENDIF
      IF (IAX.EQ.1) CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT)
C
C  that's it folks
      RETURN
C 
10123 FORMAT('Warning: step(',I1,') of father frame = ',G10.5,
     +       ' result may be corrupted ...')
C 
      END

      SUBROUTINE SUBMAP
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine SUBMAP		version 1.00	881028
C  K. Banse                  	ESO - Garching
C
C.KEYWORDS
C  bulk data frame, mapping
C
C.PURPOSE
C  map an image via another one (which acts like a LUT)
C
C.ALGORITHM
C  straight forward
C
C.INPUT/OUTPUT
C  the following keywords are used:
C
C  IN_A/C/1/80				name of input frame
C  IN_B/C/1/80				name of map frame
C  OUT_A/C/1/80				name of result frame
C  INPUTI/I/1/2				map flag:
C					(1) = 0, indices outside map space go to
C					         start + end
C					(1) = 1, indices outside avoid mapping,
C					         i.e output pixel = input pixel
C					(2) = 0, map frame is 1-dim and has
C						 equidistant coord space
C					(2) = 1, map frame is 2-dim and has
C						 pixel coords stored in 1. line,
C						 pixel intensities in 2. line
C
C.VERSIONS
C
C  1.00		from version 1.00 s of 871112	
C 		move to FORTRAN 77 + new ST Interfaces
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      N,IAV,NAXISA,NAXISB
      INTEGER*8    PNTRA,PNTRB,PNTRC
      INTEGER      IMNOA,IMNOB,IMNOC,STAT
      INTEGER      NPIXA(3),NPIXB(3),KMAP(2)
      INTEGER      UNI(1),NULO,MADRID(1)
C
      DOUBLE PRECISION STARTA(3),STEPA(3),STARTB(3),STEPB(3)
C
      CHARACTER    CUNITA*64,IDENTA*72,CUNITB*10,IDENTB*10
      CHARACTER    FRAMEA*80,FRAMEB*80,FRAMEC*80
C
      COMMON  /VMR/  MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA    IDENTA  /' '/,  CUNITA  /' '/
C
C  get name of input image, map image + output image
      CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT)
      CALL STKRDC('IN_B',1,1,80,IAV,FRAMEB,UNI,NULO,STAT)
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT)
      CALL STKRDI('INPUTI',1,2,IAV,KMAP,UNI,NULO,STAT)
C
      DO 400 N=1,3
         NPIXA(N) = 1
         STARTA(N) = 0.D0
         STEPA(N) = 1.D0
         NPIXB(N) = 1
         STARTB(N) = 0.D0
         STEPB(N) = 1.D0
400   CONTINUE
C
C  map input frame + map frame for reading
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,PNTRA,IMNOA,STAT)
      CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXISB,NPIXB,STARTB,STEPB,IDENTB,
     +            CUNITB,PNTRB,IMNOB,STAT)
      IF (KMAP(2).GT.0) THEN
         IF (NPIXB(2).LE.1)
     +      CALL STETER(1,'we need 2-dim map for this option...')
      ENDIF
C
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,PNTRC,IMNOC,STAT)
C
C  execute algorithm in a subroutine
      CALL GXIZAMAP(NPIXA,NPIXB,STARTB,STEPB,
     +            MADRID(PNTRA),MADRID(PNTRB),MADRID(PNTRC),KMAP)
C
C  that's it folks...
      RETURN
      END

      SUBROUTINE SUBMAT
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  subroutine SUBMAT		version 1.00	881028
C  K. Banse                  	ESO - Garching
C
C.KEYWORDS
C  bulk data frame, matrix ops
C
C.PURPOSE
C  do the usual matrix operations on an image
C
C.ALGORITHM
C  straight forward
C
C.INPUT/OUTPUT
C  the following keywords are used:
C
C  IN_A/C/1/80				name of input frame
C  ACTION/C/3/2				action
C					= 'TR ', transpose a matrix
C 					         around major diagonal
C					= 'LC ', transpose a matrix
C 					         around minor diagonal
C					         i.e. change lines into columns
C  OUT_A/C/1/80				name of result frame
C  INPUTI/I/1/2				xsize,ysize for 'LC'
C
C.VERSIONS
C  1.00		from version 1.20 as of 870105
C	
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER*8    PNTRA,PNTRC
      INTEGER      IAV,NAXIS,STAT
      INTEGER      IMNOA,IMNOC
      INTEGER      NPIX(2),NPIXC(2),ISIZE(2)
      INTEGER      UNI(1),NULO,MADRID(1)
C
      DOUBLE PRECISION START(2),STEP(2),DTMP
C
      CHARACTER    CUNIT*48,IDENT*72
      CHARACTER    FRAMEA*80,FRAMEC*80,ACTION*2
C
      COMMON  /VMR/  MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA    IDENT  /' '/,  CUNIT  /' '/
C
C  get name of input, ouput image + action flags
      CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT)
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT)
      CALL STKRDC('ACTION',1,3,2,IAV,ACTION,UNI,NULO,STAT)
      CALL UPCAS(ACTION,ACTION)
C
C  map input frame
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXIS,NPIX,START,STEP,IDENT,
     +            CUNIT,PNTRA,IMNOA,STAT)
C
C  check dimensions...
      IF (NAXIS.NE.2)
     +   CALL STETER(1,'input must be a 2-dim frame ...')
C
C  exchange NPIX, adjust START, STEP + map output frame
      NPIXC(1) = NPIX(2)
      NPIXC(2) = NPIX(1)
      IF (ACTION.EQ.'TR') THEN
         DTMP = START(1) + (NPIX(1)-1) * STEP(1)          !end in X
         START(1) = START(2) + (NPIX(2)-1) * STEP(2)      !end in Y
         START(2) = DTMP
         DTMP = - STEP(1)
         STEP(1) = - STEP(2)
         STEP(2) = DTMP
      ELSE
         DTMP = START(1)
         START(1) = START(2)
         START(2) = DTMP
         DTMP = STEP(1)
         STEP(1) = STEP(2)
         STEP(2) = DTMP
      ENDIF
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXIS,NPIXC,START,STEP,IDENT,
     +            CUNIT,PNTRC,IMNOC,STAT)
C
C  execute algorithm in a subroutine
      IF (ACTION.EQ.'TR') THEN
         CALL GXMATRX(MADRID(PNTRA),NPIX,MADRID(PNTRC))
      ELSE
         CALL STKRDI('INPUTI',1,2,IAV,ISIZE,UNI,NULO,STAT)
         CALL LINCOL(MADRID(PNTRA),NPIX,ISIZE,MADRID(PNTRC))
      ENDIF
      CALL DSCUPT(IMNOA,IMNOC,' ',STAT)
C
C  that's it folks...
      RETURN
      END
