* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      SUBROUTINE RESMAP8(NODIN,LHEAD1,NGRID1,NCOLS1,NROWS1,DUMMYMIN,
     $DUMMYMAX,UNITFACT1,XINC1,YINC1,BMSEP,
     $THETADEG,UNBAL,STRING,ISTART,ISTAT,MASK,RESTORE_MASK,
     $RESTORE_FILTER,Z1,WORK8,SCAN8,CSCAN8,FFUN8,RFUN8,NSINCS)
      INCLUDE 'const.inc'
*
C Parameters for subroutine RESMAP8
*
      INTEGER*4 NODIN, LHEAD1, NGRID1, NCOLS1,NROWS1
      REAL*8    DUMMYMIN, DUMMYMAX
      REAL*8    UNITFACT1, XINC1, YINC1, BMSEP,THETADEG, UNBAL
      INTEGER   NSINCS
      CHARACTER*80 STRING
      INTEGER   ISTART, ISTAT
      REAL*8    SCAN8(NCOLS1), CSCAN8(NCOLS1),
     $FFUN8(2*NCOLS1), RFUN8(2*NCOLS1), Z1(NCOLS1*NROWS1),
     $WORK8(NCOLS1*NROWS1)
      INTEGER   MASK(1)
      LOGICAL   RESTORE_MASK, RESTORE_FILTER
*
* Local data
      REAL*8    ROTMAT(9), PTSEP,
     $XINC1DEG, YINC1DEG, FNORM, RNORM, SMALLIM
      INTEGER*4 LPROMT, LBLANK, LEN1, LEN2,
     $ICENTR, I, IROW, NGRIDLI, NCOLSLI, NROWSLI,
     $MODLRMAP, MODLRSCN1, MODLRSCN2
      CHARACTER*80  PROMT, BLANK
*
c This is a program to take a map made with a dual-beam system, and to
c re-weight the spatial fourier frequencies via a convolution along the
c line joining the beams, to simulate a map made with a single beam ("restore"
c the map). The output map will be normalised to the stronger of the two
c beams.
c SCAN8 is a vector to hold one row of the map Zl to be convolved (restored),
c CSCAN8 ls a vector holding the row after its convolution,
c RFUN8 is the restoring function and
c FFUN8 is a filter function used to filter out from the original data
c spatial frequencies which cannot be astronomical - i.e. wavelengths
c equal to and harmonics of the original beam spacing BMSEP.
*
c Written by Darrel Emerson and Chris Salter ?
c real*8 version 08.07.1993, Robert Zylka, MPIfR Bonn
*
      LPROMT = 7
      PROMT(1:LPROMT+2) = 'RESMAP8> '
      LBLANK = LPROMT + 2
      DO I = 1, LBLANK
         BLANK(I:I) = ' '
      ENDDO
*
      SMALLIM = 1.0D-35
*
      IF (DABS(XINC1) .LE. SMALLIM) THEN
         PRINT *, PROMT(1:LPROMT), '> XINC = O !'
         ISTAT = 1
         RETURN
      ENDIF
      IF (DABS(YINC1) .LE. SMALLIM) THEN
         PRINT *, PROMT(1:LPROMT), '> YINC = O !'
         ISTAT = 2
         RETURN
      ENDIF
      IF (DABS(BMSEP) .LE. SMALLIM) THEN
         PRINT *, PROMT(1:LPROMT), '> Beam separation = 0 !'
         ISTAT = 3
         RETURN
      ENDIF
*
c      CALL CHKDUMMY8(PROMT,LPROMT,BLANK,LBLANK,Z1,LHEAD1,NGRID1,
c     $DUMMYMIN,DUMMYMAX,ISTAT)
c      IF (ISTAT .NE. 0) RETURN
*
c      IF (NODIN .GE. 1) CALL DUMPMT(nco*nli,Z1,AMAT)
*
      NCOLSLI = NCOLS1                   ! number of columns/rows in working matrix WO
      NROWSLI = NROWS1                   ! Can be changed in LRMAP8 if THETADEG .ne. O
      NGRIDLI = NGRID1
      XINC1DEG = XINC1*UNITFACT1 / 3600.0D0
      YINC1DEG = YINC1*UNITFACT1 / 3600.0D0
*
      MODLRMAP = 1
      CALL LRMAP8(BLANK,LBLANK,ROTMAT,NODIN,LHEAD1,NCOLSLI,NROWSLI,
     $NGRIDLI,DUMMYMIN,DUMMYMAX,XINC1DEG,YINC1DEG,THETADEG,
     $MODLRMAP,ISTAT,Z1,WORK8)
      IF (ISTAT .NE. 0) RETURN
*
c MODLRMAP=1 for LRMAP8 uses THETADEG; MODLRMAP=2 uses the matrices set
c already in Z1 and WORK8 (WORK8 is set during MODLRMAP=1). In both cases if
c THETADEG = 0 Z1 is copied directly to WORK8 or vice versa.
c If THETADEG .ne. 0,
c subroutine LRMAP8 uses MRETAB to rotate the map angle THETADEG clockwise.
c Note that the XINC1 (z1(45)*3600/UNITFACTl) and YINC1 (zl(46)*3600/UNITFACTl)
c after rotation are not changed in the map header, but NCOLSlI
c and NROWSlI may be different due to XTENT and YTENT being calculated by
c subroutine BULGE, called by LRMAP8 in MODLRMAP=1.
*
*
c ICENTR is the centre element of the convolution function, which has size
c LEN2.
      LEN1 = NCOLSLI
      LEN2 = 2*NCOLSLI - 1
      ICENTR = NCOLSLI
      PTSEP = DABS(BMSEP) / XINC1
*
c If the two beams are equal in amplitude the input maps will be pre-filtered
c to remove spatial frequencies which contain only noise, i.e. spatial
c frequencies corresponding to the original beam spacing and harmonics of
c this.
*
      IF (DABS(DABS(UNBAL) - 1.0D0) .LE. SMALLIM) THEN
         CALL CALFFUN8(LEN1,LEN2,ICENTR,PTSEP,PI,ISTAT,FFUN8)
         IF (ISTAT .NE. 0) RETURN
      ENDIF
*
      CALL CALRFUN8(LEN1,LEN2,ICENTR,PTSEP,UNBAL,PI,ISTAT,RFUN8,NSINCS)
      IF (ISTAT .NE. 0) RETURN
      MODLRSCN1 = 1
      MODLRSCN2 = 2
      RNORM = 2.0D0
*
c FNORM is the normalisation factor for the filtering and RNORM for the
c restoring function, which were calculated in the subroutines CALFFUN8 and
c CALRFUN8 respectively.
*
c The convolution is along a row, asuming the two beams were scanned
c along a longitude.
c LRSCAN8 dumps (MODLRSCN1=1) the irow'th row of the map WORK8 into the vector
c SCAN8.
c CONV8 convolves a vector with the particular convolution function. The array
c sizes are LEN1 and LEN2. ICENTR is the centre element of the particular
c convolution function used.
c LRSCAN8 dumps (MODLRSCN2=2) a vector into the irow'th row of the map WORK8.
*
      IF (DABS(DABS(UNBAL)-1D0).LE.SMALLIM.AND.RESTORE_FILTER) THEN
c         PRINT *, BLANK(1:LBLANK), 'Filtering and restoring the map ...'
         FNORM = 1.0D0
         DO IROW = 1, NROWSLI
            CALL LRSCAN8(SCAN8,LHEAD1,NCOLSLI,NROWSLI,IROW,MODLRSCN1,
     $      WORK8)
            CALL CONV8(SCAN8,CSCAN8,LEN1,FFUN8,LEN2,ICENTR,FNORM,
     $      DUMMYMIN,DUMMYMAX,MASK((IROW-1)*NCOLSLI+1),.FALSE.)
            CALL CONV8(CSCAN8,SCAN8,LEN1,RFUN8,LEN2,ICENTR,RNORM,
     $      DUMMYMIN,DUMMYMAX,MASK((IROW-1)*NCOLSLI+1),RESTORE_MASK)
            CALL LRSCAN8(SCAN8,LHEAD1,NCOLSLI,NROWSLI,IROW,MODLRSCN2,
     $      WORK8)
         ENDDO
      ELSE
c         PRINT *, BLANK(1:LBLANK), 'Restoring the map ...'
         DO IROW = 1, NROWSLI
            CALL LRSCAN8(SCAN8,LHEAD1,NCOLSLI,NROWSLI,IROW,MODLRSCN1,
     $      WORK8)
            CALL CONV8(SCAN8,CSCAN8,LEN1,RFUN8,LEN2,ICENTR,RNORM,
     $      DUMMYMIN,DUMMYMAX,MASK((IROW-1)*NCOLSLI+1),RESTORE_MASK)
            CALL LRSCAN8(CSCAN8,LHEAD1,NCOLSLI,NROWSLI,IROW,MODLRSCN2,
     $      WORK8)
         ENDDO
      ENDIF
*
      MODLRMAP = 2
      CALL LRMAP8(BLANK,LBLANK,ROTMAT,NODIN,LHEAD1,NCOLS1,NROWS1,NGRID1,
     $DUMMYMIN,DUMMYMAX,XINC1DEG,YINC1DEG,THETADEG,MODLRMAP,
     $ISTAT,Z1,WORK8)
*
c MODLRMAP=2 in LRMAP8 means load Z1 back from WORK8, using the matrices set
c already in the maps. THETADEG is only used to check for a finite rotation
c in this case.
*
ccc      if (bmsep .lt. 0.0d0) call FLIPINX(lhead1,ncols1,nrows1)
*
c This is a facility to flip the map left-to-right, to allow for maps which
c may have been made in az-el co-ordinates, and need to be reversed.
c Now retrieve the 9-element matrix stored for safety in ANAT at the beginning.
*
      IF (NODIN .GE. 1) THEN
c         CALL UNDUMPMT(IDIMZ1,Z1,AMAT)
c Now write the restoring parameters in the map title
c         Z1(21) = 'res '
c         READ(STRING(ISTART:ISTART+7),'(2a4)',ERR=10) Z1(22), Z1(23)
      ENDIF
*
      RETURN
      END
*
      SUBROUTINE CHKDUMMY8(PROMT,LPROMT,BLANK,LBLANK,SCANMAP,LHEAD,
     $NOOFDAT,DUMMYMIN,DUMMYMAX,ISTAT)
      IMPLICIT NONE
      INTEGER*4 LPROMT, LBLANK, LHEAD, NOOFDAT, ISTAT, I
      REAL*8    SCANMAP(LHEAD+NOOFDAT), DUMMYMIN, DUMMYMAX
      CHARACTER*80 PROMT, BLANK
*
c Subroutine CHKDUMMY cheks, if blanked/masked points dummys are still in
c the data.
c Written by Robert Zylka, MPIfR Bonn
c real*8 version 06.07.1993, R. Zylka
*
      PRINT *, PROMT(1:LPROMT), '> Searching for blanked data ...'
      DO I = LHEAD + 1, LHEAD + NOOFDAT
         IF (SCANMAP(I).GE.DUMMYMIN .AND. SCANMAP(I).LE.DUMMYMAX) THEN
            PRINT *, BLANK(1:LBLANK),
     $      'There are still blanked data (DUMMYS) in the map!'
            ISTAT = 4
            RETURN
         ENDIF
      ENDDO
*
      ISTAT = 0
*
      RETURN
      END
*
      SUBROUTINE DUMPMT(IDIMAMAP,AMAP,AMAT)
*
      INTEGER*4  IDIMAMAP
      REAL*8     AMAP(IDIMAMAP), AMAT(9)
      INTEGER*4  I
*
c DUMPMT dumps the co-ordinate matrix from the map for safety. Routine
c UNDUMPMT retrieves it. Note that DUMPMT can only store a single 9-element
c matrix.
c real*8 version 13.06.1993, Robert Zylka, MPIfR Bonn
*
      DO I = 1, 9
         AMAT(I) = AMAP(I+31)
      ENDDO
*
      RETURN
      END
*
      SUBROUTINE UNDUMPMT(IDIMAMAP,AMAP,AMAT)
      INTEGER*4 IDIMAMAP, I
      REAL*8 AMAP(IDIMAMAP), AMAT(9)
*
c DUMPMT dumps the co-ordinate matrix from the map for safety. Routine
c UNDUMPMT retrieves it. Note that DUMPMT can only store a single 9-element
c matrix.
c real*8 version 13.06.1993, Robert Zylka, MPIfR Bonn
*
      DO I = 1, 9
         AMAP(I+31) = AMAT(I)
      ENDDO
      RETURN
      END
*
*
      SUBROUTINE LRMAP8(BLANK,LBLANK,ROTMAT,NODIN,LHEAD,NCOLS,NROWS,
     $NGRID,DUMMYMIN,DUMMYMAX,XINCDEG,YINCDEG,
     $THETADEG,MODE,ISTAT,Z1,WORK8)
*
*
C Parameters for restore
*
      REAL*8    Z1(1),WORK8(1)
*
*
      REAL*8    ROTMAT(9), DUMMYMIN,
     $DUMMYMAX, XINCDEG, YINCDEG, THETADEG, SMALLIM
*
      INTEGER*4 LBLANK, NODIN, LHEAD, NCOLS, NROWS, NGRID, MODE, ISTAT,I
      CHARACTER*80 BLANK
*
c Subroutine LRMAP8 sets up the working array WORK8.
c If the angle THETADEG [degrees] .ne. 0, LRMAP8 rotates the map AMAP through
c THETADEG clockwise into WORK8. With MODE=l, THETADEG is used and the matrix
c is set in WORK8 using MSETUP, then the coordinate-matrix is transferred to
c ROTMAT for direct use by MRETAB. With MODE=2, the matrix already set in
c ROTMAT is transposed to give TROTMAT directly.
c real*8 version 12.10.1993, R. Zylka
*
      SMALLIM = 1.0D-35
*
      IF (MODE .EQ. 1) THEN
cc         if (DABS(thetadeg) .lt. smallim) then
         DO I = LHEAD + 1, LHEAD + NGRID
            WORK8(I) = Z1(I)
         ENDDO
cc         else
cc            call ROTMAP8(blank,lblank,idimz1,z1,idimzw,work8,rotmat,
cc     1                   nodin,lhead,ncols,nrows,ngrid,xincdeg,yincdeg,
cc     2                   thetadeg,mode,istat)
cc            if (istat .ne. 0) RETURN
cc
cc            do i = lhead + 1, lhead + ngrid
cc               if (work8(i).ge.dummymin .and. work8(i).le.dummymax) then
cc                  work8(i) = 0.0d0
cc               endif
cc            enddo
cc         endif
      ELSE
cc         if (DABS(thetadeg) .lt. smallim) then
         DO I = LHEAD + 1, LHEAD + NGRID
            Z1(I) = WORK8(I)
         ENDDO
cc         else
cc            call ROTMAP8(blank,lblank,idimzw,work8,idimz1,z1,rotmat,
cc     1                   nodin,lhead,ncols,nrows,ngrid,xincdeg,yincdeg,
cc     2                   thetadeg,mode,istat)
cc            if (istat .ne. 0) RETURN
cc         endif
      ENDIF
*
      ISTAT = 0
*
      RETURN
      END
*
      SUBROUTINE BULGE(XTENT,YTENT,THETADEG)
      INCLUDE 'const.inc'
      REAL*8 XTENT, YTENT, THETADEG, ABSXTENT, ABSYTENT,
     $ABSTHETA, ANGL, ANG2, ANG3, ANG4, DIAG, SMALLIM
*
c real*8 version 02.07.1993, Robert Zylka, MPIfR Bonn
*
*
      SMALLIM = 1.0D-35
      ABSXTENT = DABS(XTENT)
      ABSYTENT = DABS(YTENT)
      ABSTHETA = DABS(THETADEG)
      DIAG = DSQRT(XTENT*XTENT + YTENT*YTENT)
      IF (XTENT .GE. SMALLIM) THEN
         ANGL = DATAN(YTENT/XTENT)/DEG_TO_RAD
      ELSE
         ANGL = 90.0D0
      ENDIF
      ANG2 = -ANGL
      ANG3 = ANGL - ABSTHETA
      ANG4 = ANG2 - ABSTHETA
      YTENT = DABS(DIAG*DSIN(ANG4*DEG_TO_RAD))
      XTENT = DABS(DIAG*DCOS(ANG3*DEG_TO_RAD))
      IF (XTENT .LT. ABSXTENT) XTENT = ABSXTENT
      IF (YTENT .LT. ABSYTENT) YTENT = ABSYTENT
c This checks that XTENT and YTENT have not actually been made smaller, as
c may have happened for very small values of THETADEG and resulting from
c rounding errors in the computer.
*
      RETURN
      END
*
      SUBROUTINE TRANSP(AMAT,TAMAT)
      REAL*8 AMAT(9), TAMAT(9)
*
c Subroutine TRANSP puts the transpose of the 3x3 matrix AMAT into TAMAT.
c real*8 version 24.08.1992, Robert Zylka, MPIfR Bonn
*
      TAMAT(1) = AMAT(1)
      TAMAT(5) = AMAT(5)
      TAMAT(9) = AMAT(9)
      TAMAT(2) = AMAT(4)
      TAMAT(3) = AMAT(7)
      TAMAT(4) = AMAT(2)
      TAMAT(6) = AMAT(8)
      TAMAT(7) = AMAT(3)
      TAMAT(8) = AMAT(6)
*
      RETURN
      END
*
      SUBROUTINE LRSCAN8(SCAN8,LHEAD,NCOLS,NROWS,IROW,MODE,WORK8)
*
      INTEGER   NCOLS
      REAL*8    SCAN8(NCOLS)
      INTEGER*4 LHEAD, NROWS, IROW, MODE, ISTART, IEND,
     $K, I
      REAL*8    WORK8(1)
*
c Subroutine LRSCAN8 dumps for MODE=l row IROW from the map WORK8 into the
c vector SCAN8, for MODE=2 loads vector into row IROW of map WORK8.
c WORK8 is set up in LRMAP8 during MODE=l.
c real*8 version 07.07.1993, Robert Zylka, MPIfR Bonn
*
      ISTART = LHEAD + (IROW-1)*NCOLS + 1
      IEND = ISTART + NCOLS - 1
      K = 0
      IF (MODE .EQ. 1) THEN
         DO I = ISTART, IEND
            K = K + 1
            SCAN8(K) = WORK8(I)
         ENDDO
      ELSE
         DO I = ISTART, IEND
            K = K + 1
            WORK8(I) = SCAN8(K)
         ENDDO
      ENDIF
      RETURN
      END
*
      SUBROUTINE CALFFUN8(LEN1,LEN2,ICENTR,PTSEP,PI,ISTAT,FFUN8)
*
      REAL*8   PTSEP, PI, SINCS, SINCSNEG, SUMT, DX,
     $SN, SUM, PHI_IJ, SINC8
      INTEGER*4 LEN1, LEN2, ICENTR, NSINCS, ISTAT,
     $NSINCS2, I, J
      REAL*8    FFUN8(1)
*
*
c SUMT is to be a normalizing factor for the convolution function.
c Subroutine CALFFUN8 is a modified version of CONF2. CALFUN8 calculates a
c convolution function which will filter out spatial frequencies of
c wavelength equal to the spacing of the original two-beam system, and
c harmonics of this frequency. If any of these spatial frequencies exist in
c the map, they must result from noise or interference, since the two-beam
c observation should have acted as a spatial frequency filter. Hopefully if
c a map is restored after the application of filtering by this convolution
c function, the resultant signal-to-noise ratio and any effects of incorrect
c restoration will be substantially improved.
c Note that CALFFUN8 calculates an odd number (=2*nsincs+1) of delta
c functions, while CALRFUN8 works on an even number. CALFFUN8 produces a
c symetric function, CALRFUN8 an anti-symetric function.
c SUMT is to be a normalizing factor for the convolution function.
c Written by Darrel Emerson ?
c real*8 version 08.07.1993, Robert Zylka, MPIfR Bonn
*
      NSINCS = INT(DFLOAT(LEN1) / PTSEP + 0.5D0)
      NSINCS2 = 2*NSINCS + 1
      SINCS = DFLOAT(NSINCS + 1)
      SINCSNEG = -SINCS
      SUMT = 0.0D0                       ! a normalizing factor for the filtering func
      DO I = 1, LEN2
         DX = DFLOAT(I - LEN1)
         SN = SINCSNEG*PTSEP
         SUM = 0.0D0
         DO J = 1, NSINCS2
            SN = SN + PTSEP
            PHI_IJ = (DX-SN)*PI
c                      PHI_IJ is the phase angle at point i relative to SINC8 j
            SUM = SUM - SINC8(PHI_IJ)
         ENDDO
         FFUN8(I) = SUM
         SUMT = SUMT + SUM
      ENDDO
*
      SUMT = -SUMT*DFLOAT(LEN1) / DFLOAT(LEN2)
c This is to make the normalizing factor sumt for the array of sincs such
c that the sum over the map size length is equal to -1.
      DO I = 1, LEN2
         FFUN8(I) = FFUN8(I) / SUMT
      ENDDO
*
      FFUN8(ICENTR) = FFUN8(ICENTR) + 1.0D0
c The centre element of FFUN8 is supposed to be 1.0 greater than it would be
c from the sum of the various sincs.
*
      RETURN
      END
*
      SUBROUTINE CALRFUN8(LEN1,LEN2,ICENTR,PTSEP,UNBAL,PI,ISTAT,RFUN8,
     $NSINCS)
*
      REAL*8    PTSEP, UNBAL, PI, SINCS, ABSUNBAL, DX, ALPHA,
     $SN, SUM, PHI_IJ, SINC8, ANORM, X, SMALLIM
      INTEGER*4 LEN1, LEN2, ICENTR, NSINCS, ISTAT, NSINCS2, I, J
      REAL*8    RFUN8(1)
*
*
c Subroutine CALRFUN8 works out a convolution function for the 2-beam
c restoration program, where the 2-beam separation need not be an integral
c number of map sampling intervals. This subroutine interpolates to a more
c general grid using a SINC8 function sampled at information-theory limit.
c UNBAL is the ratio in amplitude of the two beams (e.g. the ratio left-hand
c beam amplitude over right-hand beam amplitude). Note that if the two beams
c are not equal in size, then the pre-filtering usually carried out in
c routine FILTSCAN or FILTMAP is suppressed.
c Written by Darrel Emerson ?
c real*8 version 08.07.1993, Robert Zylka, MPIfR Bonn
*
*
      SMALLIM = 1.0D-35
      IF (NSINCS.EQ.0) THEN
         NSINCS = INT(DFLOAT(LEN1) / PTSEP + 0.5D0)
      ENDIF
      NSINCS=MIN(INT(DFLOAT(LEN1) / PTSEP + 0.5D0),NSINCS)
      NSINCS2 = 2*NSINCS
      SINCS = DFLOAT(NSINCS + 1)
      ABSUNBAL = DABS(UNBAL)
      IF (ABSUNBAL .GT. 1.0D0) ABSUNBAL = 1.0D0 / ABSUNBAL
c******* Chris Salter modified this on Jun 13th 1985.****
c The necessary function is calculated with a factor ABSUNBAL less than one
c and corrected if necessary afterwards.
*
      DO I = 1, LEN2
         DX = DFLOAT(I - LEN1)
         ALPHA = 1.0D0
         SN = (0.5D0 - SINCS)*PTSEP
         SUM = 0.0D0
c ALPHA is a scaling factor for each of the delta functions in the convolutio
c function, depending on the ratio of amplitudes of the two beams.
         DO J = 1, NSINCS2
            IF (J .EQ. NSINCS+1) THEN
               ALPHA = 1.0D0
            ELSE
               ALPHA = ALPHA*ABSUNBAL
            ENDIF
            SN = SN + PTSEP
            PHI_IJ = (DX-SN)*PI
c                      PHI_IJ is the phase angle at point i relative to SINC8 j
            SUM = SUM - DSIGN(ALPHA,UNBAL)*SINC8(PHI_IJ)*DSIGN(1.0D0,SN)
         ENDDO
         RFUN8(I) = SUM
      ENDDO
*
*      open(unit=7,file="restore1.dat")
*      do i=1,len2
*         write(7,*)i,rfun8(i)
*      enddo
*      close(unit=7)
      IF (DABS(ABSUNBAL - 1.0D0) .LT. SMALLIM) RETURN
c Now normalise the whole function according to the values of the two centre
c delta functions.
      ANORM = 2.0D0 / (1.0D0 + ABSUNBAL**(NSINCS+1))
      DO I = 1, LEN2
         RFUN8(I) = RFUN8(I)*ANORM
      ENDDO
*      open(unit=7,file="restore2.dat")
*      do i=1,len2
*         write(7,*)i,rfun8(i)
*      enddo
*      close(unit=7)
      IF (DABS(UNBAL) .LT. 1.0D0) RETURN
*
c Now correct for the case where UNBAL was greater than one, by reversing
c the order of the array RFUN8, and changing its sign.
c Note that we only swop points up to exactly the middle polnt.
*
      DO I = 1, ICENTR
         J = LEN2 + 1 - I
         X = RFUN8(J)
         RFUN8(J) = -RFUN8(I)
         RFUN8(I) = -X
      ENDDO
*      open(unit=7,file="restore3.dat")
*      do i=1,len2
*         write(7,*)i,rfun8(i)
*      enddo
*      close(unit=7)
*
      RETURN
      END
*
      SUBROUTINE CONV8(SCAN8,CSCAN8,LEN1,CFUN8,LEN2,ICENTR,DNORM,
     $DUMMYMIN,DUMMYMAX,MASK,MASK_LOGI)
*
      INCLUDE 'parameter.inc'
*
      INTEGER*4 LEN1, LEN2, ICENTR
      REAL*8 SCAN8(LEN1), CSCAN8(LEN1), CFUN8(LEN2), DNORM
      REAL*8 DUMMYMIN,DUMMYMAX
      INTEGER MASK(1)
      LOGICAL MASK_LOGI                  ! consider mask or not
* Local variables
      REAL*8  SUM
      INTEGER I,J,I1,DI,NSUM
      LOGICAL TEST_BLANKING
*
c CONV8 convolves vector SCAN8 with vector CFUN8, putting the result into
c vector CSCAN8. LEN1 and LEN2 are the sizes of SCAN8 and CFUN8; the centre
c element of CSCAN8 is at element ICENTR. The convolution is normalized by
c dividing the answer in CSCAN8 by the factor DNORM.
c real*8 version 09.06.1993, Robert Zylka MPIfR Bonn
*
*
      TEST_BLANKING = .FALSE.
      DO I = 1,LEN1
         IF (SCAN8(I).LT.DUMMYMAX.AND.SCAN8(I).GT.DUMMYMIN) THEN
            TEST_BLANKING = .TRUE.
         ENDIF
      ENDDO
* DB:03-99 :  if (.NOT.MASK_LOGI.and.not.test_blanking) : to save time if
*             no mask and no blanking value
      IF (.NOT.MASK_LOGI.AND..NOT.TEST_BLANKING) THEN!if no mask & no blanking
         DO I = 1, LEN1
            SUM = 0
            NSUM = 0
            DI = ICENTR -I               ! 0 <= di <= len1-1
            DO J = 1+DI, LEN1+DI
               I1 = J-DI
               SUM = SUM + SCAN8(I1)*CFUN8(J)
               NSUM = NSUM + 1
            ENDDO
            CSCAN8(I) = SUM / DNORM
            IF (NSUM.EQ.0) CSCAN8(I)=D_BLANKING
         ENDDO
*
      ELSEIF  (MASK_LOGI.AND..NOT.TEST_BLANKING) THEN! test mask
         DO I = 1, LEN1
            SUM = 0
            NSUM = 0
            DI = ICENTR -I
            DO J = 1+DI, LEN1+DI
               I1 = J-DI
               IF (MASK(I1).EQ.0) THEN   ! if inside polygon
                  SUM = SUM + SCAN8(I1)*CFUN8(J)
                  NSUM = NSUM + 1
               ENDIF
            ENDDO
            CSCAN8(I) = SUM / DNORM
            IF (NSUM.EQ.0) CSCAN8(I)=D_BLANKING
         ENDDO
*
      ELSE                               ! test mask and blanking value
         DO I = 1, LEN1
            SUM = 0
            NSUM = 0
            DI = ICENTR -I
            DO J = 1+DI, LEN1+DI
               I1 = J-DI
               IF (SCAN8(I1).GT.DUMMYMAX.OR.SCAN8(I1).LT.DUMMYMIN) THEN
                  IF (MASK(I1).EQ.0) THEN! if inside polygon
                     SUM = SUM + SCAN8(I1)*CFUN8(J)
                     NSUM = NSUM + 1
                  ENDIF
               ENDIF
            ENDDO
            CSCAN8(I) = SUM / DNORM
            IF (NSUM.EQ.0) CSCAN8(I)=D_BLANKING
         ENDDO
      ENDIF
      RETURN
      END
*
      REAL*8 FUNCTION SINC8(PHI)
      IMPLICIT NONE
      REAL*8 PHI, SMALLIM
c real*8 version 07.02.1992, Robert Zylka MPIfR Bonn
*
      SMALLIM = 1.0D-35
*
      IF (DABS(PHI) .LT. SMALLIM) THEN
         SINC8 = 1.0D0
      ELSE
         SINC8 = DSIN(PHI) / PHI
      ENDIF
      RETURN
      END
*
*
