* Last processed by NICE on 12-Jun-2000 15:51:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      PROGRAM FILL_CUBE
C----------------------------------------------------------------------
C GILDAS	Stand alone program
C	Resample a partially filled input cube to a different
C	finer filled one...
C Subroutines
C	(IMAGE),(SIC)
C	REGRID
C----------------------------------------------------------------------
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:errcod.inc'
      INCLUDE 'inc:memory.inc'
      REAL*8 PI
      PARAMETER (PI=3.141592653589793D0)
*
      CHARACTER*80 NAMEY,NAMEX
      CHARACTER*4 METHOD
      INTEGER I,J,N,NDP,IER,LENC,SIC_GETVM
      LOGICAL ERROR,AUTO,EXTRAPOL,LOGRAN(7)
      REAL*8 CONV(3,2)
      INTEGER DIM(2),NXY,BLC(4),TRC(4)
      INTEGER IPX, IPY, IPI, IPW, ADDR, XD,YD,ZD,IA,JA, POINTER
*
      DATA LOGRAN/7*.FALSE./, BLC/4*0/, TRC/4*0/
*
      CALL GILDAS_OPEN
      CALL GILDAS_CHAR('Y_NAME$',NAMEY)
      CALL GILDAS_CHAR('X_NAME$',NAMEX)
      CALL GILDAS_INTE('PIXELS$',DIM,2)
      CALL GILDAS_LOGI('MODE$',AUTO,1)
      IF (.NOT.AUTO) THEN
         CALL GILDAS_DBLE('AXIS_1$',CONV(1,1),3)
         CALL GILDAS_DBLE('AXIS_2$',CONV(1,2),3)
      ENDIF
      CALL GILDAS_CHAR('METHOD$',METHOD)
      CALL GILDAS_LOGI('EXTRAPOLATE$',EXTRAPOL,1)
      CALL GILDAS_CLOSE
*
      N = LENC(NAMEY)
      CALL SIC_PARSEF(NAMEY(1:N),Y_FILE,' ','.gdf')
      CALL GDF_GEIS (Y_ISLO,ERROR)
      IF (.NOT.ERROR) CALL
     $GDF_REIS (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) 'F-FILL_CUBE,  Cannot read input file'
         GOTO 100
      ENDIF
      IF (Y_FORM.NE.FMT_R4) THEN
         WRITE(6,*) 'E-FILL_CUBE,  Only real images supported'
         GOTO 100
      ENDIF
      CALL GDF_READY(Y_ISLO,ERROR)
*
* Now define output
      CALL GDF_CHYX
      N = LENC(NAMEX)
      CALL SIC_PARSEF(NAMEX(1:N),X_FILE,' ','.gdf')
      IF (DIM(1).LT.1 .OR. DIM(2).LT.1) THEN
         WRITE(6,*) 'F-FILL_CUBE,  Invalid dimensions'
         GOTO 100
      ELSE
         X_DIM(1) = DIM(1)
         X_DIM(2) = DIM(2)
      ENDIF
      IF (AUTO) THEN
         WRITE(6,*) 'W-FILL_CUBE,  Computing conversion formulae'
         X_REF1 = 0.5
         X_VAL1 = (0.5-Y_REF1)*Y_INC1 + Y_VAL1 ! Pixel 0.5 of input
         X_INC1 = Y_DIM(1)*Y_INC1/X_DIM(1)     ! Total scale
         X_REF2 = 0.5
         X_VAL2 = (0.5-Y_REF2)*Y_INC2 + Y_VAL2 ! Pixel 0.5 of input
         X_INC2 = Y_DIM(2)*Y_INC2/X_DIM(2)     ! Total scale
      ELSE
         DO J=1,2
            DO I=1,3
               X_CONV(I,J) = CONV(I,J)
            ENDDO
         ENDDO
      ENDIF
*
* Method of solution
      IF (METHOD.EQ.'SLOW') THEN
         WRITE(6,*) 'W-FILL_CUBE,  Using SLOW general method'
      ELSE
         METHOD = 'FAST'
         WRITE(6,*)
     $   'W-FILL_CUBE,  Using FAST method, input blanking ignored'
      ENDIF
*
* Create output image
      ERROR = .FALSE.
      X_SIZE = X_DIM(1)*X_DIM(2)*X_DIM(3)
      CALL GDF_GEIS (X_ISLO,ERROR)
      IF (.NOT.ERROR) CALL GDF_WRITX(X_ISLO,ERROR)
      IF (.NOT.ERROR) CALL
     $GDF_CRIS (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) 'F-FILL_CUBE,  Cannot create X file'
         WRITE(6,*) X_FILE
         GOTO 100
      ENDIF
*
      IF (METHOD.EQ.'FAST') THEN
         DO J=1,X_DIM(4)
            BLC(4) = J
            TRC(4) = J
            DO I=1,X_DIM(3)
               BLC(3) = I
               TRC(3) = I
               CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
               IF (ERROR) GOTO 100
               IPX = POINTER(X_ADDR,MEMORY)
               CALL GDF_GEMS (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
               IF (ERROR) GOTO 100
               IPY = POINTER(Y_ADDR,MEMORY)
               CALL FILL_CUBE_FAST (
     $         MEMORY(IPY),Y_DIM(1),Y_DIM(2),Y_CONV,
     $         MEMORY(IPX),X_DIM(1),X_DIM(2),X_CONV,
     $         X_BVAL)
               CALL GDF_FRMS(X_MSLO,ERROR)
               CALL GDF_FRMS(Y_MSLO,ERROR)
            ENDDO
         ENDDO
      ELSE
         NXY = Y_DIM(1)*Y_DIM(2)
         IER = SIC_GETVM(5*NXY,ADDR)
         IF (IER.NE.1) GOTO 99
         XD = POINTER(ADDR,MEMORY)
         YD = XD + NXY
         ZD = YD + NXY
         IA = ZD + NXY
         JA = IA + NXY
         X_EVAL = MAX(X_EVAL,ABS(X_BVAL*1E-7))
*
         BLC(4) = MAX(1,Y_DIM(4)/2)
         TRC(4) = BLC(4)
         BLC(3) = MAX(1,Y_DIM(3)/2)
         TRC(3) = BLC(3)
         CALL GDF_GEMS (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
         IF (ERROR) GOTO 100
         IPY = POINTER(Y_ADDR,MEMORY)
*
         CALL INIT_CUBE_SLOW (MEMORY(IPY),
     $   Y_DIM(1),Y_DIM(2),Y_CONV,
     $   X_DIM(1),X_DIM(2),X_CONV,
     $   X_BVAL,X_EVAL,
     $   MEMORY(XD),MEMORY(YD),MEMORY(ZD),MEMORY(IA),MEMORY(JA),
     $   NDP,ERROR)
         IF (ERROR) GOTO 100
         CALL GDF_FRMS(Y_MSLO,ERROR)
         LOGRAN(3) = .FALSE.
         LOGRAN(7) = EXTRAPOL
*
         IER = SIC_GETVM (13*NDP,ADDR)
         IF (IER.NE.1) GOTO 99
         IPW = POINTER(ADDR,MEMORY)
         IER = SIC_GETVM (35*NDP,ADDR)
         IF (IER.NE.1) GOTO 99
         IPI = POINTER(ADDR,MEMORY)
*
         DO J=1,Y_DIM(4)
            BLC(4) = J
            TRC(4) = J
            DO I=1,Y_DIM(3)
               BLC(3) = I
               TRC(3) = I
               CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
               IF (ERROR) GOTO 100
               IPX = POINTER(X_ADDR,MEMORY)
               CALL GDF_GEMS (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
               IF (ERROR) GOTO 100
               IPY = POINTER(Y_ADDR,MEMORY)
               CALL FILL_CUBE_SLOW (
     $         MEMORY(IPY),Y_DIM(1),Y_DIM(2),
     $         MEMORY(IPX),X_DIM(1),X_DIM(2),
     $         X_BVAL,
     $         MEMORY(XD),MEMORY(YD),MEMORY(ZD),MEMORY(IA),MEMORY(JA),
     $         MEMORY(IPI),MEMORY(IPW),NDP,
     $         LOGRAN,ERROR)
               IF (ERROR) GOTO 100
               CALL GDF_FRMS(X_MSLO,ERROR)
               CALL GDF_FRMS(Y_MSLO,ERROR)
            ENDDO
         ENDDO
      ENDIF
      CALL GDF_FRIS(X_ISLO,ERROR)
      CALL GDF_FRIS(Y_ISLO,ERROR)
      STOP 'E-FILL_CUBE,  Successful completion'
*
99    WRITE(6,*) 'F-FILL_CUBE,  Memory allocation failure'
100   CALL SYSEXI(FATALE)
      END
*<FF>
      SUBROUTINE INIT_CUBE_SLOW (A,MX,MY,ACONV,
     $NX,NY,BCONV,BLANK,EBLANK,
     $XD,YD,ZD,IA,JA,NDP,ERROR)
C----------------------------------------------------------------------
C GILDAS	Stand alone subroutine
C	Regrid an input cube with blanking values to an output
C	regular cube on a finer grid.
C Arguments
C	A		R*4(*)	Input map of dimensions MX MY
C	MX,MY		I
C	ACONV		R*8(6)	Pixel conversion formulae
C			CONV(1) = Xref, CONV(2)=Xval, CONV(3)=Xinc
C	NX,NY		I	Output dimensions
C	BCONV		R*8(6)
C	BLANK		R*4	Blanking value
C	EBLANK		R*4	Blanking tolerance
C	XD,YD,ZD	R*4	Work arrays of dimension MX*MY
C	IA,JA		I*4	Work arrays of dimensions MX*MY
C	NDP		I*4	Number of selected points
C----------------------------------------------------------------------
      INTEGER NX,NY,MX,MY,NDP
      REAL*4 A(MX,MY), BLANK, EBLANK
      REAL*4 XD(NX*NY),YD(NX*NY),ZD(NX*NY)
      REAL*8 ACONV(6), BCONV(6)
      INTEGER IA(MX*MY),JA(MX*MY)
      LOGICAL ERROR
*
      INTEGER I,J,NCP
      REAL*8 AXREF,AXVAL,AXINC, AYREF,AYVAL,AYINC
      REAL*8 BXREF,BXVAL,BXINC, BYREF,BYVAL,BYINC
*
* Preliminary processing
      NCP = 4
      AXREF = ACONV(1)
      AXVAL = ACONV(2)
      AXINC = ACONV(3)
      AYREF = ACONV(4)
      AYVAL = ACONV(5)
      AYINC = ACONV(6)
      BXREF = BCONV(1)
      BXVAL = BCONV(2)
      BXINC = BCONV(3)
      BYREF = BCONV(4)
      BYVAL = BCONV(5)
      BYINC = BCONV(6)
*
* Load XD,YD,ZD arrays with current plane (or mask array)
      NDP = 0
      DO J=1,MY
         DO I=1,MX
            IF (ABS(A(I,J)-BLANK).GT.EBLANK) THEN
               NDP = NDP +1
               XD(NDP) = (I-AXREF)*AXINC + AXVAL
               YD(NDP) = (J-AYREF)*AYINC + AYVAL
               ZD(NDP) = A(I,J)
               IA(NDP) = I
               JA(NDP) = J
            ENDIF
         ENDDO
      ENDDO
      IF (NDP.LE.0) THEN
         WRITE(6,*) 'F-SLOW,  No data point in input'
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Initialize Gridding
      CALL GRIDINI(NX,BXREF,BXVAL,BXINC, NY,BYREF,BYVAL,BYINC)
      END
*<FF>
      SUBROUTINE FILL_CUBE_SLOW (A,MX,MY,
     $B,NX,NY,BLANK,
     $XD,YD,ZD,IA,JA,
     $IPI,IPW,NDP,
     $LOGRAN,ERROR)
C----------------------------------------------------------------------
C GILDAS	Stand alone subroutine
C	Regrid an input cube with blanking values to an output
C	regular cube on a finer grid.
C Arguments
C	A		R*4(*)	Input map of dimensions MX MY
C	MX,MY		I
C	B		R*4(*)	Output map of dimensions NX,NY
C	NX,NY		I
C	BLANK		R*4	Blanking value
C	XD,YD,ZD	R*4	Work arrays of dimension MX*MY
C	IA,JA		I*4	Work arrays of dimensions MX*MY
C	IPI		I*4	Work arrays
C	IPW		I*4	Work arrays
C	NDP		I*4	Number of good points
C	LOGRAN		L(7)	Control interpolation behaviour
C----------------------------------------------------------------------
      INTEGER NX,NY,MX,MY,NDP
      REAL*4 A(MX,MY), B(NX,NY), BLANK
      REAL*4 XD(NX*NY),YD(NX*NY),ZD(NX*NY)
      INTEGER IPI(*),IPW(*)
      INTEGER IA(MX*MY),JA(MX*MY)
      LOGICAL LOGRAN(7),ERROR
*
      EXTERNAL DUMMY
      INTEGER I,NCP
*
      NCP = 4
      CALL GRIDSET (NCP,BLANK,LOGRAN)
*
* Process plane
      DO I=1,NDP
         ZD(I) = A(IA(I),JA(I))
      ENDDO
      CALL GRIDRAN (XD,YD,ZD,NDP,IPW,IPI,B
     $,DUMMY,DUMMY,DUMMY,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) 'F-FILL_CUBE_SLOW,  Failure in GRIDRAN ...'
         RETURN
      ENDIF
      LOGRAN(3) = .TRUE.                 ! Repeat for next one
      END
*
      SUBROUTINE DUMMY
      END
*<FF>
      SUBROUTINE FILL_CUBE_FAST (A,MX,MY,ACONV,
     $B,NX,NY,BCONV,BLANK)
C----------------------------------------------------------------------
C GILDAS 	Stand alone subroutine
C	Regrid an input data cube on a finer one for the 2 first input
C	dimensions.
C Arguments
C	A		R*4(*)	Input map of dimensions MX MY
C	MX,MY		I
C	ACONV		R*8(6)	Pixel conversion formulae
C			CONV(1) = Xref, CONV(2)=Xval, CONV(3)=Xinc
C	B		R*4(*)	Output map of dimensions NX,NY
C	NX,NY		I
C	BCONV		R*8(6)
C	BLANK	R*4	Blanking value
C----------------------------------------------------------------------
      INTEGER NX,NY,MX,MY
      REAL*4 A(MX,MY), B(NX,NY), BLANK
*
      INTEGER IB,JB,IA,JA
      REAL*8 AXREF,AXVAL,AXINC, AYREF,AYVAL,AYINC, ACONV(6)
      REAL*8 BXREF,BXVAL,BXINC, BYREF,BYVAL,BYINC, BCONV(6)
      REAL*8 XA,YA
      REAL*4 APLUS,AMOIN,AZERO,XR,YR
*
* Preliminary processing
      AXREF = ACONV(1)
      AXVAL = ACONV(2)
      AXINC = ACONV(3)
      AYREF = ACONV(4)
      AYVAL = ACONV(5)
      AYINC = ACONV(6)
      BXREF = BCONV(1)
      BXVAL = BCONV(2)
      BXINC = BCONV(3)
      BYREF = BCONV(4)
      BYVAL = BCONV(5)
      BYINC = BCONV(6)
*
* Loop over Output data points
      DO JB = 1,NY
         DO IB = 1,NX
*
* Find pixel coordinate in input map
            XA = (IB-BXREF)*BXINC + BXVAL
            YA = (JB-BYREF)*BYINC + BYVAL
            XR = (XA-AXVAL)/AXINC + AXREF
            YR = (YA-AYVAL)/AYINC + AYREF
            IA = NINT (XR)
            JA = NINT (YR)
            XR = XR - IA
            YR = YR - JA
*
* Avoid edges
            IF (IA.LE.1 .OR. IA.GE.MX .OR. JA.LE.1 .OR. JA.GE.MY) THEN
               B(IB,JB) = BLANK
            ELSE
*
* Interpolate (X or Y first, does not matter in this case)
               APLUS = ( (A(IA+1,JA+1)+A(IA-1,JA+1)-2.*A(IA,JA+1))
     $         *XR + A(IA+1,JA+1)-A(IA-1,JA+1) )*XR*0.5
     $         + A(IA,JA+1)
               AZERO = ( (A(IA+1,JA)+A(IA-1,JA)-2.*A(IA,JA))
     $         *XR + A(IA+1,JA)-A(IA-1,JA) )*XR*0.5
     $         + A(IA,JA)
               AMOIN = ( (A(IA+1,JA-1)+A(IA-1,JA-1)-2.*A(IA,JA-1))
     $         *XR + A(IA+1,JA-1)-A(IA-1,JA-1) )*XR*0.5
     $         + A(IA,JA-1)
* Then Y (or X)
               B(IB,JB) = ( (APLUS+AMOIN-2.*AZERO)
     $         *YR + APLUS-AMOIN )*YR*0.5 + AZERO
            ENDIF
         ENDDO
      ENDDO
      END
