* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* combine.f
*
      SUBROUTINE  COMBINE(SCAN,RADEC_DATA,RADEC_WEIGHT,IFUN,ERROR,
     $MAP_UNIT)
*
* input:
*  - scan()         : list of scan numbers to be combined
*  - radec_data()   : data*weight
*  - radec_weight() : weight
* output:
*  - radec_data()   : data normalized by weight
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:memory.inc'
      INCLUDE 'const.inc'
      INCLUDE 'convert.inc'
      INCLUDE 'parameter.inc'
*
      INTEGER SCAN(1)
      REAL*4  RADEC_DATA(1),RADEC_WEIGHT(1),FACTOR,FACTOR_WEIGHT
*
      REAL*8     W_MIN,W_MAX
      REAL*8     ALOW,ALIM,DIAMETER
*
C Shared memory parameters
      INTEGER*4  POINTER
      INTEGER*4  IPMAP,IPWEIGHT
*
      INTEGER     ISCAN
      INTEGER     NX,NY
      REAL*4      VAL
*
      LOGICAL   ERROR
      CHARACTER*40 FILE
      INTEGER   BLC(4), TRC(4)
*
*
      REAL*8  XINV,YINV
      REAL*8  FUNWX,FUNWY,XINVFX,XINVFY,TOTWX,TOTWY
      REAL*8  XLIMM,YLIMM
      REAL*8  DRA,DDEC
      REAL*8  XPTS,YPTS
      REAL*8  SOFFX,SOFFY
      REAL*8  X,Z,ZZ,MAGI
      REAL*8  RA_CENTER,DEC_CENTER,RA,DEC,DELTA_RA,DELTA_DEC
      REAL*8  MAP_PIXEL
      REAL*4  POINT_WEIGHT
*
      INTEGER NPTS,MPTS,IFUN
      INTEGER I,J,JJ,ICONO,ICO,KK
      INTEGER IFIRST_X, ILAST_X, IFIRST_Y, ILAST_Y
      CHARACTER*12 MAP_UNIT
*
      DATA BLC /0,0,0,0/, TRC /0,0,0,0/
      ERROR = .FALSE.
*
* Initializations
      XINV=1.D0/(XINC_RADEC/3600.)       ! unit=1/degrees
      YINV=1.D0/(YINC_RADEC/3600.)
*
      MAGI = 0
      IF (IFUN.EQ.1) MAGI = 1.6806492007530D0  ! Parabolic
      IF (IFUN.EQ.2) MAGI = 1.6982249549743D0  ! Gaussian
      IF (IFUN.EQ.3) MAGI = 1.8363197305085D0  ! Linear
      IF (MAGI.EQ.0) MAGI = 1.6574002064571D0  ! No taper
*
      DIAMETER = 30                      ! pico veleta
      IF (TELESCOPE(1:3).EQ.'CSO') DIAMETER = 10
      IF (TELESCOPE(1:3).EQ.'HHT') DIAMETER = 10
      FUNWX = VLIGHT/FREQUENCY /(MAGI*DIAMETER)
     $/(XINC_RADEC*SEC_TO_RAD)           ! lambda/(1.67*diam_telescope)/xinc
      FUNWY = VLIGHT/FREQUENCY /(MAGI*DIAMETER)
     $/(YINC_RADEC*SEC_TO_RAD)           ! lambda/(1.67*diam_telescope)/yinc
*
      IF(FUNWX.LT.1.D0)FUNWX=1.D0
      IF(FUNWY.LT.1.D0)FUNWY=1.D0
      XINVFX =1.D0/FUNWX
      XINVFY =1.D0/FUNWY
      TOTWX  = 4.0D0 * FUNWX
      TOTWY  = 4.0D0 * FUNWY
      XLIMM = NX_RADEC-1.D0+TOTWX
      YLIMM = NY_RADEC-1.D0+TOTWY
*
      NPTS = NX_RADEC
      MPTS = NY_RADEC
*
      DO I=1, NX_RADEC*NY_RADEC
         RADEC_DATA(I) = 0
         RADEC_WEIGHT(I) = 0
      ENDDO
*
*  Read each scan and combine them in the new grid
*
      DO ISCAN =1,NSCAN
         SCAN_NUMBER = SCAN(ISCAN)
*
*      Read file "scan_number".cnv
         FILE = ' '
         WRITE (FILE(1:4),'(i4.4)') SCAN_NUMBER
         CALL SIC_PARSEF(FILE,Y_FILE,' ','.cnv')
         CALL GDF_GEIS  (Y_ISLO,ERROR)   ! reserve slot
         CALL GDF_REIS  (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
         IF (ERROR) THEN
            WRITE(6,*) ' error on reading file :',Y_FILE
            RETURN
         ENDIF
         CALL GDF_READY  (Y_ISLO,ERROR)
         CALL GDF_GEMS  (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
         IPMAP = POINTER (Y_ADDR,MEMORY)
         IPWEIGHT = IPMAP + Y_DIM(1)*Y_DIM(2)
         NX = Y_DIM(1)
         NY = Y_DIM(2)
         WRITE(6,'(a,i4.4,a,a)') '   Reading map ',SCAN_NUMBER,
     $   '.cnv ; Map unit is ', Y_UNIT
         FACTOR = 1.
         IF (ISCAN.EQ.1) THEN
            MAP_UNIT = Y_UNIT
            MAP_PIXEL = Y_INC1
         ELSE
            IF (Y_UNIT.NE.MAP_UNIT) THEN
               IF (INDEX(MAP_UNIT,'mJy').NE.0
     $         .AND.INDEX(Y_UNIT,' Jy').NE.0) THEN
                  FACTOR = 1000.
               ELSEIF (INDEX(MAP_UNIT,' Jy').NE.0
     $         .AND.INDEX(Y_UNIT,'mJy').NE.0) THEN
                  FACTOR = 0.001
               ELSE
                  WRITE(6,'(/,2a)')
     $            ' *** Scans have different units. ',
     $            ' Combine is aborted *** '
                  STOP
               ENDIF
            ENDIF
            IF (ABS(Y_INC1/MAP_PIXEL -1.).GT.1.D-06) THEN
               WRITE(6,'(/,a,/,a,a)')
     $         ' *** WARNING: Your maps have different pixel size.',
     $         ' *** You may have to adjust weight maps in a ',
     $         ' convenient way.'
            ENDIF
         ENDIF
*         FACTOR_WEIGHT = 1/FACTOR**2
         FACTOR_WEIGHT = 1.
         RA  = Y_A0 / DEG_TO_RAD
         DEC = Y_D0 / DEG_TO_RAD
*
* Calculate delta_ra and delta_dec between 2 maps
         RA_CENTER  = SLAM               ! DELCEN needs REAL*8
         DEC_CENTER = SBET
         CALL DELCEN(RA_CENTER,DEC_CENTER,RA,DEC,DELTA_RA,DELTA_DEC)
*
* Compute rms
c         CALL COMPUTE_RMS_MAP(MEMORY(IPMAP),RMS,NX*NY)
c         IF (RMS.EQ.0) RETURN
c         MAP_WEIGHT = 1.E6/RMS
*
* Add scan data to radec_data and radec_weight
         DO J = 1, NY
            DO I = 1, NX
               CALL R4TOR4(MEMORY(IPMAP+(J-1)*NX+I-1),VAL,1)
               CALL R4TOR4(MEMORY(IPWEIGHT+(J-1)*NX+I-1),POINT_WEIGHT,1)
               POINT_WEIGHT = POINT_WEIGHT*FACTOR_WEIGHT
               IF (ABS(VAL-BLANKING).GT.0.5) THEN
                  DRA =  Y_VAL1 + (I-Y_REF1)*Y_INC1
     $            - DELTA_RA*DEG_TO_RAD  ! radians
                  DDEC = Y_VAL2 + (J-Y_REF2)*Y_INC2
     $            - DELTA_DEC*DEG_TO_RAD ! radians
                  XPTS = (DRA - X_VAL1) / X_INC1 + X_REF1
                  YPTS = (DDEC- X_VAL2) / X_INC2 + X_REF2
                  IF(XPTS.LE.-TOTWX.OR.XPTS.GE.XLIMM) GO TO 8
                  IF(YPTS.LE.-TOTWY.OR.YPTS.GE.YLIMM) GO TO 8
                  SOFFY = YPTS
*
                  IFIRST_Y = SOFFY - TOTWY + 1
                  ILAST_Y  = SOFFY + TOTWY
                  IF (IFIRST_Y.LT.1)    IFIRST_Y = 1
                  IF (ILAST_Y .GT.MPTS) ILAST_Y  = MPTS
*
                  DO JJ = IFIRST_Y, ILAST_Y
                     X = (DFLOAT(JJ) - SOFFY) * XINVFY
                     X = DABS(X)
                     ICONO = INT(5D1*X) + 1
                     Z = FUNC(ICONO)*POINT_WEIGHT
c                    write(2,*)j,i,icono, FUNC(ICONO)
c                    write(3,*) IFIRST_Y, ILAST_Y,map_weight,
c     $                    point_weight,z,x,icono,func(icono)
                     ICO = (JJ-1) * NPTS ! index of beginning of line
                     SOFFX = XPTS
                     IFIRST_X = SOFFX - TOTWX + 1
                     ILAST_X  = SOFFX + TOTWX
                     IF (IFIRST_X.LT.1)    IFIRST_X = 1
                     IF (ILAST_X .GT.NPTS) ILAST_X  = NPTS
                     DO  KK = IFIRST_X, ILAST_X
                        X = (DFLOAT(KK)-SOFFX) * XINVFX
                        X = DABS(X)
                        ICONO = INT(5D1*X) + 1
                        ZZ = Z * FUNC(ICONO)
*                           ZZ=y_weight* x_weight*point_weight
*
                        RADEC_DATA(ICO+KK) = RADEC_DATA(ICO+KK) +
     $                  VAL*ZZ*FACTOR
                        RADEC_WEIGHT(ICO+KK) = RADEC_WEIGHT(ICO+KK)
     $                  + ZZ
                     ENDDO
                  ENDDO
               ENDIF
8              CONTINUE
            ENDDO
         ENDDO
         CALL GDF_FRIS(Y_ISLO,ERROR)
      ENDDO
*
*
* --- Compute maximum and minimum weight
*
      CALL WMAX(RADEC_DATA,RADEC_WEIGHT,W_MAX,W_MIN,BLANKING,
     $NX_RADEC*NY_RADEC)
      IF (W_MAX.LT.1.D-08) THEN
         WRITE(6,*) ' '
         WRITE(6,*) ' *** No point in the selected field. ',
     $   'Convert task is aborted*** '
         ERROR = .TRUE.
         RETURN
      ENDIF
*
*
* --- Normalize :
*        divide by weight
*        flagge data if weight<0.5*w_min or weight<0.01*w_max
      ALIM = CUTOFF*W_MAX                ! original value: 0.2*w_max
      ALOW = 0.5*W_MIN                   ! useful if W_MIN negative ???
      CALL MNORM(RADEC_DATA,RADEC_WEIGHT,ALIM,ALOW,BLANKING,
     $NX_RADEC*NY_RADEC,1.)
      CALL NORM_WEIGHT(RADEC_WEIGHT,NX_RADEC*NY_RADEC,IFUN,XINC_RADEC,
     $YINC_RADEC,TELESCOPE)
*
      RETURN
      END
*
*
      SUBROUTINE COMPUTE_RMS_MAP(X,RMS,N)
*
      INCLUDE 'parameter.inc'
      REAL    X(1),RMS
      INTEGER N,NP,I
      REAL MEAN
*
      MEAN = 0
      NP = 0
      DO I=1,N
         IF (ABS(X(I)-BLANKING).GT.0.5) THEN
            MEAN = MEAN + X(I)
            NP = NP+1
         ENDIF
      ENDDO
      IF (NP.EQ.0) RETURN
      MEAN = MEAN/NP
*
      RMS = 0
      DO I=1,N
         IF (ABS(X(I)-BLANKING).GT.0.5) THEN
            RMS = RMS + (X(I)-MEAN)**2
         ENDIF
      ENDDO
      RMS = SQRT(RMS/NP)
      RETURN
      END
*
*
      SUBROUTINE COMBINE_LIMITS(SCAN)
*-------------------------------------------------------------------
*
*  Find automatic limits for combined maps
*
*-------------------------------------------------------------------
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:memory.inc'
      INCLUDE 'const.inc'
      INCLUDE 'convert.inc'
      INCLUDE 'parameter.inc'
*
      INTEGER   SCAN(1)
*
      REAL      YMAX
      REAL      X1,X2,Y1,Y2
      LOGICAL   ERROR
      CHARACTER*40 FILE
      INTEGER   ISCAN
      REAL      DX,DY
      INTEGER   IX,IY
      INTEGER   BLC(4), TRC(4)
      DATA BLC /0,0,0,0/, TRC /0,0,0,0/
*
      XMIN =  1.E32
      XMAX = -1.E32
      YMIN =  1.E32
      YMAX = -1.E32
      DO ISCAN =1,NSCAN
         SCAN_NUMBER = SCAN(ISCAN)
*
*       Read file "scan_number".cnv
         FILE = ' '
         WRITE (FILE(1:4),'(i4.4)') SCAN_NUMBER
         CALL SIC_PARSEF(FILE,Y_FILE,' ','.cnv')
         CALL GDF_GEIS  (Y_ISLO,ERROR)   ! reserve slot
         CALL GDF_REIS  (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
         IF (ERROR) THEN
            WRITE(6,*) ' error on reading file :',Y_FILE
            RETURN
         ENDIF
         CALL GDF_READY  (Y_ISLO,ERROR)
         CALL GDF_GEMS  (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
         X1 = Y_VAL1 + (1 - Y_REF1)*Y_INC1 + Y_A0
         X2 = Y_VAL1 + (Y_DIM(1) - Y_REF1)*Y_INC1 + Y_A0
         Y1 = Y_VAL2 + (1 - Y_REF2)*Y_INC2 + Y_D0
         Y2 = Y_VAL2 + (Y_DIM(2) - Y_REF2)*Y_INC2 + Y_D0
         IF (X1.LT.XMIN) XMIN = X1
         IF (X2.LT.XMIN) XMIN = X2
         IF (X1.GT.XMAX) XMAX = X1
         IF (X2.GT.XMAX) XMAX = X2
         IF (Y1.LT.YMIN) YMIN = Y1
         IF (Y2.LT.YMIN) YMIN = Y2
         IF (Y1.GT.YMAX) YMAX = Y1
         IF (Y2.GT.YMAX) YMAX = Y2
         CALL GDF_FRIS(Y_ISLO,ERROR)
      ENDDO
*
      DX  = (XMAX - XMIN)/SEC_TO_RAD/XINC_RADEC/2. + 0.9999
      DY =  (YMAX - YMIN)/SEC_TO_RAD/YINC_RADEC/2. + 0.9999
      IX = DX
      IY = DY
      NX_RADEC = 2*IX + 1
      NY_RADEC = 2*IY + 1
*
      X_VAL1 = 0
      X_INC1 = -XINC_RADEC*SEC_TO_RAD
      X_REF1 = NINT((SLAM*DEG_TO_RAD-(XMAX+XMIN)/2)/X_INC1)
     $+ IX + 1
*
      X_INC2 = YINC_RADEC*SEC_TO_RAD
      X_VAL2 = 0
      X_REF2 = NINT((SBET*DEG_TO_RAD-(YMAX+YMIN)/2)/X_INC2)
     $+ IY + 1
*
      X_A0   = SLAM*DEG_TO_RAD           ! longitude of projection center
      X_D0   = SBET*DEG_TO_RAD           ! latitude of projection center
      END
*
      SUBROUTINE COMBINE_INIT(SCAN)
*-------------------------------------------------------------------
*
*  Find some parameters in first file
*
*-------------------------------------------------------------------
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:memory.inc'
      INCLUDE 'const.inc'
      INCLUDE 'convert.inc'
      INCLUDE 'parameter.inc'
*
      INTEGER   SCAN(1)
*
      LOGICAL   ERROR
      CHARACTER*40 FILE
      INTEGER   BLC(4), TRC(4)
      DATA BLC /0,0,0,0/, TRC /0,0,0,0/
*
      SCAN_NUMBER = SCAN(1)
*
*       Read file "scan_number".cnv
      FILE = ' '
      WRITE (FILE(1:4),'(i4.4)') SCAN_NUMBER
      CALL SIC_PARSEF(FILE,Y_FILE,' ','.cnv')
      CALL GDF_GEIS  (Y_ISLO,ERROR)      ! reserve slot
      CALL GDF_REIS  (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) ' error on reading file :',Y_FILE
         RETURN
      ENDIF
      CALL GDF_READY  (Y_ISLO,ERROR)
      CALL GDF_GEMS  (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
      SNAM = Y_NAME
      SLAM = Y_RA/DEG_TO_RAD
      SBET = Y_DEC/DEG_TO_RAD
      CALL GDF_FRIS(Y_ISLO,ERROR)
      END
