* Last processed by NICE on 12-Jun-2000 15:51:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      PROGRAM REPROJECT
C----------------------------------------------------------------------
C GDF	Stand alone program
C	Resample an input image to a different projection system
C       Now accepts EPOCH conversion for Equatorials. Forced interpola-
C       tion to bilinear.
C----------------------------------------------------------------------
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:zpar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:errcod.inc'
      INCLUDE 'inc:pi.inc'
      INCLUDE 'inc:memory.inc'
*
      CHARACTER*80 NAMEX,NAMEY,NAMEZ,CENTER_1,CENTER_2
      LOGICAL ERROR
      INTEGER NKEY,LENC,SYS_CODE,MVOC
      PARAMETER (MVOC=8)
      CHARACTER*16 VOCAB(MVOC),KEY,PROJ,SYSTEM
      CHARACTER*40 SYSTCHAIN
      CHARACTER*132 CHAIN
*
      INTEGER N, DIM(2), IL, IB
      REAL ARR(2),NEWEPOCH
      REAL*8 CENTER(2),ANGLE,CONV(6)
      LOGICAL CHANGE,TEMPLATE
      INTEGER*4 ADDR1,ADDR2, IPX, IPY, IP1, IP2, POINTER
      INTEGER*4 BLC(4), TRC(4)
*
      DATA BLC /4*0/, TRC/4*0/, CONV/6*0/
      DATA VOCAB/'NONE','GNOMONIC','ORTHOGRAPHIC','AZIMUTHAL',
     $'STEREOGRAPHIC','LAMBERT','AITOFF','RADIO'/
*
      CALL GILDAS_OPEN
      CALL GILDAS_CHAR('Z_NAME$',NAMEZ)
      CALL GILDAS_CHAR('X_NAME$',NAMEX)
      CALL GILDAS_LOGI('TEMPLATE$',TEMPLATE,1)
      IF (TEMPLATE) THEN
         CALL  GILDAS_CHAR('Y_NAME$',NAMEY)
      ELSE
         CALL GILDAS_CHAR('PROJECTION$',PROJ)
         CALL GILDAS_CHAR('SYSTEM$',SYSTCHAIN)
         CALL GILDAS_CHAR('CENTER_1$',CENTER_1)
         CALL GILDAS_CHAR('CENTER_2$',CENTER_2)
         CALL GILDAS_DBLE('ANGLE$',ANGLE,1)
         CALL GILDAS_INTE('DIMENSIONS$',DIM,2)
         CALL GILDAS_DBLE('AXIS_1$',CONV,3)
         CALL GILDAS_DBLE('AXIS_2$',CONV(4),3)
      ENDIF
      CALL GILDAS_LOGI('CHANGE$',CHANGE,1)
      IF (CHANGE) CALL GILDAS_REAL('BLANKING$',ARR,2)
      CALL GILDAS_CLOSE
*
      N = LENC(NAMEZ)
      CALL SIC_PARSEF(NAMEZ(1:N),Z_FILE,' ','.gdf')
      CALL GDF_GEIS (Z_ISLO,ERROR)
      CALL GDF_REIS (Z_ISLO,Z_TYPE,Z_FILE,Z_FORM,Z_SIZE,ERROR)
      IF (ERROR) THEN
         CALL GAGOUT('F-REPROJECT,  Cannot read input file')
         GOTO 100
      ENDIF
      IF (Z_FORM.NE.FMT_R4) THEN
         CALL GAGOUT('E-REPROJECT,  Only real images supported')
         GOTO 100
      ENDIF
      CALL GDF_READZ (Z_ISLO,ERROR)
*
      CALL GDF_CHZY
*
* Now define output
      IF (TEMPLATE) THEN
         N = LENC(NAMEY)
         CALL SIC_PARSEF(NAMEY(1:N),Y_FILE,' ','.gdf')
         CALL GDF_GEIS (Y_ISLO,ERROR)
         CALL GDF_REIS (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
         IF (ERROR) THEN
            CALL GAGOUT('F-REPROJECT,  Cannot read input Template file')
            GOTO 100
         ENDIF
         CALL GDF_READY (Y_ISLO,ERROR)
      ENDIF
      CALL GDF_CHYX
*
      N = LENC(NAMEX)
      CALL SIC_PARSEF(NAMEX(1:N),X_FILE,' ','.gdf')
      IF (.NOT.TEMPLATE) THEN
         IF (PROJ.NE.'UNCHANGED') THEN
            CALL SIC_UPPER(PROJ)
            N = LENC(PROJ)
            IF (N.NE.0) THEN
               CALL SIC_AMBIGS('PROJECTION',PROJ,KEY,NKEY,
     $         VOCAB,MVOC,ERROR)
               IF (ERROR) GOTO 98
               X_PTYP = NKEY-1
            ENDIF
         ENDIF
*
* Coordinate system. test presence of an epoch
*
         NEWEPOCH=Y_EPOC
         IL = LENC(SYSTCHAIN)
         CALL SIC_BLANC(SYSTCHAIN,IL)
         IB = INDEX(SYSTCHAIN,' ')
         IF (IB.NE.0.AND.IL.GT.IB) THEN
            READ(SYSTCHAIN(IB+1:IL),*,ERR=1055) NEWEPOCH
            GOTO 1056
1055        CALL GAGOUT('W-REPROJECT, Error decoding Epoch')
            NEWEPOCH=Y_EPOC
1056        SYSTEM=SYSTCHAIN(1:IB-1)
         ELSE
            SYSTEM=SYSTCHAIN(1:IL)
         ENDIF
         CALL SIC_UPPER(SYSTEM)
         IF (SYSTEM.NE.'UNCHANGED') THEN
            IF (SYSTEM.EQ.'GALACTIC') THEN
               X_SYST = SYSTEM
               IF (Y_SYST.EQ.'UNKNOWN') THEN
                  GOTO 102
               ELSEIF (Y_SYST.EQ.'EQUATORIAL') THEN
                  SYS_CODE = 1
               ELSE
                  SYS_CODE = 0
               ENDIF
            ELSEIF (SYSTEM.EQ.'EQUATORIAL') THEN
               X_SYST = SYSTEM
               X_EPOC = NEWEPOCH
               IF (Y_SYST.EQ.'UNKNOWN') THEN
                  GOTO 102
               ELSEIF (Y_SYST.EQ.'GALACTIC') THEN
                  SYS_CODE = -1
               ELSEIF (Y_EPOC.NE.X_EPOC) THEN
                  WRITE(CHAIN,1002) Y_EPOC,X_EPOC
                  CALL GAGOUT(CHAIN)
                  SYS_CODE = -2
               ELSE
                  SYS_CODE = 0
               ENDIF
            ELSEIF (SYSTEM.EQ.'UNKNOWN') THEN
               X_SYST = SYSTEM
               IF (Y_SYST.EQ.'GALACTIC') THEN
                  GOTO 102
               ELSEIF (SYSTEM.EQ.'EQUATORIAL') THEN
                  GOTO 102
               ENDIF
               SYS_CODE = 0
            ENDIF
         ELSE
            SYS_CODE = 0
         ENDIF
         CALL SIC_UPPER(CENTER_1)
         CALL SIC_UPPER(CENTER_2)
* we could guess center also:
         IF(CENTER_1.EQ.'UNCHANGED'.OR.CENTER_2.EQ.'UNCHANGED')THEN
            IF (X_SYST.EQ.'EQUATORIAL') THEN
               IF(Y_SYST.EQ.'EQUATORIAL') THEN
                  CALL EQU_EQU(Y_A0,Y_D0,Y_EPOC,X_A0,X_D0,X_EPOC,1)
               ELSEIF(Y_SYST.EQ.'GALACTIC')THEN
                  CALL GAL_EQU(Y_A0,Y_D0,X_A0,X_D0,X_EPOC,1)
               ELSE
                  CHAIN='F-REPROJECT, '
     $            //'UNCHANGED Center not permitted here'
                  CALL GAGOUT(CHAIN)
                  GOTO 100
               ENDIF
            ELSEIF(X_SYST.EQ.'GALACTIC')THEN
               IF(Y_SYST.EQ.'EQUATORIAL') THEN
                  CALL EQU_GAL(Y_A0,Y_D0,Y_EPOC,X_A0,X_D0,1)
               ELSEIF(Y_SYST.EQ.'GALACTIC')THEN
* do nothing
               ELSE
                  CHAIN='F-REPROJECT, '
     $            //'UNCHANGED Center not permitted here'
                  CALL GAGOUT(CHAIN)
                  GOTO 100
               ENDIF
            ENDIF
         ELSE
            IF (X_SYST.EQ.'EQUATORIAL') THEN
               CALL SIC_DECODE(CENTER_1,CENTER(1),24,ERROR)
               IF (ERROR) GOTO 101
            ELSE
               CALL SIC_DECODE(CENTER_1,CENTER(1),360,ERROR)
               IF (ERROR) GOTO 101
            ENDIF
            CALL SIC_DECODE(CENTER_2,CENTER(2),360,ERROR)
            IF (ERROR) GOTO 101
            X_A0 = CENTER(1)
            X_D0 = CENTER(2)
         ENDIF
         IF (X_PTYP.EQ.7.OR.X_PTYP.EQ.6 .AND. ANGLE.NE.0.0) THEN
            CHAIN='W-REPROJECT, '
     $      //'Angle ignored in This Projection'
            CALL GAGOUT(CHAIN)
            ANGLE=0.0
         ENDIF
         X_PANG = ANGLE*PI/180.0
*
         IF (DIM(1).EQ.0.OR.DIM(2).EQ.0.OR.CONV(3).EQ.0
     $   .OR.CONV(6).EQ.0) THEN
            CALL GET_VM(2*Z_DIM(1)*Z_DIM(2),ADDR1,*99)
            IP1 = POINTER(ADDR1,MEMORY)
            CALL GET_VM(2*Z_DIM(1)*Z_DIM(2),ADDR2,*99)
            IP2 = POINTER(ADDR2,MEMORY)
            CALL GDF_GEMS (Z_MSLO,Z_ISLO,BLC,TRC,Z_ADDR,Z_FORM,ERROR)
            IPY = POINTER(Z_ADDR,MEMORY)
* "no blanking" should be explicit:
            IF (Z_BLAN.EQ.0) THEN
               Z_BVAL = 0.0
               Z_EVAL = -1.0
            ENDIF
            CALL FINDGRID(
     $      MEMORY(IPY),Z_DIM(1),Z_DIM(2),Z_DIM(3)*Z_DIM(4),
     $      Z_CONV,Z_PTYP,Z_A0,Z_EPOC,Z_BVAL,Z_EVAL,
     $      X_PTYP,X_A0,X_EPOC,SYS_CODE,
     $      MEMORY(IP1),MEMORY(IP2),CONV,DIM)
            CALL FREE_VM (2*Z_DIM(1)*Z_DIM(2),ADDR1)
            CALL FREE_VM (2*Z_DIM(1)*Z_DIM(2),ADDR2)
            CALL GDF_FRMS(Z_ISLO,ERROR)
         ENDIF
         CALL R8TOR8 (CONV,X_CONV,6)
         X_DIM(1) = DIM(1)
         X_DIM(2) = DIM(2)
         X_DIM(3) = MAX(1,X_DIM(3))
         X_DIM(4) = MAX(1,X_DIM(4))
         X_SIZE = X_DIM(1)*X_DIM(2)*X_DIM(3)*X_DIM(4)
*in case everything is not present in the input file:
         X_GENE = 29
         X_BLAN = 2
         X_EXTR = 10
         X_DESC = 18
         X_POSI = 12
         X_PROJ = 9
      ENDIF
      IF (Z_BLAN.EQ.0) THEN
         Z_BVAL = 0.0
         Z_EVAL = -1.0
         X_BLAN = 2
         X_BVAL = ARR(1)
         X_EVAL = ARR(2)
      ELSEIF (CHANGE) THEN
         X_BLAN = 2
         X_BVAL = ARR(1)
         X_EVAL = ARR(2)
      ENDIF
      IF (TEMPLATE)X_FORM=Z_FORM
*
* Method of solution
      CALL GAGOUT('I-REPROJECT,  Using bilinear gridding method')
      IF (Z_EVAL.LT.0.0) THEN
         CALL GAGOUT('              with no input blanking...')
      ELSE
         CALL GAGOUT('              taking input blanking in account')
      ENDIF
*
* Create output image
      ERROR = .FALSE.
      CALL GDF_GEIS (X_ISLO,ERROR)
      CALL GDF_WRITX(X_ISLO,ERROR)
      CALL GDF_CRIS (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
      IF (ERROR) THEN
         CALL GAGOUT('F-REPROJECT,  Cannot create X file:')
         CALL GAGOUT(X_FILE)
         GOTO 100
      ENDIF
*
      CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
      CALL GDF_GEMS (Z_MSLO,Z_ISLO,BLC,TRC,Z_ADDR,Z_FORM,ERROR)
      IPX = POINTER(X_ADDR,MEMORY)
      IPY = POINTER(Z_ADDR,MEMORY)
      CALL GET_VM(2*X_DIM(1)*X_DIM(2),ADDR1,*99)
      IP1 = POINTER(ADDR1,MEMORY)
      CALL GET_VM(2*X_DIM(1)*X_DIM(2),ADDR2,*99)
      IP2 = POINTER(ADDR2,MEMORY)
      CALL GAGOUT ('Start computing')
      CALL GRIDLIN(
     $MEMORY(IPY),Z_DIM(1),Z_DIM(2),Z_DIM(3)*Z_DIM(4),
     $Z_CONV,Z_PTYP,Z_A0,Z_EPOC,Z_BVAL,
     $MEMORY(IPX),X_DIM(1),X_DIM(2),
     $X_CONV,X_PTYP,X_A0,X_EPOC,X_BVAL,
     $SYS_CODE,MEMORY(IP1),MEMORY(IP2))
*
      CALL FREE_VM (2*Z_DIM(1)*Z_DIM(2),ADDR1)
      CALL FREE_VM (2*Z_DIM(1)*Z_DIM(2),ADDR2)
      CALL GDF_FRIS (X_ISLO,ERROR)
      CALL GDF_FRIS (Z_ISLO,ERROR)
      IF (TEMPLATE) CALL GDF_FRIS(Y_ISLO,ERROR)
*
      STOP 'S-REPROJECT,  Successful completion'
98    CALL GAGOUT('F-REPROJECT,  Invalid projection system')
      CALL SYSEXI (FATALE)
99    CALL GAGOUT('F-REPROJECT,  Lost my mind')
100   CALL SYSEXI (FATALE)
101   WRITE(6,1000)
1000  FORMAT('F-REPROJECT,  Bad syntax in center coordinates')
      CALL SYSEXI (FATALE)
102   WRITE(CHAIN,1001) Y_SYST, X_SYST
      CALL GAGOUT(CHAIN)
      CALL SYSEXI (FATALE)
1001  FORMAT ('F-REPROJECT, Cannot convert from ',A,' to ',A)
1002  FORMAT ('Converting from epoch ',F8.2,' to ',F8.2)
      END
*<FF>
      SUBROUTINE FINDGRID(A,MX,MY,MZ,ACONV,ATYPE,APROJ,AEPOC,
     $ABVAL,AEVAL,BTYPE,BPROJ,BEPOC,CODE,DXD,DYD,BCONV,BDIM)
C----------------------------------------------------------------------
C GILDAS	Internal routine
C Arguments
C       A       R*4     pseudo-map
C	MX,MY	I       size of pseudo-map
C	ACONV	R*8(6)	Pixel conversion formulae
C			CONV(1) = Xref, CONV(2)=Xval, CONV(3)=Xinc
C	ATYPE	I	Type of projection
C	APROJ	R*8(3)	Projection constants
C			PROJ(1) = A0,	PROJ(2)=D0,	PROJ(3)=Angle
C       AEPOC   R*4     Epoch if Needed
C       ABVAL   R*4     Blank val
C       AEVAL   R*4     Blank errval
C	BTYPE	I
C	BPROJ	R*8(3)
C       BEPOC   R*4     Epoch if Needed
C	CODE	I	(-1) GA --> EQU	et (1) EQ --> GAL (-2) EQ->EQ
C       DXD     R*8(*)  1-D arrays of projected coordinates (work)
C       DYD     R*8(*)  id.
C	BCONV	R*8(6)  output best conv
C	BDIM	I(2)    output best DIM
C----------------------------------------------------------------------
      INCLUDE 'inc:pi.inc'
      INTEGER MX,MY,MZ,ATYPE,BTYPE,CODE,BDIM(2)
      REAL*4 A(MX,MY,MZ),ABVAL, AEVAL
      REAL*4 AEPOC,BEPOC
      REAL*8 DXD(*),DYD(*)
      REAL*8 ACONV(6), APROJ(6), BCONV(6), BPROJ(6), BOUNDARY
*
      INTEGER I,J,NDP
      REAL*8 AXREF,AXVAL,AXINC, AYREF,AYVAL,AYINC
      REAL*8 DXMIN,DXMAX,DYMIN,DYMAX
      REAL SIZX,SIZY,RATIO
      CHARACTER*80 CHAIN
      LOGICAL WARN
*
* Preliminary processing
      AXREF = ACONV(1)
      AXVAL = ACONV(2)
      AXINC = ACONV(3)
      AYREF = ACONV(4)
      AYVAL = ACONV(5)
      AYINC = ACONV(6)
*
* Load X,Y,Z arrays
*
* Setup Input projection
      CALL GR8_PROJEC (ATYPE,APROJ)
*
* Convert Input projection to Absolute coordinates and load data points
      NDP=0
      DO J=1,MY
         DO I=1,MX
            IF (ABS(A(I,J,1)-ABVAL).GE.AEVAL) THEN
               NDP = NDP+1
               DXD(NDP) = (I-AXREF)*AXINC + AXVAL
               DYD(NDP) = (J-AYREF)*AYINC + AYVAL
            ENDIF
         ENDDO
      ENDDO
      CALL GAGOUT('I-REPROJECT,  Evaluating Output Image Header...')
      WRITE(CHAIN,1000) NDP
      CALL GAGOUT(CHAIN)
      IF (NDP.EQ.0) STOP ' No valid data points in input image'
      CALL REL_TO_ABS(DXD,DYD,DXD,DYD,NDP)
*
* Change coordinate system. use new epoch also
      IF (CODE.EQ.1) THEN
         CALL EQU_GAL(DXD,DYD,AEPOC,DXD,DYD,NDP)
      ELSEIF (CODE.EQ.-1) THEN
         CALL GAL_EQU(DXD,DYD,DXD,DYD,BEPOC,NDP)
      ELSEIF (CODE.EQ.-2) THEN
         CALL EQU_EQU(DXD,DYD,AEPOC,DXD,DYD,BEPOC,NDP)
      ENDIF
*
* Setup Output projection
      CALL GR8_PROJEC (BTYPE,BPROJ)
*
* Convert Absolute coordinates to input projection
      CALL ABS_TO_REL(DXD,DYD,DXD,DYD,NDP)
*
      CALL BESTBOX(DXD,DYD,NDP,DXMIN,DXMAX,DYMIN,DYMAX)
* Test extreme cases of unbound projections:
      IF (BTYPE.EQ.1)THEN                !GNOMONIC
         BOUNDARY=1.5D0
      ELSEIF(BTYPE.EQ.4)THEN             !STEREOGRAPHIC
         BOUNDARY=1.0D0
      ENDIF
      IF (BTYPE.EQ.1.OR.BTYPE.EQ.4)THEN
         IF (DXMIN.LT.-BOUNDARY) THEN
            WARN=.TRUE.
            DXMIN=-BOUNDARY
         ENDIF
         IF (DYMIN.LT.-BOUNDARY) THEN
            WARN=.TRUE.
            DYMIN=-BOUNDARY
         ENDIF
         IF (DXMAX.GT.BOUNDARY) THEN
            WARN=.TRUE.
            DXMAX=BOUNDARY
         ENDIF
         IF (DYMAX.GT.BOUNDARY) THEN
            WARN=.TRUE.
            DYMAX=BOUNDARY
         ENDIF
         IF(WARN) THEN
            CHAIN='W-REPROJECT, Output Projection is unbound'
            CALL GAGOUT(CHAIN)
            CHAIN='Will be limited to one hemisphere'
            CALL GAGOUT(CHAIN)
         ENDIF
      ENDIF
*
      SIZX=DXMAX-DXMIN
      SIZY=DYMAX-DYMIN
      RATIO=SIZX/SIZY
* If user gave number of pixels, use it
      IF (BDIM(2).EQ.0) THEN
         BDIM(2)=NINT(SQRT(FLOAT(MX*MY)*2.0/RATIO))
      ENDIF
      IF (BDIM(1).EQ.0) THEN
         BDIM(1)=NINT(FLOAT(BDIM(2))*RATIO)
      ENDIF
      BCONV(1) = 1
      BCONV(4) = 1
      BCONV(2) = DXMIN
      BCONV(5) = DYMIN
* If user gave increment, recompute BDIM
      IF (BCONV(3).EQ.0.0) THEN
         BCONV(3) = SIZX/BDIM(1)
*keep input map conventions, if any.
         IF(ACONV(3).LT.0.0) THEN
            BCONV(2) = DXMAX
            BCONV(3)=-BCONV(3)
         ENDIF
      ELSEIF (BCONV(3).LT.0.0) THEN
         BDIM(1) = -1.0*SIZX/BCONV(3)
         BCONV(2) = DXMAX
      ELSE
         BDIM(1) = SIZX/BCONV(3)
      ENDIF
      IF (BCONV(6).EQ.0.0) THEN
         BCONV(6) = SIZY/BDIM(2)
         IF(ACONV(6).LT.0.0) THEN
            BCONV(5) = DYMAX
            BCONV(6)=-BCONV(6)
         ENDIF
      ELSEIF (BCONV(6).LT.0.0) THEN
         BDIM(2) = -1.0*SIZY/BCONV(6)
         BCONV(5) = DYMAX
      ELSE
         BDIM(2) = SIZY/BCONV(6)
      ENDIF
* Check size:
      IF (BDIM(1)*BDIM(2).GT.2048*2048) THEN
         CALL GAGOUT('Computed Output Image size has too many pixels,')
         CALL GAGOUT('Please Specify Output Image Header. STOP.')
         STOP
      ELSE
         WRITE(*,100)
         WRITE(*,101) BDIM(1),(BCONV(J),J=1,3)
         WRITE(*,102) BDIM(2),(BCONV(J),J=4,6)
      ENDIF
*
      RETURN
100   FORMAT ('Output Image header values:',
     $/,'Axis Npix     Ref          Val            Inc')
101   FORMAT ('  1 ',I4,1X,3(G12.4E2,1X))
102   FORMAT ('  2 ',I4,1X,3(G12.4E2,1X))
1000  FORMAT ('              Found ',I8,' Valid points')
      END
*<FF>
      SUBROUTINE BESTBOX(DXD,DYD,NDP,DXMIN,DXMAX,DYMIN,DYMAX)
      REAL*8 DXMIN,DXMAX,DYMIN,DYMAX,DXD(*),DYD(*)
      INTEGER NDP,I
*
      DXMIN=DXD(1)
      DXMAX=DXD(1)
      DYMIN=DYD(1)
      DYMAX=DYD(1)
      DO I=2,NDP
         DXMIN=MIN(DXMIN,DXD(I))
         DYMIN=MIN(DYMIN,DYD(I))
         DXMAX=MAX(DXMAX,DXD(I))
         DYMAX=MAX(DYMAX,DYD(I))
      ENDDO
      RETURN
      END
*<FF>
      SUBROUTINE GRIDLIN(A,MX,MY,MZ,ACONV,ATYPE,APROJ,AEPOC,ABLANK,
     $B,NX,NY,BCONV,BTYPE,BPROJ,BEPOC,BLANK,CODE,XB,YB)
C----------------------------------------------------------------------
C GREG	Stand alone subroutine
C	Regrid an input map in a given projection to
C	an output map in another projection, bilinear interpolation
C Arguments
C	A	R*4(*)	Input map of dimensions MX MY MZ
C	MX,MY,MZ   I	Size of A
C	ACONV	R*8(6)	Pixel conversion formulae
C			CONV(1) = Xref, CONV(2)=Xval, CONV(3)=Xinc
C	ATYPE	I	Type of projection
C	APROJ	R*8(3)	Projection constants
C			PROJ(1) = A0,	PROJ(2)=D0,	PROJ(3)=Angle
C       AEPOC   R*4     Epoch if Needed
C	B	R*4(*)	Output map of dimensions NX,NY,MZ
C	NX,NY	I
C	BCONV	R*8(6)
C	BTYPE	I
C	BPROJ	R*8(3)
C       BEPOC   R*4     Epoch if Needed
C	BLANK	R*4	Blanking value
C	CODE	I	Galactic to Equatorial (-1)
C			or Equatorial to Galactic (1)
C	XB,YB	R*8	Work arrays of dimension MX,MY
C----------------------------------------------------------------------
      INTEGER NX,NY,MX,MY,MZ,ATYPE,BTYPE,CODE
      REAL*4 A(MX,MY,MZ), B(NX,NY,MZ), BLANK, ABLANK(2), BVAL,EVAL
      REAL*4 AEPOC,BEPOC
      REAL*8 XB(NX,NY), YB(NX,NY),
     $ACONV(6), APROJ(6), BCONV(6), BPROJ(6)
*
      INTEGER IB,JB,IA,JA,K,I
      REAL*8 AXREF,AXVAL,AXINC, AYREF,AYVAL,AYINC
      REAL*8 BXREF,BXVAL,BXINC, BYREF,BYVAL,BYINC
      REAL*8 XA,YA
      REAL*4 XR,YR
*
* Preliminary processing
      BVAL=ABLANK(1)
      EVAL=ABLANK(2)
*
* I force the conversion formula for input to be adapted to the
* NINT function that will be used below. Now Xref=1.0
      AXREF = 1.0                        ! was ACONV(1)
      AXVAL = (1.0-ACONV(1))*ACONV(3)+ACONV(2)
      AXINC = ACONV(3)
      AYREF = 1.0                        ! was ACONV(4)
      AYVAL = (1.0-ACONV(4))*ACONV(6)+ACONV(5)
      AYINC = ACONV(6)
* the output is unchanged.
      BXREF = BCONV(1)
      BXVAL = BCONV(2)
      BXINC = BCONV(3)
      BYREF = BCONV(4)
      BYVAL = BCONV(5)
      BYINC = BCONV(6)
*
* Setup Output projection
      CALL GR8_PROJEC (BTYPE,BPROJ)
*
* Convert Output projection to Absolute coordinates
      DO JB=1,NY
         DO IB=1,NX
            XB(IB,JB) = (IB-BXREF)*BXINC + BXVAL
            YB(IB,JB) = (JB-BYREF)*BYINC + BYVAL
         ENDDO
      ENDDO
      CALL REL_TO_ABS(XB,YB,XB,YB,NX*NY)
*
* Change coordinate system. use new epoch also. Warning: codes are
* reversed since we start form output map!
      IF (CODE.EQ.1) THEN
         CALL GAL_EQU(XB,YB,XB,YB,AEPOC,NX*NY)
      ELSEIF (CODE.EQ.-1) THEN
         CALL EQU_GAL(XB,YB,BEPOC,XB,YB,NX*NY)
      ELSEIF (CODE.EQ.-2) THEN
         CALL EQU_EQU(XB,YB,BEPOC,XB,YB,AEPOC,NX*NY)
      ENDIF
*
* Setup Input projection
      CALL GR8_PROJEC (ATYPE,APROJ)
*
* Convert Absolute coordinates to input projection
      CALL ABS_TO_REL(XB,YB,XB,YB,NX*NY)
*
* Interpolate
*
* Take Blanking into account, version #1
      IF (EVAL.LT.0.0) THEN
* Loop over Output data points
         DO K = 1,MZ
            DO JB = 1,NY
               DO IB = 1,NX
*
* Find pixel coordinate in input map
                  XA = XB(IB,JB)
                  YA = YB(IB,JB)
                  XR = (XA-AXVAL)/AXINC+AXREF
                  YR = (YA-AYVAL)/AYINC+AYREF
                  IA = INT(XR)
                  JA = INT(YR)
                  XR = XR-FLOAT(IA)
                  YR = YR-FLOAT(JA)
*
* Avoid edges
                  IF (IA.LT.1 .OR. IA.GE.MX .OR.
     $            JA.LT.1 .OR. JA.GE.MY) THEN
                     B(IB,JB,K) = BLANK
                  ELSE
* Interpolate
                     B(IB,JB,K) = (1-XR)*(1-YR)*A(IA,JA,K)+
     $               XR*(1-YR)*A(IA+1,JA,K)+
     $               XR*YR*A(IA+1,JA+1,K)+
     $               (1-XR)*YR*A(IA,JA+1,K)
                  ENDIF
               ENDDO
            ENDDO
         ENDDO
      ELSE
* Loop over Output data points
         DO K = 1,MZ
            DO JB = 1,NY
               DO IB = 1,NX
                  XA = XB(IB,JB)
                  YA = YB(IB,JB)
                  XR = (XA-AXVAL)/AXINC+AXREF
                  YR = (YA-AYVAL)/AYINC+AYREF
                  IA = INT(XR)
                  JA = INT(YR)
                  XR = XR-FLOAT(IA)
                  YR = YR-FLOAT(JA)
                  IF (IA.LT.1 .OR. IA.GE.MX .OR.
     $            JA.LT.1 .OR. JA.GE.MY) THEN
                     B(IB,JB,K) = BLANK
                  ELSEIF (ABS(A(IA,JA,K)-BVAL).LE.EVAL) THEN
                     B(IB,JB,K) = BLANK
                  ELSEIF (ABS(A(IA+1,JA,K)-BVAL).LE.EVAL) THEN
                     B(IB,JB,K) = BLANK
                  ELSEIF (ABS(A(IA,JA+1,K)-BVAL).LE.EVAL) THEN
                     B(IB,JB,K) = BLANK
                  ELSEIF (ABS(A(IA+1,JA+1,K)-BVAL).LE.EVAL) THEN
                     B(IB,JB,K) = BLANK
                  ELSE
                     B(IB,JB,K) = (1-XR)*(1-YR)*A(IA,JA,K)+
     $               XR*(1-YR)*A(IA+1,JA,K)+
     $               XR*YR*A(IA+1,JA+1,K)+
     $               (1-XR)*YR*A(IA,JA+1,K)
                  ENDIF
               ENDDO
            ENDDO
         ENDDO
      ENDIF
      END
