* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* Subroutines for fit command
*
      SUBROUTINE  COMP_GAUSS(XX,YY,XPAR,DOUBLEBEAM_FIT)
C-----------------------------------------------------------------------
C Compute a model function.
C
C	xx     Real*8    x coordinate in arcsec             Input
C	yy     Real*8    y coordinate in arcsec             Input
C	Xpar   Real*8    Parameters                         Input
C	Y      Real*8    Function (real, imag.)             Output
C	DY     Real*8    Partial derivatives.               Output
C-----------------------------------------------------------------------
C Global:
C ARGUMENTS:
*      INTEGER IFUNC, KFUNC
      REAL*8    XX, YY, XPAR(6)
      REAL*8    REF1, INC1, VAL1, REF2, INC2, VAL2, WIDTH
      REAL*8    Y, DY(8)
      LOGICAL   DOUBLEBEAM_FIT
* Common
      INTEGER*4 IPDATA
      INTEGER*4 NX,NY, MPAR,NDATA,NITE
      LOGICAL   PARAM_OPT, LEVEL_OPT
      REAL      LEVEL
      COMMON/FIT_COMMON/REF1, INC1, VAL1, REF2, INC2, VAL2,
     $WIDTH,Y,DY,IPDATA,NX,NY,MPAR,NDATA,PARAM_OPT,NITE,
     $LEVEL_OPT,LEVEL
      SAVE /FIT_COMMON/
C Local:
      REAL*8  Q1, Q2, A, PI, DPI, DEG, A1, A2, AE, CST, FAR
      REAL*8  Q2_POS,Q2_NEG,A_POS,A_NEG,A1_POS,A1_NEG,A2_POS,A2_NEG
      REAL*8  CA ,SA, XP4T, XP5T
      PARAMETER (PI = 3.14159265358979323846D0, DPI=2D0*PI,
     $CST=2.772588722,                   !4.d0*ln(2.d0)
     $DEG=PI/180D0,
     $FAR=100)                           ! maximum exponent
*
*
*     model = flux* exp (-4*log(2)*(r/b)**2)
*     4 pars (x0, y0, flux, b)
*
      IF (.NOT.DOUBLEBEAM_FIT .AND. MPAR.EQ.2) THEN
         Q2 = (XX-XPAR(1))**2+ (YY-XPAR(2))**2
         A1 = -CST*Q2
         A2 = A1/WIDTH**2
         A = 0
         IF (-A2.LT.100) A = DEXP(A2)
         Y   = A
         DY(1) = -CST*2*(XPAR(1)-XX)/WIDTH**2 * Y
         DY(2) = -CST*2*(XPAR(2)-YY)/WIDTH**2 * Y
*
      ELSEIF (.NOT.DOUBLEBEAM_FIT .AND. MPAR.EQ.4) THEN
         Q2 = (XX-XPAR(1))**2+ (YY-XPAR(2))**2
         A1 = -CST*Q2
         A2 = A1/XPAR(4)**2
         A = 0
         IF (-A2.LT.100) A = DEXP(A2)
         Y   = XPAR(3)*A
         DY(1) = -CST*2*(XPAR(1)-XX)/XPAR(4)**2 * Y
         DY(2) = -CST*2*(XPAR(2)-YY)/XPAR(4)**2 * Y
         DY(3) = A
         DY(4) = -2*A1/XPAR(4)**3 *Y
*
      ELSEIF (.NOT.DOUBLEBEAM_FIT .AND. MPAR.EQ.6) THEN
         XP4T = XPAR(4)**2
         XP5T = XPAR(5)**2
         CA = COS(XPAR(6))
         SA = SIN(XPAR(6))
         Q1 =  CA*(XX-XPAR(1)) + SA*(YY-XPAR(2))
         Q2 = -SA*(XX-XPAR(1)) + CA*(YY-XPAR(2))
         A1 = -CST*Q1**2
         A2 = -CST*Q2**2
         A  = A1*XP4T + A2*XP5T
         AE = 0
         IF (-A.LT.FAR) AE = DEXP(A)
         Y   = XPAR(3)*AE
         DY(1) = 2*CST*(Q1*CA*XP4T-Q2*SA*XP5T) * Y
         DY(2) = 2*CST*(Q1*SA*XP4T+Q2*CA*XP5T) * Y
         DY(3) = AE
         DY(4) = 2*A1*XPAR(4) *Y
         DY(5) = 2*A2*XPAR(5) *Y
         DY(6) = 2*CST*Q1*Q2*(-XP4T + XP5T) * Y
*
      ELSEIF (DOUBLEBEAM_FIT .AND. MPAR.EQ.4) THEN
         Q2_POS  = (XX-XPAR(1))**2+ (YY-XPAR(2))**2
         A1_POS  = -CST*Q2_POS
         A2_POS  = A1_POS/WIDTH**2
         A_POS = 0
         IF (-A2_POS.LT.100) A_POS = DEXP(A2_POS)
         Q2_NEG  = (XX-XPAR(3))**2+ (YY-XPAR(4))**2
         A1_NEG  = -CST*Q2_NEG
         A2_NEG  = A1_NEG/WIDTH**2
         A_NEG = 0
         IF (-A2_NEG.LT.100) A_NEG = -DEXP(A2_NEG)
         Y   = A_POS + A_NEG
         DY(1) = -CST*2*(XPAR(1)-XX)/WIDTH**2 * A_POS
         DY(2) = -CST*2*(XPAR(2)-YY)/WIDTH**2 * A_POS
         DY(3) = -CST*2*(XPAR(3)-XX)/WIDTH**2 * A_NEG
         DY(4) = -CST*2*(XPAR(4)-YY)/WIDTH**2 * A_NEG
*
*
      ELSEIF (DOUBLEBEAM_FIT .AND. MPAR.EQ.5) THEN
         Q2_POS  = (XX-XPAR(1))**2+ (YY-XPAR(2))**2
         A1_POS  = -CST*Q2_POS
         A2_POS  = A1_POS/XPAR(5)**2
         A_POS = 0
         IF (-A2_POS.LT.100) A_POS = DEXP(A2_POS)
         Q2_NEG  = (XX-XPAR(3))**2+ (YY-XPAR(4))**2
         A1_NEG  = -CST*Q2_NEG
         A2_NEG  = A1_NEG/XPAR(5)**2
         A_NEG = 0
         IF (-A2_NEG.LT.100) A_NEG = -DEXP(A2_NEG)
         Y   = A_POS + A_NEG
         DY(1) = -CST*2*(XPAR(1)-XX)/XPAR(5)**2 * A_POS
         DY(2) = -CST*2*(XPAR(2)-YY)/XPAR(5)**2 * A_POS
         DY(3) = -CST*2*(XPAR(3)-XX)/XPAR(5)**2 * A_NEG
         DY(4) = -CST*2*(XPAR(4)-YY)/XPAR(5)**2 * A_NEG
         DY(5) = -2*A1_POS/XPAR(5)**3 * A_POS
     $   -2*A1_NEG/XPAR(5)**3 * A_NEG
*
      ELSEIF (DOUBLEBEAM_FIT .AND. MPAR.EQ.6) THEN
         Q2_POS  = (XX-XPAR(1))**2+ (YY-XPAR(2))**2
         A1_POS  = -CST*Q2_POS
         A2_POS  = A1_POS/WIDTH**2
         A_POS = 0
         IF (-A2_POS.LT.100) A_POS =  DEXP(A2_POS)
         Q2_NEG  = (XX-XPAR(3))**2+ (YY-XPAR(4))**2
         A1_NEG  = -CST*Q2_NEG
         A2_NEG  = A1_NEG/WIDTH**2
         A_NEG = 0
         IF (-A2_NEG.LT.100) A_NEG = DEXP(A2_NEG)
         Y   = XPAR(5)*A_POS + XPAR(6)*A_NEG
         DY(1) = -CST*2*(XPAR(1)-XX)/WIDTH**2 * A_POS *XPAR(5)
         DY(2) = -CST*2*(XPAR(2)-YY)/WIDTH**2 * A_POS* XPAR(5)
         DY(3) = -CST*2*(XPAR(3)-XX)/WIDTH**2 * A_NEG* XPAR(6)
         DY(4) = -CST*2*(XPAR(4)-YY)/WIDTH**2 * A_NEG* XPAR(6)
         DY(5) = A_POS
         DY(6) = A_NEG
      ENDIF
      RETURN
      END
*
      SUBROUTINE DIFF(X,Y,Z)
*-----------------------------------------------------------------------
*
* Compute difference: Z = X-Y . Useful when variables are in allocated area
*
* Input : X,Y
* Output: Z
*-----------------------------------------------------------------------
      REAL*8 X,Z
      REAL*4 Y
      INCLUDE 'parameter.inc'
*
      IF (Y.NE.BLANKING) THEN
         Z = X-Y
      ELSE
         Z = 0
      ENDIF
      RETURN
      END
*
      SUBROUTINE SCALE_DATA(X,Z,N,XMAX,IMAX,IMIN,DOUBLEBEAM_FIT,
     $CNV_MAP)
*-----------------------------------------------------------------------
*
* Find maximum and minimum in an array X(N). Scale X(N) with 1/XMAX
*    Input:  X(N)  = input array
*            N     = size of X array
*
*    Output: Y(N)  = X(N)/XMAX
*            XMAX  = maximum value
*            IMAX  = index of maximum value
*            IMIN  = index of minimum value
*            NDATA = number of unblanked values
*-----------------------------------------------------------------------
*
      INCLUDE 'parameter.inc'
      REAL*4  X(1),Z(1)
      INTEGER N,IMAX,IMIN
      REAL*8  XMAX
      LOGICAL DOUBLEBEAM_FIT,CNV_MAP
* Common data
      INTEGER*4 IPDATA
      REAL*8    REF1, INC1, VAL1, REF2, INC2, VAL2, WIDTH
      INTEGER*4 NX,NY, MPAR,NDATA,NITE
      LOGICAL   PARAM_OPT,LEVEL_OPT
      REAL      LEVEL
      REAL*8    Y,DY(8)
      COMMON/FIT_COMMON/REF1, INC1, VAL1, REF2, INC2, VAL2,
     $WIDTH,Y,DY,IPDATA,NX,NY,MPAR,NDATA,PARAM_OPT,NITE,
     $LEVEL_OPT,LEVEL
      SAVE /FIT_COMMON/
* Local data
      REAL*8 XMIN
      INTEGER*4 I
*
      NDATA = 0
*
* Find maximum and minimum (if param_opt.eq..false)
      IF (.NOT.PARAM_OPT) THEN
         XMIN = 1.D32
         XMAX = -1.D32
         DO I=1,N
            IF (X(I).NE.BLANKING) THEN
               IF (X(I).GT.XMAX) THEN
                  XMAX=X(I)
                  IMAX = I
               ENDIF
               IF (X(I).LT.XMIN) THEN
                  XMIN=X(I)
                  IMIN = I
               ENDIF
            ENDIF
         ENDDO
      ELSE
         XMIN =0
      ENDIF
* Divide by xmax or xmin
      IF (ABS(XMIN).GT.ABS(XMAX) .AND. .NOT.DOUBLEBEAM_FIT .AND.
     $.NOT.CNV_MAP) THEN
         XMAX = XMIN
         IMAX = IMIN
      ENDIF
      IF (ABS(XMAX).GT.1D-32) THEN       ! do not divide by 0 !
         IF (.NOT.LEVEL_OPT) THEN
            DO I=1,N
               IF (X(I).NE.BLANKING) THEN
                  Z(I) = X(I) / XMAX
                  NDATA = NDATA+1
               ELSE
                  Z(I) = BLANKING
               ENDIF
            ENDDO
         ELSE
            DO I=1,N
               IF (X(I).NE.BLANKING .AND. ABS(X(I)/XMAX).GT.LEVEL) THEN
                  Z(I) = X(I) / XMAX
                  NDATA = NDATA+1
               ELSE
                  Z(I) = BLANKING
               ENDIF
            ENDDO
         ENDIF
      ENDIF
      RETURN
      END
*
      LOGICAL FUNCTION TEST_BLANKING(X)
*-----------------------------------------------------------------------
*
* Test blanking value (in allocated area)
*
*-----------------------------------------------------------------------
      INCLUDE 'parameter.inc'
      REAL*4  X
*
      TEST_BLANKING = .FALSE.
      IF (X.EQ.BLANKING) TEST_BLANKING = .TRUE.
      RETURN
      END
*
*
*
      SUBROUTINE REMOVE_FIT_VARIABLE
*-----------------------------------------------------------------------
*
*  Remove previous FIT variable and free associated memory
*
*-----------------------------------------------------------------------
      INCLUDE 'nic.inc'
      LOGICAL ERROR
      CALL SIC_DELVARIABLE('fit_az_pos',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_az_neg',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_el_pos',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_el_neg',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_amp_pos',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_amp_neg',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_fwhm',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_az_off',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_el_off',.FALSE.,ERROR)
*
      CALL SIC_DELVARIABLE('fit_az',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_el',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_RA',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_DEC',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_amp',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_major',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_minor',.FALSE.,ERROR)
      CALL SIC_DELVARIABLE('fit_angle',.FALSE.,ERROR)
*
* Free memory
      IF (IADDR_FIT.NE.0) CALL FREE_VM(LENGTH_FIT,IADDR_FIT)
      IADDR_FIT = 0
      RETURN
      END
*
      SUBROUTINE CREATE_FIT_VARIABLE
*-----------------------------------------------------------------------
*
*  Initialize memory which will be associated to nic variables by
*  command VAR FIT
*
*-----------------------------------------------------------------------
      INCLUDE 'nic.inc'
      INTEGER  SIC_GETVM,POINTER
*
* local variables
      INTEGER IER,I
* common variables
      INTEGER*4 IP_AZ_POS, IP_EL_POS,  IP_AZ_NEG,IP_EL_NEG
      INTEGER*4 IP_WIDTH,  IP_AMP_POS, IP_AMP_NEG
      INTEGER*4 IP_WIDTH1, IP_WIDTH2,  IP_ANGLE
      INTEGER*4 IP_AZ_OFF, IP_EL_OFF
      COMMON /IP_FIT/ IP_AZ_POS, IP_EL_POS,  IP_AZ_NEG,IP_EL_NEG,
     $IP_WIDTH,  IP_AMP_POS, IP_AMP_NEG,
     $IP_WIDTH1, IP_WIDTH2,  IP_ANGLE, IP_AZ_OFF, IP_EL_OFF
      SAVE /IP_FIT/
*
* Allocate memory
      IF (DOUBLEBEAM_FIT) THEN
         LENGTH_FIT = NCHAN_FIT*(7+2)    ! 7 parameters + (az,el) offset compute
      ELSE
         LENGTH_FIT = NCHAN_FIT*(6+2)    ! 6 parameters + (az,el) offset compute
      ENDIF
      IER = SIC_GETVM(LENGTH_FIT, IADDR_FIT)   ! lwork
      IF (IER.NE.1) THEN
         CALL MESSAGE (2,1,'FIT','Memory allocation failure')
         RETURN
      ENDIF
*
* Raz memory
      IP_AZ_POS  = POINTER (IADDR_FIT,MEMORY)
      DO I=0,LENGTH_FIT-1                ! raz memory
         CALL R4TOR4(0,MEMORY(IP_AZ_POS+I),1)
      ENDDO
* Initialize pointers
      IF (DOUBLEBEAM_FIT) THEN
         IP_EL_POS  = IP_AZ_POS  + NCHAN_FIT
         IP_AZ_NEG  = IP_EL_POS  + NCHAN_FIT
         IP_EL_NEG  = IP_AZ_NEG  + NCHAN_FIT
         IP_WIDTH   = IP_EL_NEG  + NCHAN_FIT
         IP_AMP_POS = IP_WIDTH   + NCHAN_FIT
         IP_AMP_NEG = IP_AMP_POS + NCHAN_FIT
         IP_AZ_OFF  = IP_AMP_NEG + NCHAN_FIT
         IP_EL_OFF  = IP_AZ_OFF  + NCHAN_FIT
      ELSE
         IP_EL_POS  = IP_AZ_POS  + NCHAN_FIT
         IP_AMP_POS = IP_EL_POS  + NCHAN_FIT
         IP_WIDTH1  = IP_AMP_POS + NCHAN_FIT
         IP_WIDTH2  = IP_WIDTH1  + NCHAN_FIT
         IP_ANGLE   = IP_WIDTH2  + NCHAN_FIT
         IP_AZ_OFF  = IP_ANGLE   + NCHAN_FIT
         IP_EL_OFF  = IP_AZ_OFF  + NCHAN_FIT
      ENDIF
      RETURN
      END
*
      SUBROUTINE UPDATE_FIT_VARIABLE(XPARAM,ICHAN)
*-----------------------------------------------------------------------
*
*  Update variable
*
*-----------------------------------------------------------------------
      INCLUDE 'nic.inc'
      INCLUDE 'const.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INTEGER ICHAN
      REAL*8  XPARAM(7)
*
* common variables
      INTEGER*4 IP_AZ_POS, IP_EL_POS,  IP_AZ_NEG,IP_EL_NEG
      INTEGER*4 IP_WIDTH,  IP_AMP_POS, IP_AMP_NEG
      INTEGER*4 IP_WIDTH1, IP_WIDTH2,  IP_ANGLE
      INTEGER*4 IP_AZ_OFF, IP_EL_OFF
      COMMON /IP_FIT/ IP_AZ_POS, IP_EL_POS,  IP_AZ_NEG,IP_EL_NEG,
     $IP_WIDTH,  IP_AMP_POS, IP_AMP_NEG,
     $IP_WIDTH1, IP_WIDTH2,  IP_ANGLE, IP_AZ_OFF, IP_EL_OFF
      SAVE /IP_FIT/
* local variables
      REAL    MEAN_ELE, COSELE, SINELE, OFF_AZ, OFF_EL
      REAL*8  VAL1, VAL2
*
      IF (DOUBLEBEAM_FIT) THEN
         CALL R8TOR4(XPARAM(1),MEMORY(IP_AZ_POS +ICHAN-1),1)
         CALL R8TOR4(XPARAM(2),MEMORY(IP_EL_POS +ICHAN-1),1)
         CALL R8TOR4(XPARAM(3),MEMORY(IP_AZ_NEG +ICHAN-1),1)
         CALL R8TOR4(XPARAM(4),MEMORY(IP_EL_NEG +ICHAN-1),1)
         CALL R8TOR4(XPARAM(5),MEMORY(IP_AMP_POS+ICHAN-1),1)
         CALL R8TOR4(XPARAM(6),MEMORY(IP_AMP_NEG+ICHAN-1),1)
         CALL R8TOR4(XPARAM(7),MEMORY(IP_WIDTH  +ICHAN-1),1)
      ELSE
         CALL R8TOR4(XPARAM(1),MEMORY(IP_AZ_POS +ICHAN-1),1)
         CALL R8TOR4(XPARAM(2),MEMORY(IP_EL_POS +ICHAN-1),1)
         CALL R8TOR4(XPARAM(3),MEMORY(IP_AMP_POS+ICHAN-1),1)
         CALL R8TOR4(XPARAM(4),MEMORY(IP_WIDTH1 +ICHAN-1),1)
         CALL R8TOR4(XPARAM(5),MEMORY(IP_WIDTH2 +ICHAN-1),1)
         CALL R8TOR4(XPARAM(6),MEMORY(IP_ANGLE  +ICHAN-1),1)
      ENDIF
*
* Compute channel offsets (if elevation is defined)
      IF (COMPUTE_OFF_FIT) THEN
*
* Old version takes mean elevation => irecord = nrecord/2
*
*         IRECORD = NRECORD/2
*         CALL R4TOR4(MEMORY(IP_EL+IRECORD),MEAN_ELE,1)
*         COSELE = COS(MEAN_ELE*DEG_TO_RAD)
*         SINELE = SIN(MEAN_ELE*DEG_TO_RAD)
*
         IF (DOUBLEBEAM_FIT) THEN
            VAL1 = (XPARAM(1)+XPARAM(3))/2
            VAL2 = (XPARAM(2)+XPARAM(4))/2
            CALL GET_ELEVATION (VAL1, VAL2, MEAN_ELE,
     $      MEMORY(IP_SCAN_COORD),MEMORY(IP_SCAN_COORD+NRECORD),
     $      MEMORY(IP_EL))
            COSELE = COS(MEAN_ELE*DEG_TO_RAD)
            SINELE = SIN(MEAN_ELE*DEG_TO_RAD)
*            off_az =  val1*cosele - val2*sinele
*            off_el = -val1*sinele - val2*cosele
            OFF_AZ = VAL1*COSELE - VAL2*SINELE
            OFF_EL = VAL1*SINELE + VAL2*COSELE
            IF (TELESCOPE(1:3).EQ.'CSO') THEN
               OFF_AZ = VAL1
               OFF_EL = VAL2
            ENDIF
         ELSE
            CALL GET_ELEVATION (XPARAM(1), XPARAM(2), MEAN_ELE,
     $      MEMORY(IP_SCAN_COORD),MEMORY(IP_SCAN_COORD+NRECORD),
     $      MEMORY(IP_EL))
            COSELE = COS(MEAN_ELE*DEG_TO_RAD)
            SINELE = FIELD_ROTATION*SIN(MEAN_ELE*DEG_TO_RAD)
*            off_az =  xparam(1)*cosele - xparam(2)*sinele
*            off_el = -xparam(1)*sinele - xparam(2)*cosele
            OFF_AZ =  XPARAM(1)*COSELE - XPARAM(2)*SINELE
            OFF_EL =  XPARAM(1)*SINELE + XPARAM(2)*COSELE
            IF (TELESCOPE(1:3).EQ.'CSO') THEN
               OFF_AZ = XPARAM(1)
               OFF_EL = XPARAM(2)
            ENDIF
         ENDIF
         CALL R4TOR4(OFF_AZ,MEMORY(IP_AZ_OFF+ICHAN-1),1)
         CALL R4TOR4(OFF_EL,MEMORY(IP_EL_OFF+ICHAN-1),1)
      ENDIF
      RETURN
      END
*
      SUBROUTINE UPDATE_BEAMSEP(FIRST_CHAN,LAST_CHAN)
*-----------------------------------------------------------------------
*
*  Update variable
*
*-----------------------------------------------------------------------
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
* common variables
      INTEGER*4 IP_AZ_POS, IP_EL_POS,  IP_AZ_NEG,IP_EL_NEG
      INTEGER*4 IP_WIDTH,  IP_AMP_POS, IP_AMP_NEG
      INTEGER*4 IP_WIDTH1, IP_WIDTH2,  IP_ANGLE
      INTEGER*4 IP_AZ_OFF, IP_EL_OFF, FIRST_CHAN, LAST_CHAN
      COMMON /IP_FIT/ IP_AZ_POS, IP_EL_POS,  IP_AZ_NEG,IP_EL_NEG,
     $IP_WIDTH,  IP_AMP_POS, IP_AMP_NEG,
     $IP_WIDTH1, IP_WIDTH2,  IP_ANGLE, IP_AZ_OFF, IP_EL_OFF
      SAVE /IP_FIT/
*
* local variables
      INTEGER ICHAN,N
      REAL*4  AZ_POS,AZ_NEG
      REAL*8  BEAMSEP
      LOGICAL ERROR
      CHARACTER*40 STRING
*
      BEAMSEP = 0
      N = 0
      DO ICHAN=FIRST_CHAN,LAST_CHAN
         IF (.NOT.RECEIVER_FLAG(ICHAN)) THEN
            CALL R4TOR4(MEMORY(IP_AZ_POS+ICHAN-1),AZ_POS,1)
            CALL R4TOR4(MEMORY(IP_AZ_NEG+ICHAN-1),AZ_NEG,1)
            IF (AZ_POS.NE.0 .OR. AZ_NEG.NE.0) THEN   ! test available channel
               BEAMSEP = BEAMSEP + ABS(AZ_POS - AZ_NEG)
               N = N+1
            ENDIF
         ENDIF
      ENDDO
      IF (N.NE.0) THEN
         BEAMSEP = BEAMSEP/N
      ELSE
         CALL MESSAGE (2,1,'FIT','No channel available')
      ENDIF
      CALL SIC_LET_DBLE('BEAMSEP',BEAMSEP,ERROR)
      WRITE (STRING,'(a,f8.2,a)') 'BEAMSEP adjusted to',BEAMSEP,
     $' arcseconds'
      CALL MESSAGE (2,1,'FIT',STRING)
      RETURN
      END
*
      SUBROUTINE GET_ELEVATION (AZ_FIT, EL_FIT, MEAN_ELEV,
     $AZ_POS, EL_POS, ELEVATION)
*
*   Get elevation to a fraction of a degree
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      REAL*8  AZ_FIT, EL_FIT
      REAL    MEAN_ELEV, ELEVATION(NRECORD)
      REAL    AZ_POS(NRECORD), EL_POS(NRECORD)
      REAL    BEST, RADIUS, OFFSET, AZ_PEAK, EL_PEAK
      INTEGER I, IPOS_A, IPOS_B
*
      IF (TELESCOPE(1:3).EQ.'CSO') THEN
         MEAN_ELEV = ELEVATION(NRECORD/2)
         RETURN
      ENDIF
*
      AZ_PEAK = AZ_FIT*10
      EL_PEAK = EL_FIT*10
*
      BEST = 1E19
      DO I=1,NRECORD
         RADIUS = (AZ_POS(I)-AZ_PEAK)**2+(EL_POS(I)-EL_PEAK)**2
         IF (RADIUS.LT.BEST) THEN
            IPOS_A = I
            BEST = RADIUS
         ENDIF
      ENDDO
*
* skip one row (= subscan)
*
      OFFSET = ABS(SINT)*10
      IF (EL_POS(IPOS_A).GT.EL_PEAK) OFFSET = -OFFSET
      EL_PEAK = EL_PEAK+OFFSET
      BEST = 1E19
      DO I=1,NRECORD
         RADIUS = (AZ_POS(I)-AZ_PEAK)**2+(EL_POS(I)-EL_PEAK)**2
         IF (RADIUS.LT.BEST) THEN
            IPOS_B = I
            BEST = RADIUS
         ENDIF
      ENDDO
*
      EL_PEAK = EL_PEAK-OFFSET
      MEAN_ELEV = ELEVATION(IPOS_A)*ABS(EL_POS(IPOS_B)-EL_PEAK)+
     $ELEVATION(IPOS_B)*ABS(EL_POS(IPOS_A)-EL_PEAK)
      IF (IPOS_A.EQ.IPOS_B) THEN
         MEAN_ELEV = MEAN_ELEV/2
      ELSE
         MEAN_ELEV = MEAN_ELEV/ABS(EL_POS(IPOS_A)-EL_POS(IPOS_B))
      ENDIF
      RETURN
      END
