* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      PROGRAM RESTORE_TASK
*
*
* ==============
*
* This is the main program of the gildas task RESTORE.
*
*
      INCLUDE 'inc:memory.inc'
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'const.inc'
*
C Shared memory parameters
      INTEGER*4  POINTER
      INTEGER*4  IADDR_DATA,IADDR        ! address of shared memory area
      INTEGER*4  IP_RMAP,IP_DBMAP,IP
      INTEGER*4  IP_CHAN_FLAG
      INTEGER*4  IP_MASK,IP_SUPPORT_MDB,IP_SUPPORT_CHAN
      INTEGER*4  IP_SUPPORT_EKH
      INTEGER*4  IP_Z1,IP_WORK8,IP_SCAN8,IP_CSCAN8
      INTEGER*4  IP_FFUN8,IP_RFUN8
      INTEGER*4  IP_PARAM
      INTEGER*4        NLI,NCO,NRECORD,NSUBSCAN,NSUP_MDB,NSUP_EKH
      INTEGER          MAX_SEG
      PARAMETER        (MAX_SEG=100)
*
C General parameters
      INTEGER          NCHAN             ! number of channels (7)
      INTEGER          SCAN_NUMBER       ! filename = scan number
      REAL*4           MAXMIN(8)
      CHARACTER        SNAM*40           ! source name
      CHARACTER        FILE*40
      INTEGER          BLC(4), TRC(4)
*
      INTEGER          SIC_GETVM
      INTEGER          IER
      INTEGER          I
      INTEGER          DATA_SHMID
      INTEGER          N_RECEIVERS,NPARAM_RECEIVER,NPARAM_SUBSCAN
*
C Parameters for restore
*
      REAL*8 DUMMYMIN, DUMMYMAX,
     $UNITFACT1, XINC1, YINC1, BMSEP, THETADEG, UNBAL
      INTEGER*4 NODIN, LHEAD1, NGRID1, NCOLS1,
     $NROWS1, ISTART, ISTAT
*
      CHARACTER*80 STRING
*
      LOGICAL*4    ERROR,FLAG,BLANK_DATA
      REAL*4       VAL,XVAL,YVAL,XREF,YREF
      LOGICAL      RESTORE_MASK, RESTORE_FILTER
*
      INTEGER*4    ICHAN,ICHAR0,OFFMAP
      INTEGER      NPIXEL
      INTEGER      COMB_TEETH
      DATA   BLC /0,0,0,0/, TRC /0,0,0,0/
* initialisation.
      BLANK_DATA = .TRUE.
*
* Read gildas parameters
*
      CALL GILDAS_OPEN
      CALL GILDAS_INTE('DATA_SHMID$',DATA_SHMID,1)
      CALL GILDAS_INTE('NRECORD$',NRECORD,1)
      CALL GILDAS_INTE('NSUBSCAN$',NSUBSCAN,1)
      CALL GILDAS_INTE('SCAN_NUMBER$',SCAN_NUMBER,1)
      CALL GILDAS_CHAR('SNAM$',SNAM)
      CALL GILDAS_DBLE('BEAMSEP$',BMSEP,1)
*
      CALL GILDAS_DBLE('THETADEG$',THETADEG,1)
      CALL GILDAS_DBLE('UNBAL$',UNBAL,1)
      CALL GILDAS_LOGI('RESTORE_FILTER$',RESTORE_FILTER,1)
      CALL GILDAS_LOGI('RESTORE_MASK$',RESTORE_MASK,1)
      CALL GILDAS_INTE('COMB_TEETH$',COMB_TEETH,1)
      CALL GILDAS_CLOSE
*
* Read file xxxx.mdb
      FILE = ' '
      WRITE (FILE(1:4),'(i4.4)') SCAN_NUMBER
      CALL SIC_PARSEF(FILE,Y_FILE,' ','.mdb')
      CALL GDF_GEIS  (Y_ISLO,ERROR)      ! reserve slot
      CALL GDF_REIS  (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
      CALL GDF_READY  (Y_ISLO,ERROR)
      CALL GDF_GEMS  (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
      IP_DBMAP = POINTER (Y_ADDR,MEMORY)
      NCO    = Y_DIM(1)
      NLI    = Y_DIM(2)
      NCHAN  = Y_DIM(3)
*
C+WIN32
c      CALL OPEN_SHARED_MEMORY(DATA_SHMID)
C-WIN32
      CALL ATTACH_MEMORY(DATA_SHMID,IADDR)
      IP = POINTER(IADDR,MEMORY)
      CALL R4TOR4(MEMORY(IP),N_RECEIVERS,1)
      CALL R4TOR4(MEMORY(IP+1),NPARAM_RECEIVER,1)
      CALL R4TOR4(MEMORY(IP+2),NPARAM_SUBSCAN,1)
      IP_PARAM = IP +3
      IP_SUPPORT_MDB = IP_PARAM +
     $NPARAM_RECEIVER*NCHAN + NSUBSCAN*NPARAM_SUBSCAN
     $+ (5+N_RECEIVERS)*NRECORD
      CALL R4TOR4(MEMORY(IP_SUPPORT_MDB),NSUP_MDB,1) ! mdb support
      IP_SUPPORT_EKH = IP_SUPPORT_MDB +  4*NCHAN*MAX_SEG+1
      CALL R4TOR4(MEMORY(IP_SUPPORT_EKH),NSUP_EKH,1)
      IF (RESTORE_MASK) THEN
         IF (NSUP_MDB.EQ.0) THEN
            WRITE(6,'(/,a)')
     $      '*** Mask not available => restore without mask  ***'
            RESTORE_MASK = .FALSE.
         ENDIF
      ENDIF
      IF (RESTORE_MASK .AND. RESTORE_FILTER) THEN
         WRITE(6,'(/,a,a,/,a,a)')
     $   '*** If you are trying to mask a strong source,',
     $   'you''d better remove filter ***',
     $   '*** Use the following command :',
     $   ' let restore_filter no ***'
      ENDIF
*
      NGRID1 = NLI*NCO
      NCOLS1 = NCO
      NROWS1 = NLI
      LHEAD1 = 0
*
      XINC1 = Y_INC1/SEC_TO_RAD          ! arcseconds
      YINC1 = Y_INC2/SEC_TO_RAD
      XVAL  = Y_VAL1/SEC_TO_RAD
      YVAL  = Y_VAL2/SEC_TO_RAD
      XREF  = Y_REF1
      YREF  = Y_REF2
*
      ISTART = 0
*
      DUMMYMIN = D_BLANKING - 0.5D0
      DUMMYMAX = D_BLANKING + 0.5D0
*
      UNITFACT1 = 1                      ! used by LRMAP8 (rotation of the map)
* 1,       if sint in arseconds
* 3600,    if sint in degrees
      NODIN = 0                          ! nod2 information (if > 0)
*
* allocate memory for working, for results and for masking area
      IER = SIC_GETVM(NCO*NLI*(NCHAN+5)+12*NCO,IADDR_DATA)
      IF (IER.NE.1) THEN
         WRITE(6,*) ' Unable to allocate memory'
         STOP
      ENDIF
      IP_Z1    = POINTER(IADDR_DATA,MEMORY)
      IP_WORK8 = IP_Z1    + 2*NCO*NLI
      IP_SCAN8 = IP_WORK8 + 2*NCO*NLI
      IP_CSCAN8= IP_SCAN8 + 2*NCO
      IP_FFUN8 = IP_CSCAN8+ 2*NCO
      IP_RFUN8 = IP_FFUN8 + 4*NCO
*
      IP_RMAP  = IP_RFUN8 + 4*NCO
      IP_MASK  = IP_RMAP  + NCO*NLI*NCHAN
*
* Channel flags
      IP_CHAN_FLAG = IP_PARAM+2*NCHAN
*
* initialization of image
      ICHAR0 = ICHAR('0')
      DO ICHAN=1,NCHAN
         CALL R4TOR4(MEMORY(IP_CHAN_FLAG+ICHAN-1),FLAG,1)
         IF (.NOT.FLAG) THEN
            IF (RESTORE_MASK) THEN
               IP_SUPPORT_CHAN = IP_SUPPORT_MDB +
     $         2*NSUP_MDB*(ICHAN-1)+1
               CALL INIT_MASK (MEMORY(IP_SUPPORT_CHAN),
     $         MEMORY(IP_SUPPORT_CHAN +NSUP_MDB),
     $         NSUP_MDB,MEMORY(IP_MASK),NCO,NLI,
     $         XVAL,XINC1,YVAL,YINC1,XREF,YREF)
            ENDIF
            CALL R4TOR8(MEMORY(IP_DBMAP+(ICHAN-1)*NCO*NLI),
     $      MEMORY(IP_Z1),NCO*NLI)
***
            CALL  RESMAP8(NODIN,LHEAD1,NGRID1,NCOLS1,NROWS1,DUMMYMIN,
     $      DUMMYMAX,UNITFACT1,XINC1,YINC1,BMSEP,
     $      THETADEG,UNBAL,STRING,ISTART,ISTAT,MEMORY(IP_MASK),
     $      RESTORE_MASK,RESTORE_FILTER,
     $      MEMORY(IP_Z1),MEMORY(IP_WORK8),MEMORY(IP_SCAN8),
     $      MEMORY(IP_CSCAN8),MEMORY(IP_FFUN8),MEMORY(IP_RFUN8),
     $      COMB_TEETH)
***
            CALL R8TOR4(MEMORY(IP_Z1),
     $      MEMORY(IP_RMAP+(ICHAN-1)*NCO*NLI),NCO*NLI)
            IF (BLANK_DATA) THEN
               OFFMAP = (ICHAN-1)*NCO*NLI
               DO I=1,NCO*NLI
                  CALL R4TOR4(MEMORY(IP_DBMAP+OFFMAP+I-1),VAL,1)
                  IF (VAL.EQ.BLANKING)
     $            CALL R4TOR4(BLANKING,MEMORY(IP_RMAP+OFFMAP+I-1),1)
               ENDDO
            ENDIF
*
            IF (RESTORE_MASK) THEN
               IP_SUPPORT_CHAN = IP_SUPPORT_EKH +
     $         2*NSUP_EKH*(ICHAN-1)+1
               CALL INIT_MASK (MEMORY(IP_SUPPORT_CHAN),
     $         MEMORY(IP_SUPPORT_CHAN +NSUP_EKH),
     $         NSUP_EKH,MEMORY(IP_MASK),NCO,NLI,
     $         XVAL,XINC1,YVAL,YINC1,XREF,YREF)
               CALL MASK_APPLY(MEMORY(IP_RMAP+(ICHAN-1)*NCO*NLI),
     $         BLANKING,MEMORY(IP_MASK),NCO*NLI)
            ENDIF
**DB:03-96   mask beamsep/2.
            NPIXEL = BMSEP/XINC1/2.
            CALL MASK_EDGES(MEMORY(IP_RMAP+(ICHAN-1)*NCO*NLI),
     $      BLANKING,NCO,NLI,NPIXEL)
**DB:03-96
         ELSE
            DO I=1,NCO*NLI
               CALL R4TOR4(BLANKING,
     $         MEMORY(IP_RMAP+(ICHAN-1)*NCO*NLI+I-1),1)
            ENDDO
         ENDIF
      ENDDO
*
C --- make GILDAS image file
*
      CALL MAXMIN_IMAGE (MEMORY(IP_RMAP),NCO,NLI,NCHAN,MAXMIN)
      X_NDIM   = 3
      X_UNIT   = Y_UNIT
      X_DIM(1) = NCO
      X_DIM(2) = NLI
      X_DIM(3) = NCHAN
      X_DIM(4) = 1
      X_SIZE   = X_DIM(1)*X_DIM(2)*X_DIM(3)
*
      X_REF1 = 1
      X_VAL1 = Y_VAL1
      X_INC1 = Y_INC1
*
      X_REF2 = 1
      X_VAL2 = Y_VAL2
      X_INC2 = Y_INC2
*
      X_REF3 = 0
      X_VAL3 = 0
      X_INC3 = 1
*
      X_TYPE = 'GILDAS_IMAGE'
      X_NAME = SNAM
      X_FORM = FMT_R4
      X_GENE = 29
      X_BLAN = 2
      X_BVAL = BLANKING
      X_EVAL = 0.
      X_EXTR = 10
      X_RMIN = MAXMIN(1)
      X_RMAX = MAXMIN(2)
      X_MIN1 = MAXMIN(3)
      X_MIN2 = MAXMIN(4)
      X_MIN3 = MAXMIN(5)
      X_MAX1 = MAXMIN(6)
      X_MAX2 = MAXMIN(7)
      X_MAX3 = MAXMIN(8)
      X_MIN4 = 1
      X_MAX4 = 1
*
      X_DESC = 18
      X_POSI = 12
      X_RA   = Y_RA
      X_DEC  = Y_DEC
      X_PROJ = 9
      X_SPEC = 12
      X_RESO = 3
      X_SYST = 'AZIMUTHAL'
      X_FREQ = Y_FREQ
*
      SNAM = ' '
      WRITE (SNAM(1:4),'(i4.4)') SCAN_NUMBER
      CALL SIC_PARSEF(SNAM,X_FILE,' ','.ekh')
      CALL GDF_GEIS  (X_ISLO,ERROR)      ! reserve slot
      CALL GDF_WRITX (X_ISLO,ERROR)
      CALL GDF_CRIS  (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
      CALL GDF_GEMS  (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
      IP = POINTER (X_ADDR,MEMORY)
      CALL R4TOR4(MEMORY(IP_RMAP),MEMORY(IP),NCO*NLI*NCHAN)
      CALL GDF_FRIS (X_ISLO,ERROR)
      CALL GDF_FRIS (Y_ISLO,ERROR)
      END
*
*
c --------------------------------------------------------------------------
*
*
      SUBROUTINE INIT_MASK(XU,YU,PNT,MASK_MAP,NX,NY,
     $VAL1,INC1,VAL2,INC2,REF1,REF2)
*
      INCLUDE 'parameter.inc'
      INTEGER NX,NY,PNT
      INTEGER MASK_MAP(NX,NY)
      REAL*4  XU(1),YU(1)
      REAL*8  INC1,INC2
      REAL*4  VAL1,VAL2,REF1,REF2
*
      INTEGER MAX_SEG
      PARAMETER (MAX_SEG = 100)
      REAL*8  GONS(MAX_SEG,4,2), XCOO, YCOO, BOUND(5,2)
      INTEGER I,J,NPOLY,PNT_POLY(2),IPOL,II
      LOGICAL IN,GR8_IN
*
* Initialization
      DO I=1,NX
         DO J=1,NY
            MASK_MAP(I,J) = 0
         ENDDO
      ENDDO
*
* Test if any support defined
      IF (ABS(XU(1)-BLANKING).LT.1) RETURN
*
* Test if 1 or 2 polygons
      NPOLY = 1
      PNT_POLY(1) = PNT
      DO I=2,PNT-1
         IF ((ABS(XU(I)-XU(1)).LT.0.001).AND.
     $   (ABS(YU(I)-YU(1)).LT.0.001)) THEN
            NPOLY = 2
            PNT_POLY(1) = I
            GOTO 100
         ENDIF
      ENDDO
100   PNT_POLY(2) = PNT - PNT_POLY(1)
*
      DO IPOL=1,NPOLY
         DO I=1,PNT_POLY(IPOL)
            II = I
            IF (IPOL.EQ.2) II = I+PNT_POLY(1)
            GONS(I,1,IPOL) = XU(II)
            GONS(I,2,IPOL) = YU(II)
            IF (I.LT.PNT_POLY(IPOL)) THEN
               GONS(I,3,IPOL) = XU(II+1)-XU(II)
               GONS(I,4,IPOL) = YU(II+1)-YU(II)
            ENDIF
         ENDDO
*     Initialize BOUND (for GR8_IN)
*     bound(1) don't need to be initialized
         BOUND(2,IPOL)= GONS(1,1,IPOL)
         BOUND(3,IPOL)= GONS(1,1,IPOL)
         BOUND(4,IPOL)= GONS(1,2,IPOL)
         BOUND(5,IPOL)= GONS(1,2,IPOL)
         DO I=2,PNT_POLY(IPOL)
            BOUND(2,IPOL) = MIN(BOUND(2,IPOL),GONS(I,1,IPOL))    ! xmin
            BOUND(3,IPOL) = MAX(BOUND(3,IPOL),GONS(I,1,IPOL))    ! xmax
            BOUND(4,IPOL) = MIN(BOUND(4,IPOL),GONS(I,2,IPOL))    ! ymin
            BOUND(5,IPOL) = MAX(BOUND(5,IPOL),GONS(I,2,IPOL))    ! ymax
         ENDDO
      ENDDO
*
      DO I=1,NX
         DO J=1,NY
            XCOO = VAL1+INC1*(I-REF1)
            YCOO = VAL2+INC2*(J-REF2)
            IN = GR8_IN (XCOO,YCOO,MAX_SEG,PNT_POLY(1)-1,GONS(1,1,1),
     $      BOUND(1,1))
            IF (.NOT.IN .AND. NPOLY.EQ.2)
     $      IN = GR8_IN (XCOO,YCOO,MAX_SEG,PNT_POLY(2)-1,GONS(1,1,2),
     $      BOUND(1,2))
            IF (.NOT.IN) MASK_MAP(I,J) = 1
         ENDDO
      ENDDO
      RETURN
      END
c --------------------------------------------------------------------------
*
*
      SUBROUTINE MASK_APPLY(DATA,BLANKING,MASK,NPTS)
*
      INTEGER MASK(1),NPTS
      REAL    DATA(1),BLANKING
*
      INTEGER I
*
      DO I=1,NPTS
         IF (MASK(I).EQ.1) DATA(I) = BLANKING
      ENDDO
      END
c --------------------------------------------------------------------------
*
*
      SUBROUTINE MASK_EDGES(DATA,BLANKING,NCO,NLI,NPTS)
*
      INTEGER NCO,NLI,NPTS
      REAL    DATA(NCO,NLI), BLANKING
*
      INTEGER I,J
*
      DO J=1,NLI
         DO I=1,NPTS
            DATA(I,J) = BLANKING
         ENDDO
         DO I=NCO-NPTS+1,NCO
            DATA(I,J) = BLANKING
         ENDDO
      ENDDO
      END