C @(#)perspec.for	13.1.1.1 (ES0-DMD) 06/02/98 18:20:55
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
      PROGRAM PERSP
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program PERSP 			version 1.00    920103
C  K. Banse                     	ESO - Garching
C  1.10  930512		2.00	940104
C
C.KEYWORDS
C  perspective, rebinning
C
C.INPUT/OUTPUT
C  the following keywords are used:
C
C  IN_A/C/1/60			input frame
C  OUT_A/C/1/60			output frame
C  INPUTR/1/2			rotation angle in degrees, in [0,360]
C				tilt angle in degrees in [0,90]
C
C.VERSIONS
C  1.00   creation
C  1.10   STKWRR had 1 param. too many !
C  2.00   add z-perspective plot
C	
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      IAV,N,NAXISA,NAXISC
      INTEGER*8    PNTRA,PNTRB,PNTRC
      INTEGER      IMNOA,IMNOB,IMNOC,STAT
      INTEGER      NPIXA(3),NPIXB(2),NPIXC(2),KZ,KNO
      INTEGER      N1,N2,NSIZ,NPERSP,DIM3C,KDIFF,PLANES(20)
      INTEGER      UNI(1),NULO,MADRID(1)
      INTEGER      ISCAL
C
      CHARACTER*60 FRAMEA,FRAMEB
      CHARACTER    CUNITA*64,IDENTA*72,STR*80
C
      DOUBLE PRECISION STEPA(3),STARTA(3),STARTB(2),STEPB(2),STEPC(2)
      DOUBLE PRECISION ST(2)
C
      REAL         RNULL
      REAL         ANGDEG,ANGLE,AC,AS,X1,X2,X3,Y1,Y2,Y3,X4,Y4
      REAL         PC, PCDEG
      REAL         PIXROT(2),ROTOFF(2),EN(2),RST(2)
      REAL         QUOT,RQUOT,R1,R2,RROT(4),CUTS(4)
      REAL         FACTO,RBUF(20),ZCUTS(2)
C
      COMMON   /VMR/  MADRID
C
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C 
      EQUIVALENCE (X2,RBUF(1)),(Y2,RBUF(2)),(X4,RBUF(3)),(Y4,RBUF(4))
      EQUIVALENCE (X3,RBUF(5)),(Y3,RBUF(6)),(X1,RBUF(7)),(Y1,RBUF(8))
C
      DATA     FACTO  /0.0174533/               !  Pi / 180.
      DATA     NPIXA  /1,1,1/
      DATA     IDENTA /' '/, CUNITA /' '/       !necessary ... !!!
C
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C 
C  get into MIDAS
      CALL STSPRO('PERSP')
C
C  get input frame + map it
      CALL STKRDC('IN_A',1,1,60,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)
      IF (NAXISA.EQ.1) 
     +   CALL STETER(1,'We need a 2-dim image ...')
      IF ((NPIXA(1).EQ.1) .OR. (NPIXA(2).EQ.1))
     +   CALL STETER(1,'We need a 2-dim image ...')
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT)
      ZCUTS(1) = CUTS(3)
      ZCUTS(2) = CUTS(4)
      RNULL = CUTS(3) - 0.01*(CUTS(4)-CUTS(3))                !null val < min
C 
      R1 = NPIXA(1) - 1.
      R2 = NPIXA(2) - 1.
      KDIFF = 12
      DIM3C = 1
C
C  get name of output frame + plane no.s (or action code)
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEB,UNI,NULO,STAT)
      CALL STKRDC('INPUTC',1,1,80,IAV,STR,UNI,NULO,STAT)
C
C  get rotation angle + perspective angle + (scaling size)
      CALL STKRDR('INPUTR',1,4,IAV,RROT,UNI,NULO,STAT)
      IF (STR(1:4).EQ.'zper') THEN
         PLANES(1) = 1
         PLANES(2) = 0
         KNO = 0
         ISCAL = NINT(RROT(3))
         IF (ISCAL.LT.1) 
     +      CALL STETER(1,'Invalid scaling size (should be > 0)...')
         ZCUTS(2) = RROT(4)
         IF (ZCUTS(2).LE.ZCUTS(1)) CALL STETER(1,
     +      'Invalid scale limit (should be > LHCUTS(3))...')
C 
C  get planes to work on
      ELSE IF ((STR(1:1).EQ.'A') .OR. (STR(1:1).EQ.'a')) THEN
         PLANES(1) = -1                    !indicate all planes
         KNO = NPIXA(3) - 1                !no. of planes - 1
      ELSE
         CALL GENCNV(STR,1,20,PLANES,QUOT,STEPC,IAV)
         IF (IAV.LE.0)
     +      CALL STETER(5,'Invalid planes specified ...')
         IF ((PLANES(1).LE.0).OR.(PLANES(1).GT.NPIXA(3)))
     +      CALL STETER(6,'specified plane outside limits... ')
         N = PLANES(1) - 1                 !update pointer
         KNO = IAV - 1                     !no. of planes - 1
         PNTRA = PNTRA + N*NPIXA(1)*NPIXA(2)
         IAV = IAV + 1
         PLANES(IAV) = 0                   !mark the end
      ENDIF
C 
      ANGDEG = RROT(1)
      PCDEG = RROT(2)
C 
      ANGLE = ANGDEG * FACTO
      N = (NPIXA(1) + 1) / 2                  !use center pixel
      PIXROT(1) = N - 1
      N = (NPIXA(2) + 1) / 2                  !use center pixel
      PIXROT(2) = N - 1
C 
      AC = COS(ANGLE)
      AS = SIN(ANGLE)
      IF ((ANGDEG.LT.0.1).OR.(ANGDEG.GT.359.9)) THEN
         NPIXB(1) = NPIXA(1)                   !no rotation
         NPIXB(2) = NPIXA(2)
         STARTB(1) = STARTA(1)
         STARTB(2) = STARTA(2)
         STEPB(1) = STEPA(1)
         STEPB(2) = STEPA(2)
         N1 = 0
C 
         RBUF(1) = 0.
         RBUF(2) = R2
         RBUF(3) = 0.
         RBUF(4) = 0.
         RBUF(5) = R1
         RBUF(6) = 0.
         RBUF(7) = R1
         RBUF(8) = R2
C 
         GOTO 3000
      ELSE
         N1 = 1                                !yes. we do
      ENDIF
C
C  to get size of output frame, rotate corners around lower left corner
C  (assumed = 0,0)
C 
      CALL ROTA1(R1,R2,AC,AS,X1,Y1)         !(NX-1,NY-1)
      CALL ROTA1(0.,R2,AC,AS,X2,Y2)         !(0,NY-1)
      CALL ROTA1(R1,0.,AC,AS,X3,Y3)         !(NX-1,0)
      X4 = 0.
      Y4 = 0.                               !rotation point stays fixed
C
C  get result frame "corners"
      ST(1) = MIN(0.,X1,X2,X3)             !rotated x-start
      EN(1) = MAX(0.,X1,X2,X3)             !rotated x-end
      ST(2) = MIN(0.,Y1,Y2,Y3)             !rotated y-start
      EN(2) = MAX(0.,Y1,Y2,Y3)             !rotated y-end
C
C  from coords. of rotation point determine shift after rotation
      QUOT = 1. - AC
      ROTOFF(1) = PIXROT(1)*QUOT + PIXROT(2)*AS
      ROTOFF(2) = PIXROT(2)*QUOT - PIXROT(1)*AS
      CALL ROTA1(PIXROT(1),PIXROT(2),AC,AS,RROT(1),RROT(2))
C
C  choose corners such, that rotation point is hit ...
      DO 2000, N=1,2
         STEPB(N) = STEPA(N)
         QUOT =  RROT(N) - ST(N)           !distance of fixed p. to rotation p.
         RQUOT = FLOAT(INT(QUOT)) + 1.     !is forced to multiples of steps
         ST(N) = RROT(N) - RQUOT
         QUOT = EN(N) - ST(N)              !distance of start to end point
         RQUOT = FLOAT(INT(QUOT)) + 1.     !in multiples of steps
         EN(N) = ST(N) + RQUOT
         NPIXB(N) = NINT(RQUOT) + 1        !no. of points in new image
2000  CONTINUE
C
C  now update start of rotated image, so that rotation point stays fixed
      RST(1) = ST(1) + ROTOFF(1)
      RST(2) = ST(2) + ROTOFF(2)
      STARTB(1) = STARTA(1) + ( RST(1)*STEPA(1) )      !we work in 0,0 space
      STARTB(2) = STARTA(2) + ( RST(2)*STEPA(2) )
      DO 2200, N=1,7,2                       
         RBUF(N) = RBUF(N) - ST(1)              !relate to new frame
         RBUF(N+1) = RBUF(N+1) - ST(2)
2200  CONTINUE
C
C  create and map intermediate frame
C 
3000  NSIZ = NPIXB(1) * NPIXB(2)
C 
C  do the rotation
      IF (N1.EQ.1) THEN
         IF (DIM3C.EQ.1) THEN
            CALL STFCRE('midpersp',D_R4_FORMAT,F_X_MODE,1,NSIZ,
     +                  IMNOB,STAT)
            CALL STFMAP(IMNOB,F_X_MODE,1,NSIZ,IAV,PNTRB,STAT)
         ENDIF
         CALL ROTA2(MADRID(PNTRA),MADRID(PNTRB),NPIXA,NPIXB,
     +              ST,AC,AS,RNULL)
      ENDIF
C 
C  create result frame
      IF (DIM3C.EQ.1) THEN
         NPIXC(1) = NPIXB(1)
         STEPC(1) = STEPB(1)
         IF (PCDEG.LE.0.1) THEN                 !no tilt at all ...
            N2 = 0
            PC = 1
            NPIXC(2) = NPIXB(2)
         ELSE IF (PCDEG.GT.89.5) THEN           !total side view => single line
            N2 = 2
            NPIXC(2) = 1
            PC = 0
         ELSE
            N2 = 1
            PC = COS(PCDEG * FACTO)
            NPIXC(2) = INT (NPIXB(2) * PC)
            NPERSP = NPIXC(2)
            STEPC(2) = STEPB(2) * PC
         ENDIF
C
C  move corners to world coords.
         DO 4100, N=1,7,2
            RBUF(N) = STARTB(1) + (RBUF(N)*STEPC(1))
            IAV = INT (RBUF(N+1)*PC)
            IF (IAV.GE.NPIXC(2)) IAV = NPIXC(2) - 1
            RBUF(N+1) = STARTB(2) + (IAV*STEPC(2))
4100     CONTINUE
         IF (KNO.GT.0) THEN                      !if more than one plane
            QUOT = STEPC(2) * ((NPERSP+KDIFF)*KNO)
            RBUF(9) = RBUF(2) + QUOT
            RBUF(10) = RBUF(4) + QUOT
            RBUF(11) = RBUF(6) + QUOT
            RBUF(12) = RBUF(8) + QUOT
            RBUF(13) = KNO + 1                   !no. of planes
            RBUF(14) = KDIFF
            RBUF(15) = STEPC(2) * KDIFF
            RBUF(16) = STEPC(2) * (NPERSP+KDIFF)
C 
            NAXISC = 2
            NPIXC(2) = (KNO+1)*NPIXC(2) + KNO*KDIFF
            RQUOT = STARTB(2) + (NPIXC(2)-1)*STEPC(2)         !y-end
            IF (STEPC(2).GT.0.) THEN
               IF (RBUF(9).GT.RQUOT) RBUF(9) = RQUOT
               IF (RBUF(10).GT.RQUOT) RBUF(10) = RQUOT
               IF (RBUF(11).GT.RQUOT) RBUF(11) = RQUOT
               IF (RBUF(12).GT.RQUOT) RBUF(12) = RQUOT
            ELSE
               IF (RBUF(9).LT.RQUOT) RBUF(9) = RQUOT
               IF (RBUF(10).LT.RQUOT) RBUF(10) = RQUOT
               IF (RBUF(11).LT.RQUOT) RBUF(11) = RQUOT
               IF (RBUF(12).LT.RQUOT) RBUF(12) = RQUOT
            ENDIF
         ELSE
            NAXISC = NAXISA
         ENDIF
C 
         CALL STKWRR('OUTPUTR',RBUF,1,20,UNI,STAT)
         IF (STR(1:4).EQ.'zper') THEN
            IDENTA(1:) = 'perspective z-plot of input frame '
            NAXISC = 2
            NPIXC(2) = NPIXC(2)+ISCAL
            NPERSP = NPIXC(2)                     !in case we reset it later on
         ELSE
            IDENTA(1:) = 
     +      'perspective plot of selected planes of input frame '
         ENDIF
         CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXISC,NPIXC,STARTB,STEPC,
     +               IDENTA,CUNITA,PNTRC,IMNOC,STAT)
      ENDIF
C
C  and do the perspective rebinning
      IF (N2.EQ.1) THEN                         !real job
         NPIXC(2) = NPERSP                      !reset to original value (3-dim)
         IF (N1.EQ.1) THEN
            IF (STR(1:4).EQ.'zper') THEN
               CALL JPERSP(MADRID(PNTRB),MADRID(PNTRC),
     +                     NPIXB,NPIXC,PC,RNULL,ZCUTS,ISCAL)
               CALL STDWRI(IMNOC,'NPIX',NPIXC(2),2,1,UNI,STAT)
            ELSE
               CALL KPERSP(MADRID(PNTRB),MADRID(PNTRC),
     +                     NPIXB,NPIXC,PC)
            ENDIF
         ELSE
            IF (STR(1:4).EQ.'zper') THEN
               CALL JPERSP(MADRID(PNTRA),MADRID(PNTRC),
     +                     NPIXB,NPIXC,PC,RNULL,ZCUTS,ISCAL)
               CALL STDWRI(IMNOC,'NPIX',NPIXC(2),2,1,UNI,STAT)
            ELSE
               CALL KPERSP(MADRID(PNTRA),MADRID(PNTRC),
     +                     NPIXB,NPIXC,PC)
            ENDIF
         ENDIF
      ELSE IF (N2.EQ.0) THEN                    !just rotated image
         IF (N1.EQ.1) THEN
            CALL COPYF(MADRID(PNTRB),MADRID(PNTRC),NSIZ)
         ELSE
            CALL COPYF(MADRID(PNTRA),MADRID(PNTRC),NSIZ)
         ENDIF
      ELSE                                      !only sideline
         IF (N1.EQ.1) THEN
            CALL COPYF(MADRID(PNTRB),MADRID(PNTRC),NPIXB(1))
         ELSE
            CALL COPYF(MADRID(PNTRA),MADRID(PNTRC),NPIXB(1))
         ENDIF
      ENDIF
C 
C  check for 3-dim case
      DIM3C = DIM3C + 1
      IF (PLANES(1).EQ.-1) THEN
         IF (DIM3C.GT.NPIXA(3)) THEN
            GOTO 9000
         ELSE
            KZ = 1
         ENDIF
      ELSE
         IF (PLANES(DIM3C).EQ.0) THEN
            GOTO 9000
         ELSE
            IF ((PLANES(DIM3C).LE.0).OR.(PLANES(DIM3C).GT.NPIXA(3)))
     +         CALL STETER(6,'specified plane outside limits... ')
            KZ = PLANES(DIM3C) - PLANES(DIM3C-1)
         ENDIF
      ENDIF
      PNTRA = PNTRA + (NPIXA(1)*NPIXA(2))*KZ
      PNTRC = PNTRC + (NPIXC(1) * (NPIXC(2) + KDIFF))
      GOTO 3000
C
9000  CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT)
      CALL STSEPI
      END

      SUBROUTINE ROTA1(X,Y,AC,AS,XOUT,YOUT)
C
      IMPLICIT NONE
C
      REAL     X,Y,AC,AS,XOUT,YOUT
C
      XOUT = X*AC - Y*AS
      YOUT = X*AS + Y*AC
C
      RETURN
      END

      SUBROUTINE KPERSP(A,B,NPIXA,NPIXB,AC)
C
C  K. Banse            920103
C
C  do "backward" perspective
C 
C  XPERS = X
C  YPERS = Y/AC
C
      IMPLICIT NONE
C
      INTEGER      INDXB,KSTX,KY
      INTEGER      NO1,NO3
      INTEGER      NX,NY,SIZEA,XDIMA,YDIMA
      INTEGER      NPIXA(2),NPIXB(2)
C
      REAL         A(*),B(*)
      REAL         AC,YFRAC
C 
      DOUBLE PRECISION  RYAC,ACINV
C
C  init
      XDIMA = NPIXA(1)
      YDIMA = NPIXA(2)
      SIZEA = XDIMA * YDIMA
C
      ACINV = 1.D0 / AC
      INDXB = 0
      RYAC = ACINV
C
C  loop through result lines
      DO 3000, NY=1,NPIXB(2)
         DO 1000, NX=1,NPIXB(1)
            INDXB = INDXB + 1              !loop along result line
            KY = INT(RYAC)
            KSTX = XDIMA*KY                !get 1. pixel in this line - 1
            NO1 = NX + KSTX                !get pixel close to integer part of
            NO3 = NO1 + XDIMA               !get low pixel in next line up
            IF (NO3.GT.SIZEA) THEN          !are we out of grid in y?
               B(INDXB) = A(NO1)
            ELSE
               YFRAC = RYAC - KY
               B(INDXB) = A(NO1) + (YFRAC * (A(NO3) - A(NO1)))
            ENDIF
1000     CONTINUE
         RYAC = RYAC + ACINV
3000  CONTINUE
C 
      RETURN
      END

      SUBROUTINE JPERSP(A,B,NPIXA,NPIXB,AC,ZNULL,FMINMA,ISCAL)
C
C
C  do "backward" z-perspective
C
C  XPERS = X
C  YPERS = Y/AC 
C 
C  Z = YPERS * ISCAL  =>      XPERS,YPERS+n   = Z
C                             XPERS,YPERS+n-1 = ...
C                                  ...
C                             XPERS,YPERS     = Min
C
      IMPLICIT NONE
C
      INTEGER      INDXB,KSTX,KY,ISCAL
      INTEGER      NO1,NO3,N
      INTEGER      NX,NY,SIZEA,XDIMA,YDIMA,XDIMB
      INTEGER      NPIXA(2),NPIXB(2)
      INTEGER      KSCAL,IZ,KNDXB,KK
C
      REAL         A(*),B(*)
      REAL         AC,FMINMA(2),ZNULL
      REAL         YFRAC,RMIN,RMAX
      REAL         Z,DZ,SCALF,CMIN,FF
      REAL         VALU
C
      DOUBLE PRECISION  RYAC,ACINV
C
C  init
      XDIMA = NPIXA(1)
      YDIMA = NPIXA(2)
      SIZEA = XDIMA * YDIMA
      XDIMB = NPIXB(1)
C
      RMIN = FMINMA(1)
      RMAX = FMINMA(2)
      ACINV = 1.D0 / AC
      INDXB = NPIXB(1)*NPIXB(2)
      KSCAL = ISCAL - 1
      SCALF = KSCAL / (RMAX-RMIN)
      CMIN = RMIN*ISCAL
      KK = NPIXB(2)-ISCAL               !reserved space for spikes of last line
C
C  loop through result lines
      DO 300, NY=NPIXB(2),KK+1,-1
         DO 100, NX=1,XDIMB
            B(INDXB) = ZNULL
            INDXB = INDXB - 1              !loop along result line
100      CONTINUE
300   CONTINUE
C 
      RYAC = KK * ACINV
      DO 3000, NY=KK,1,-1
C
         DO 1000, NX=1,XDIMB
            KY = INT(RYAC)
            KSTX = XDIMA*KY                !get 1. pixel in this line - 1
            NO1 = NX + KSTX                !get pixel close to integer part of
            NO3 = NO1 + XDIMA              !get low pixel in next line up
            IF (NO3.GT.SIZEA) THEN         !are we out of grid in y?
               Z = A(NO1)
            ELSE
               YFRAC = RYAC - KY
               Z = A(NO1) + (YFRAC * (A(NO3) - A(NO1)))
            ENDIF
            IF (Z .GT. RMAX) Z = RMAX
C 
C  get height from intensity
            DZ = Z - RMIN
            IZ = NINT(DZ*SCALF)
            IF (IZ .GE. 1) THEN
               B(INDXB) = RMIN
               FF = DZ/IZ
               KNDXB = INDXB + XDIMB
               VALU = FF + RMIN
               DO 800, N=1,IZ
                  B(KNDXB) = VALU
                  VALU = VALU + FF
                  KNDXB = KNDXB + XDIMB
800            CONTINUE
            ELSE
               B(INDXB) = Z
            ENDIF
C
            INDXB = INDXB - 1              !loop along result line
1000     CONTINUE
C
         RYAC = RYAC - ACINV
3000  CONTINUE
C
C  get rid of superfluous lines in the end
      INDXB = NPIXB(1)*NPIXB(2)
      DO 4000, NY=NPIXB(2),KK+1,-1
         DO 3300, NX=1,XDIMB
            IF (B(INDXB) .GT. ZNULL) THEN
               NPIXB(2) = NY
               RETURN
            ENDIF
            INDXB = INDXB - 1
3300     CONTINUE
4000  CONTINUE
C 
      RETURN
      END

      SUBROUTINE ROTA2(A,B,NPIXA,NPIXB,START,AC,AS,ZNULL)
C
C  K. Banse            861029, 880902, 910204
C
C  do "backward" rotation, i.e. solve the system XROT = X*AC - Y*AS
C                                       YROT = X*AS + Y*AC for X,Y
C  to yield X = XROT*AC + YROT*AS
C          Y = YROT*AC - XROT*AS  which takes AC**2 + AS**2 = 1 into account....
C
      IMPLICIT NONE
C
      INTEGER      INDXB,KX,KY,KSTX
      INTEGER      NO1,NO2,NO3,NO4
      INTEGER      NX,NY,SIZEA,XDIMA,YDIMA
      INTEGER      NPIXA(2),NPIXB(2)
C
      REAL         A(*),B(*)
      REAL         AC,AS,ZNULL
      REAL         XT,YT,XFRAC,YFRAC,XLAST,YLAST
      REAL         XF,RX
C
      DOUBLE PRECISION  START(2)
      DOUBLE PRECISION  RYAC,RYAS
C
C  init
      XF = START(1)
      XDIMA = NPIXA(1)
      YDIMA = NPIXA(2)
      XLAST = XDIMA - 1.
      YLAST = YDIMA - 1.
      SIZEA = XDIMA * YDIMA
      INDXB = 0
C
      RYAS = START(2) * AS                 !init RYAS,RYAC according to start-y
      RYAC = START(2) * AC
C 
C  loop through result lines
      DO 3000, NY=1,NPIXB(2)
         RX = XF                           !init RX to start-x
C
         DO 1000, NX=1,NPIXB(1)
            INDXB = INDXB + 1              !loop along result line
            XT = RX*AC + RYAS
            IF ((XT.LT.0.).OR.(XT.GT.XLAST)) THEN
               B(INDXB) = ZNULL
               GOTO 1000
            ENDIF
            YT = RYAC - RX*AS
            IF ((YT.LT.0.).OR.(YT.GT.YLAST)) THEN
               B(INDXB) = ZNULL
               GOTO 1000
            ENDIF
C
C  interpolate...
            XT = XT + 1.                  !adjust 0,0 to 1,1 ...
            YT = YT + 1.
            KX = INT(XT)
            KY = INT(YT)
            KSTX = XDIMA*(KY-1)          !get 1. pixel in this line - 1
            NO1 = KX + KSTX          !get pixel close to int. part of real pixel
C
            NO2 = NO1 + 1                  !get next pixel in x-direction
CC            IF ((NO2-KSTX).GT.XDIMA) THEN  !are we out of grid in x?
            IF (KX.GE.XDIMA) THEN          !are we out of grid in x?
               IF (NO2.GT.SIZEA) THEN      !also out in y?
                  B(INDXB) = A(NO1)        !yes.
               ELSE
                  NO3 = NO1 + XDIMA        !no.
                  YFRAC = YT - KY
                  B(INDXB) = A(NO1) + ( YFRAC * (A(NO3) - A(NO1)) )
               ENDIF
            ELSE
               XFRAC = XT - KX
               NO3 = NO1 + XDIMA               !get low pixel in next line up
               IF (NO3.GT.SIZEA) THEN          !are we out of grid in y?
                  B(INDXB) = A(NO1) + ( XFRAC * (A(NO2) - A(NO1)) )
               ELSE
                  NO4 = NO3 + 1
                  YFRAC = YT - KY
                  B(INDXB) = A(NO1) +
     +              XFRAC * (A(NO2) - A(NO1)) +
     +              YFRAC * (A(NO3) - A(NO1)) +
     +              XFRAC * YFRAC * (A(NO1) - A(NO2) - A(NO3) + A(NO4))
               ENDIF
            ENDIF
C
1000     RX = RX + 1.0                        !move to next x-coord
C
         RYAS = RYAS + AS                     !move to next y-coord.
         RYAC = RYAC + AC       
3000  CONTINUE
C 
      RETURN
      END