* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* convert.f
*
*
      SUBROUTINE CONVERT(LST,XCOORD,YCOORD,AZEL_DATA,
     $RADEC_DATA,RADEC_WEIGHT,RADEC_DATA_RAW,RECEIVER_FLAG,
     $RECEIVER_RMS,IFUN,ERROR)
*
      INCLUDE 'convert.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'inc:errcod.inc'
*
      REAL*4     LST(1),XCOORD(1),YCOORD(1)
      REAL*4     AZEL_DATA(1)            ! data input (AZ-EL map)
      REAL*4     RADEC_DATA(1),RADEC_WEIGHT(1) ! output
      REAL*4     RADEC_DATA_RAW(1)       ! output
      LOGICAL    RECEIVER_FLAG(1)
      REAL*4     RECEIVER_RMS(1)
*
      REAL*8     W_MIN,W_MAX
      REAL*8     ALOW,ALIM
      INTEGER    ICHAN,I,IFUN
      LOGICAL    ERROR
      ERROR   = .FALSE.
*
* --- Test RECEIVER_RMS : if 0, quit
*
      ERROR = .FALSE.
      DO ICHAN = 1,NCHANNEL
         IF (.NOT.RECEIVER_FLAG(ICHAN)) THEN
            IF (RECEIVER_RMS(ICHAN).EQ.0) THEN
               IF (.NOT.ERROR) WRITE (6,'(/)')
               ERROR = .TRUE.
               WRITE(6,'(2a,i2,a)') '** No receiver weight available',
     $         ' for channel: ',ICHAN,'    CHAN_RMS[] = 0'
            ENDIF
         ENDIF
      ENDDO
      IF (ERROR) THEN
         WRITE (6,'(/2a)')'** Either execute the BASELINE command **'
         WRITE (6,'(2a)')'** or execute the FLAG/CHANNEL command **'
         CALL SYSEXI(FATALE)
      ENDIF
*
*
* --- Raz RADEC map
      DO I=1,NX_RADEC*NY_RADEC
         RADEC_DATA(I)   = 0
         RADEC_WEIGHT(I) = 0
      ENDDO
*
*--- Initialize parameters to compute precession
      CALL SETUP
*
* --- Create RADEC map
*
      NX_AZEL = NCOL                     ! number of columns in AZ-EL map
      NY_AZEL = NROW                     ! number of raws    in AZ-EL map
      DO ICHAN = 1,NCHANNEL
         IF (.NOT.RECEIVER_FLAG(ICHAN)) THEN
            I = 1 + (ICHAN-1)*NX_AZEL*NY_AZEL
            CALL RDCON(AZEL_DATA(I),RADEC_DATA(1),RADEC_WEIGHT(1),
     $      LST(1),XCOORD(I),YCOORD(I),RECEIVER_RMS(ICHAN),IFUN)
         ENDIF
      ENDDO
*
* --- Save raw data
      DO I=1,NX_RADEC*NY_RADEC
         RADEC_DATA_RAW(I) = RADEC_DATA(I)
      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,*) ' *** Maximum weight is close to 0. '
         WRITE(6,'(a,$)')
     $   '  *** It could be that CHAN_RMS was too high or that'
         WRITE(6,*)'selected field was empty.'
         WRITE(6,*)' *** Convert task is aborted. '
         ERROR = .TRUE.
         RETURN
      ENDIF
*
*
* --- Normalize :
*        divide by weight
*        flagge data if weight < 0.01*w_max
*                    or weight < 0.5 *w_min
      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)
*
      END
*
*
      SUBROUTINE RDCON(AZEL_DATA,RADEC_DATA,RADEC_WEIGHT,
     $LST,XCOORD,YCOORD,RECEIVER_RMS,IFUN)
C********************************************************************
C   THIS ROUTINE MAKES THE TRANSFORMATION OF THE AZ-EL DATA INTO    *
C   BEAM-CENTRED RA-DEC COORDINATES AND COMBINES THEM INTO AN       *
C   RA-DEC GRID USING A CONVOLUTION PROCESS.                        *
C********************************************************************
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'convert.inc'
      INCLUDE 'const.inc'
      INCLUDE 'parameter.inc'
*
      REAL*4  AZEL_DATA(1),RADEC_DATA(1),RADEC_WEIGHT(1),
     $LST(1),XCOORD(1),YCOORD(1),RECEIVER_RMS
*
      REAL*8  ARAARA,DELDEL
      REAL*8  FUNWX,FUNWY,XINVFX,XINVFY,TOTWX,TOTWY
      REAL*8  XLIMM,YLIMM
      REAL*8  DRA,DDEC,RA50,DEC50,RA_DAY,DEC_DAY
      REAL*8  RECEIVER_WEIGHT
      REAL*8  DECCUR,OBSLAT
      REAL*8  P1,P2,P3
      REAL*8  XPTS,YPTS
      REAL*8  SOFFX,SOFFY
      REAL*8  X,Z,ZZ, MAGI, DIAMETER
*
      INTEGER IORIG,NPTS,MPTS
      INTEGER IOR, IFUN
      INTEGER I,J,JJ,ICONO,ICO,KK
      INTEGER IFIRST_X, ILAST_X, IFIRST_Y, ILAST_Y
*
*--- Compute RA and DEC differences between 2 maps
c        CALL DELCEN(AMAP,CMAP,ARAARA,DELDEL)
      ARAARA = 0
      DELDEL = 0
*
*
      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
c	      write(6,*)'totwx,xlimm ',totwx,xlimm
c	      write(6,*)'totwy,ylimm ',totwy,ylimm
c	      WRITE(6,*)'FUNWX,FUNWY',FUNWX,FUNWY
*        RECEIVER_WEIGHT = 5. * SNFAC * SNFAC
      RECEIVER_WEIGHT = 1.
      IF (RECEIVER_RMS.NE.0) RECEIVER_WEIGHT = 1./RECEIVER_RMS**2
*
      IORIG = 0
      NPTS = NX_RADEC
      MPTS = NY_RADEC
*
      IOR = 1
*
CC	CMAP(80) = CMAP(80) + AMAP(80)
*
      DRA = 0.D0
      DDEC = 0.D0
      RA50 = SLAM + OLAM
      DEC50 = SBET + OBET
CC      CALL CENTRE(AMAP,RA50,DEC50)
c        PRINT *,'RA 1950 = ',RA50,' DEC OF 1950 = ',DEC50
*
*--- Computes RA_DAY and DEC_DAY
      CALL PREC(3,EPOCHTO,RA_DAY,DEC_DAY,DRA,DDEC,RA50,DEC50)
c      PRINT *,'RA OF DATE = ',RA_DAY,' DEC OF DATE =',DEC_DAY
c
c      call precmat(epochfr,epochto,amat)   ! for test
c      call dcosin(ra50,dec50,avec)         ! for test
c      call matmul(amat,avec,bvec,-9)       ! for test
c      call dangle(ra_day,dec_day,bvec)     ! for test
c      PRINT *,'RA OF DATE = ',RA_DAY,' DEC OF DATE =',DEC_DAY
*
*--- Computes RA50 and DEC50 (input = RA_DAY,DEC_DAY )
      CALL PREC(1,EPOCHTO,RA_DAY,DEC_DAY,DRA,DDEC,RA50,DEC50)
c	TYPE *,'CURRENT COORDS = ',RA_DAY,DEC_DAY,
c	1    ' 1950 COORD = ',RA50, DEC50
*
      DECCUR = DEC_DAY * DEG_TO_RAD
      OBSLAT = SITE_LATITUDE * DEG_TO_RAD
      SINDEC = DSIN(DECCUR)
      COSDEC = DCOS(DECCUR)
      SINOL = DSIN(OBSLAT)
      COSOL = DCOS(OBSLAT)
*
      DO I = 1, NY_AZEL
         DO J = 1, NX_AZEL               ! first and last column not significati
            IOR = J + (I-1)*NX_AZEL
            P1 = XCOORD(IOR)             ! IOR : from 1 to NX_AZEL*NY_AZEL
            P2 = YCOORD(IOR)
            P3 = LST(IOR)
            IF (P3.NE.BLANKING .AND.
     $      AZEL_DATA(IOR).NE.BLANKING) THEN
*---  Computes DRA and DDEC
               CALL RADEC(P1,P2,P3,SITE_LATITUDE,RA_DAY,DEC_DAY,
     $         DRA,DDEC)
               IF (IOR.EQ.1) WRITE(6,*)DRA,DDEC
*---  Computes DRA and DDEC for epoch 1950
               CALL PREC(2,EPOCHTO,RA_DAY,DEC_DAY,DRA,DDEC,RA50,DEC50)
*
               DRA = DRA + ARAARA        ! degrees
               DDEC = DDEC + DELDEL      ! degrees
c              write(6,*)'dra,ddec',dra,ddec
               XPTS = (DRA*DEG_TO_RAD - X_VAL1) / X_INC1 + X_REF1
               YPTS = (DDEC*DEG_TO_RAD- X_VAL2) / X_INC2 + X_REF2
c	       write(6,*)'ddec,radec_lat',ddec,radec_lat
               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 = NINT(5D1*X) + 1
                  Z = FUNC(ICONO) * RECEIVER_WEIGHT
                  ICO = IORIG + (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
                     RADEC_DATA(ICO+KK) = RADEC_DATA(ICO+KK) +
     $               AZEL_DATA(IOR)*ZZ
                     RADEC_WEIGHT(ICO+KK) = RADEC_WEIGHT(ICO+KK) + ZZ
c                     IF (RADEC_WEIGHT(ICO+KK).GT.0.1) THEN
c                        WRITE(6,*)'kk,zz,x,icono,ico,ico+kk,',
c     $                  KK,ZZ,X,ICONO,ICO,ICO+KK
c                        WRITE(6,*)'RADEC_DATA(ICO+KK),
c     $                  RADEC_WEIGHT(ICO+KK)',
c     $                  RADEC_DATA(ICO+KK),RADEC_WEIGHT(ICO+KK)
c                        STOP
c                     ENDIF
                  ENDDO
               ENDDO
            ENDIF
8           CONTINUE
         ENDDO
      ENDDO
      RETURN
      END
*
      SUBROUTINE MAP_PARAMETER
C*****************************************************************
C
C Initialize parameters for the RADEC map
C
C*****************************************************************
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'convert.inc'
      INCLUDE 'const.inc'
*
      REAL*8   DX,DY
      INTEGER  IX,IY
      DX = RA_EXTENT/120.D0              ! RA_EXTENT and DEC_EXTENT in arcmin.
      DY = DEC_EXTENT/120.D0             ! DX and DY : 1/2 extent (in degrees)
      DX = DX/(XINC_RADEC/3600.D0) + .9999D0
      DY = DY/(YINC_RADEC/3600.D0) + .9999D0
      IX = DX
      IY = DY
      NX_RADEC = 2 * IX + 1
      NY_RADEC = 2 * IY + 1
*
      X_A0 = (SLAM+OLAM)*DEG_TO_RAD
      X_D0 = (SBET+OBET)*DEG_TO_RAD
      X_VAL1 = 0
      X_INC1 = -XINC_RADEC*SEC_TO_RAD
      X_REF1 = IX + 1
*
      X_INC2 = YINC_RADEC*SEC_TO_RAD
      X_VAL2 = 0
      X_REF2 = IY + 1
c	CMAP(80) = 0.D0
      RETURN
      END
*
*
      SUBROUTINE CONVERT_LIMITS(YMAX,LST,EL,OA,OE,RECEIVER_FLAG)
C*****************************************************************
C
C Initialize parameters for the RADEC map
C
C*****************************************************************
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'convert.inc'
      INCLUDE 'const.inc'
      REAL    YMAX,LST,EL
      REAL    OA(1),OE(1)
      LOGICAL RECEIVER_FLAG(1)
*
      REAL    X(4),Y(4)
      REAL*8  P1,P2,P3
      REAL*8  DRA,DDEC,RA50,DEC50,RA_DAY,DEC_DAY
      REAL*8  DECCUR,OBSLAT
      REAL*8  DX,DY
      INTEGER IX,IY
      REAL    RA_MIN,RA_MAX,DEC_MIN,DEC_MAX
      INTEGER I, ICHAN
      REAL    COS_EL,SIN_EL
      INTEGER IREF
*
      CALL SETUP
      DRA = 0.D0
      DDEC = 0.D0
      RA50 = SLAM+OLAM
      DEC50 = SBET+OBET
*
*---  Computes RA_DAY and DEC_DAY
      CALL PREC(3,EPOCHTO,RA_DAY,DEC_DAY,DRA,DDEC,RA50,DEC50)
*
*---  Compute DRA and DDEC (day coordinates)
      DECCUR = DEC_DAY * DEG_TO_RAD
      OBSLAT = SITE_LATITUDE * DEG_TO_RAD
      SINDEC = DSIN(DECCUR)
      COSDEC = DCOS(DECCUR)
      SINOL = DSIN(OBSLAT)
      COSOL = DCOS(OBSLAT)
*
      SIN_EL = SIN(FIELD_ROTATION*EL*DEG_TO_RAD)
      COS_EL = COS(FIELD_ROTATION*EL*DEG_TO_RAD)
      X(1) = XMIN +16                    ! 16 ~ beamsep/2
      Y(1) = YMIN
      X(2) = XMAX -16
      Y(2) = YMIN
      X(3) = XMIN +16
      Y(3) = YMAX
      X(4) = XMAX -16
      Y(4) = YMAX
      RA_MIN  = 1.E32
      RA_MAX  =-1.E32
      DEC_MIN = 1.E32
      DEC_MAX =-1.E32
*
      P3 = LST
      IREF = REF_CHAN
      IF (TELESCOPE(1:3).EQ.'CSO') THEN  ! no rotation
         COS_EL = 1
         SIN_EL = 0
      ENDIF
      DO I=1,4
         DO ICHAN=1,NCHANNEL
            IF (.NOT.RECEIVER_FLAG(ICHAN)) THEN
               P1 = X(I) - (OA(ICHAN)-OA(IREF))*COS_EL
     $         - (OE(ICHAN)-OE(IREF))*SIN_EL
               P2 = Y(I) + (OA(ICHAN)-OA(IREF))*SIN_EL
     $         - (OE(ICHAN)-OE(IREF))*COS_EL
               CALL RADEC(P1,P2,P3,SITE_LATITUDE,RA_DAY,DEC_DAY,
     $         DRA,DDEC)
*---  Computes DRA and DDEC for epoch 1950
               CALL PREC(2,EPOCHTO,RA_DAY,DEC_DAY,DRA,DDEC,RA50,DEC50)
               DRA = DRA*3600.
               DDEC = DDEC*3600.
               IF (DRA.LT.RA_MIN) RA_MIN=DRA
               IF (DRA.GT.RA_MAX) RA_MAX=DRA
               IF (DDEC.LT.DEC_MIN) DEC_MIN=DDEC
               IF (DDEC.GT.DEC_MAX) DEC_MAX=DDEC
            ENDIF
         ENDDO
      ENDDO
*
c      write(6,*) ra_min,ra_max,dec_min,dec_max
      DRA = ABS(RA_MIN)
      IF (ABS(RA_MAX).GT.DRA) DRA = ABS(RA_MAX)! max(ra_min,ra_max)
      DDEC = ABS(DEC_MIN)
      IF (ABS(DEC_MAX).GT.DDEC) DDEC = ABS(DEC_MAX)
*
      DX = DRA/XINC_RADEC + 0.9999D0
      DY = DDEC/YINC_RADEC + 0.9999D0
      IX = DX
      IY = DY
      NX_RADEC = 2*IX +1
      NY_RADEC = 2*IY +1
*
      X_A0 = (SLAM+OLAM)*DEG_TO_RAD
      X_D0 = (SBET+OBET)*DEG_TO_RAD
      X_VAL1 = 0
      X_INC1 = -XINC_RADEC*SEC_TO_RAD
      X_REF1 = IX + 1
*
      X_INC2 = YINC_RADEC*SEC_TO_RAD
      X_VAL2 = 0
      X_REF2 = IY + 1
*
      END
