* Last processed by NICE on 12-Jun-2000 15:53:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
*
* Written by r.neri July-1994
*
*  Calculates the pointing offsets for multibeam bolometer data
*  Normally channel 1 is used for pointings
*
      SUBROUTINE SOLVE_POINTING (LINE,ERROR)
*
*----------------------------------------------------------------------
*       Support routine for command
*       SOLVE /CHANNEL /ALL /BEAM /TRUNCATE /ITER /RADIUS
* Arguments :
*       LINE    C*(*)   Command line                    Input
*       ERROR   L       Logical error flag              Output
*
*----------------------------------------------------------------------
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      REAL*8  XCOO,YCOO,X8(2),Y8(2)
      LOGICAL ERROR, SIC_PRESENT,REMOVE,STOP_ITE
      CHARACTER*(*) LINE
      INTEGER*4 POINTER, IP_YFIT, IP_WEIAZ, IP_WEIEL
      INTEGER*4 IP_AZI, IP_ELE, IP_SA, IP_SE, IP_DA, IP_DE
      INTEGER*4 YFIT, IP_INP_SI, IP_INP_AZ, IP_INP_EL
      INTEGER*4 IP_OUT_SI, IP_OUT_AZ, IP_OUT_EL, IP
      INTEGER*4 SIGNAL_PLANE, IP_PLANE, IP_SAFIT, IP_SEFIT
      INTEGER*4 IP_DAFIT, IP_DEFIT, IP_AZFIT, IP_ELFIT
      INTEGER*4 IP_SMFIT
      INTEGER*4 IP_REFSIGNAL
      INTEGER  SIC_GETVM, IER, CHAN, I
      INTEGER  DIM(4),NMAX, ITRUNC, ITER, NITER
      PARAMETER (NMAX=2*16+1)
      LOGICAL  ALL_CHAN, GTG_CURS
      REAL     XMIN,XMAX,YMIN,YMAX,A,B,C,PLANE(3),DEFAULT_RADIUS,W,R
      REAL     XINC,YINC,CUTOFF,ELEV_MEAN,RGMIN,RGMAX,X4COO,Y4COO
      REAL     OFFX,OFFY,OFFX_POS,OFFX_NEG,DX,X
      REAL     SUPPX(NMAX),SUPPY(NMAX)
      REAL     BEAMSIZE, AZ_OFFSET
      REAL     PARAM_POINTING(4)         ! [1] : line orientation (degrees)
* [2] : line offset (arcseconds)
* [3] : slope (counts/arcseconds)
* [4] : rms
      REAL     RMS
      PARAMETER (DEFAULT_RADIUS=12)
      INTEGER  NX,NY
      INTEGER  IFUN,IREF
      CHARACTER LIMIT_STRING*40, SETBX_STRING*40
      CHARACTER PLOTS_STRING*40,LABEL_STRING*40
      CHARACTER CHANNEL*2,CODEMOUSE*1
      CHARACTER*60 DRAW_TEXT, DRAW_RELO, DRAW_LINE, TEXT
      CHARACTER GREEK*2, ITALIC*3
      SAVE PLANE
      SAVE PARAM_POINTING
      DATA X8/-100.,100./
*
      ALL_CHAN = SIC_PRESENT(3,0)
      DRAW_TEXT = 'dra tex'
      DRAW_RELO = 'dra rel'
      DRAW_LINE = 'dra lin'
      ITALIC = CHAR(92)//CHAR(92)//'i'   ! '\\'
      GREEK  = CHAR(92)//'g'
*
*       Test if file open
*
      IF (IADDRRECORD.EQ.0) THEN
         CALL MESSAGE(2,1,'SOLVE','No input file opened')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Option /ITER
*
      NITER = 0
      IF (SIC_PRESENT(6,0)) THEN
         CALL SIC_I4 (LINE,6,1,NITER,.TRUE.,ERROR)
         IF (ERROR) THEN
            CALL MESSAGE (4,1,'SOLVE','Invalid number of iteration...')
            RETURN
         ENDIF
         NITER = ABS(NITER)
      ENDIF
      STOP_ITE=.FALSE.
*
* Option /RADIUS
*
      R = DEFAULT_RADIUS
      IF (SIC_PRESENT(7,0)) THEN
         CALL SIC_R4 (LINE,7,1,R,.TRUE.,ERROR)
         IF (ERROR.OR.(R.LT.0)) THEN
            CALL MESSAGE (4,1,'SOLVE','Invalid radius')
         ENDIF
      ENDIF
      IF (ALL_CHAN) GOTO 1000
*
*
* Option /BEAM
*
      BEAMSIZE = 0
      IF (SIC_PRESENT(4,0)) THEN
         CALL SIC_R4 (LINE,4,1,BEAMSIZE,.TRUE.,ERROR)
         IF (ERROR) THEN
            CALL MESSAGE (4,1,'SOLVE','Invalid beamsize...')
            RETURN
         ENDIF
         BEAMSIZE = ABS(BEAMSIZE)
      ENDIF
*
* Option /TRUNCATE
*
      ITRUNC = 0
      IF (SIC_PRESENT(5,0)) THEN
         CALL SIC_I4 (LINE,5,1,ITRUNC,.TRUE.,ERROR)
         IF (ERROR) THEN
            CALL MESSAGE (4,1,'SOLVE','Invalid number of pixels...')
            RETURN
         ENDIF
         ITRUNC = ABS(ITRUNC)
      ENDIF
*
* Option /CHANNEL
*
      CHAN = REFERENCE_RECEIVER
      IF (SIC_PRESENT(1,1).EQV..TRUE.) THEN
         CALL SIC_I4 (LINE,1,1,CHAN,.TRUE.,ERROR)
         IF (CHAN.GT.NCHAN.OR.CHAN.LT.1) THEN
            CALL MESSAGE (2,1,'SOLVE','Bad channel number ')
            ERROR = .TRUE.
            RETURN
         ENDIF
      ENDIF
      WRITE (CHANNEL,'(i2.2)') CHAN
*
      IF (NSUBSCAN/4*4.NE.NSUBSCAN.OR.NSUBSCAN.GT.16) THEN
         CALL MESSAGE (2,1,'SOLVE','Bad number of subscans ')
         RETURN
      ENDIF
      CALL GR_EXEC  ('clear plot')
      CALL GR_EXEC1 ('set box 1.8 14.8 4 15')
      CALL GR_EXEC1 ('set exp .8')
*
      IER = SIC_GETVM (NRECORD*(2*17+3),YFIT)
      IF (IER.NE.1) THEN
         CALL MESSAGE (2,1,'SOLVE','Memory allocation failure')
         RETURN
      ENDIF
      IP_YFIT  = POINTER(YFIT,MEMORY)
      IP_WEIAZ = 2*NRECORD + IP_YFIT
      IP_WEIEL = 2*NRECORD + IP_WEIAZ
      IP_AZI   = 2*NRECORD + IP_WEIEL
      IP_ELE   = 2*NRECORD + IP_AZI
      IP_SA    = 2*NRECORD + IP_ELE
      IP_SE    = 2*NRECORD + IP_SA
      IP_DA    = 2*NRECORD + IP_SE
      IP_DE    = 2*NRECORD + IP_DA
      IP_SAFIT = 2*NRECORD + IP_DE
      IP_SEFIT = 2*NRECORD + IP_SAFIT
      IP_DAFIT = 2*NRECORD + IP_SEFIT
      IP_DEFIT = 2*NRECORD + IP_DAFIT
      IP_AZFIT = 2*NRECORD + IP_DEFIT
      IP_ELFIT = 2*NRECORD + IP_AZFIT
      IP_SMFIT = 2*NRECORD + IP_ELFIT
*
      IP_OUT_SI = IP_SMFIT + 2*NRECORD
      IP_OUT_AZ = IP_OUT_SI + NRECORD
      IP_OUT_EL = IP_OUT_AZ + NRECORD
*
      IP_INP_SI = IP_SIGNAL + (CHAN-1)*NRECORD
      IP_INP_AZ = IP_SCAN_COORD
      IP_INP_EL = IP_INP_AZ + NRECORD
*
      CALL R4TOR4 (MEMORY(IP_INP_AZ),MEMORY(IP_OUT_AZ),NRECORD)
      CALL R4TOR4 (MEMORY(IP_INP_EL),MEMORY(IP_OUT_EL),NRECORD)
      CALL R4TOR4 (MEMORY(IP_INP_SI),MEMORY(IP_OUT_SI),NRECORD)
*
      CALL FIT_POINTING (MEMORY(IP_OUT_AZ),MEMORY(IP_OUT_EL),
     $MEMORY(IP_OUT_SI),CHAN,MEMORY(IP_YFIT),
     $MEMORY(IP_WEIAZ),MEMORY(IP_WEIEL),
     $MEMORY(IP_AZI),MEMORY(IP_ELE),MEMORY(IP_SA),
     $MEMORY(IP_SE),MEMORY(IP_DA),MEMORY(IP_DE),
     $MEMORY(IP_POINTING),MEMORY(IP_SAFIT),MEMORY(IP_SEFIT),
     $MEMORY(IP_DAFIT),MEMORY(IP_DEFIT),MEMORY(IP_AZFIT),
     $MEMORY(IP_ELFIT),MEMORY(IP_SMFIT),NSUBSCAN,BEAMSIZE,
     $ITRUNC)
*
      CALL FREE_VM (NRECORD*(2*17+3),YFIT)
      CALL GR_EXEC1 ('set exp 1.')
*
* /ALL option
*
1000  IF (ALL_CHAN) THEN
         CALL MESSAGE(2,1,'SOLVE',
     $   'Reduction with several channels is experimental')
         IF (IADDR_POINTING.NE.0) THEN
            CALL SIC_DELVARIABLE('pointing_map',.FALSE.,ERROR)
            CALL FREE_VM(LENGTH_POINTING,IADDR_POINTING)
            IADDR_POINTING = 0
         ENDIF
         CALL COMPUTE_OFFSET (MEMORY(IP_SCAN_COORD),AZ_OFFSET)
*
*     Find map dimensions
*
         CALL POINTING_LIMITS(MEMORY(IP_SCAN_COORD),
     $   MEMORY(IP_SCAN_COORD+NRECORD),
     $   MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_EL),1,
     $   XMIN,XMAX,YMIN,YMAX,AZ_OFFSET)
c
         CALL SIC_GET_REAL('MAP_CELL[1]',XINC,ERROR)
         CALL SIC_GET_REAL('MAP_CELL[2]',YINC,ERROR)
         WRITE (TEXT,'(a,2f5.2)') 'Using map_cell= ',XINC, YINC
         CALL MESSAGE (2,1,'SOLVE',TEXT)
         NX = (XMAX-XMIN)/XINC
         NY = (YMAX-YMIN)/YINC
*
*     Allocate memory
*
         LENGTH_POINTING = 2*NX*NY
         IER = SIC_GETVM (2*NX*NY,IADDR_POINTING)
         IF (IER.NE.1) THEN
            CALL MESSAGE (2,1,'SOLVE','Memory allocation failure')
            RETURN
         ENDIF
         IP  = POINTER(IADDR_POINTING,MEMORY)
         DIM(1) = NX
         DIM(2) = NY
         CALL SIC_DEF_REAL('POINTING_MAP',MEMORY(IP),2,DIM,.FALSE.,
     $   ERROR)
         IER = SIC_GETVM (2*NRECORD*NCHAN,SIGNAL_PLANE)
         IF (IER.NE.1) THEN
            CALL MESSAGE (2,1,'SOLVE','Memory allocation failure')
            RETURN
         ENDIF
         IP_PLANE  = POINTER(SIGNAL_PLANE,MEMORY)
         IP_REFSIGNAL = IP_PLANE + NCHAN*NRECORD
         CALL R4TOR4(MEMORY(IP_SIGNAL),MEMORY(IP_REFSIGNAL),
     $   NCHAN*NRECORD)
*
         CALL SIC_GET_REAL('CUTOFF',CUTOFF,ERROR)
         CALL SIC_GET_INTE('TAPER',IFUN,ERROR)
         WRITE (TEXT,'(a,e10.5,a,i2)') 'Using cutoff= ',CUTOFF,
     $   ' and taper= ',IFUN
         CALL MESSAGE (2,1,'SOLVE',TEXT)
         CALL NIC_FUNGEN(FUNC,IFUN)
         IREF = 1
*
* Iteration to find the "best" support
*
         DO ITER=0,NITER
*
*          Define support
            W = (WOBBLER_THROW-1)/2
            IF (R.LT.W) THEN
               DO I=1,NMAX/2
                  SUPPX(I) = COS((I-0.5)*2*PI/(NMAX/2))*R-W
                  SUPPY(I) = SIN((I-0.5)*2*PI/(NMAX/2))*R
                  SUPPX(I+NMAX/2) = W-
     $            COS((NMAX/2-I+0.5)*2*PI/(NMAX/2))*R
                  SUPPY(I+NMAX/2) = SIN((NMAX/2-I+0.5)*2*PI/(NMAX/2))*R
               ENDDO
            ELSE
               DO I=1,NMAX/2
                  SUPPX(I) = COS((I-0.5+NMAX/4)*PI/(NMAX/2))*R-W
                  SUPPY(I) = SIN((I-0.5+NMAX/4)*PI/(NMAX/2))*R
                  SUPPX(I+NMAX/2) = W-
     $            COS((NMAX/2-I+0.5+NMAX/4)*PI/(NMAX/2))*R
                  SUPPY(I+NMAX/2) =
     $            SIN((NMAX/2-I+0.5+NMAX/4)*PI/(NMAX/2))*R
               ENDDO
            ENDIF
            SUPPX(NMAX) = SUPPX(1)
            SUPPY(NMAX) = SUPPY(1)
            OFFX = 0
            OFFY = 0
            IF (ITER.NE.0.AND..NOT.STOP_ITE) THEN
               CALL SIC_GET_REAL ('FIT_AZ_POS[1]',OFFX,ERROR)
               CALL SIC_GET_REAL ('FIT_EL_POS[1]',OFFY,ERROR)
            ENDIF
            DO I=1,NMAX
               SUPPX(I) = SUPPX(I)-AZ_OFFSET/10.+OFFX
               SUPPY(I) = SUPPY(I)+OFFY
            ENDDO
*          Plot map (on the last iteration)
            IF (ITER.EQ.NITER.OR.STOP_ITE) THEN
               CALL POINTING_MAP( MEMORY(IP_SIGNAL),
     $         MEMORY(IP),MEMORY(IP+NX*NY),
     $         MEMORY(IP_SCAN_COORD),
     $         MEMORY(IP_SCAN_COORD+NRECORD),
     $         MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_EL),
     $         MEMORY(IP_FLAG),
     $         IFUN,ERROR,NX,NY,XINC,YINC,XMIN,YMIN,CUTOFF,
     $         IREF,ELEV_MEAN,AZ_OFFSET)
               LIMIT_STRING = 'LIM -000.000 +000.000 -000.000 +000.000'
               SETBX_STRING = 'SET BOX M 000.000 000.000 000.000'
               WRITE (LIMIT_STRING(05:12),'(f8.3)') XMIN-XINC/2
               WRITE (LIMIT_STRING(14:21),'(f8.3)') XMIN+(NX-1)*XINC
               WRITE (LIMIT_STRING(23:30),'(f8.3)') YMIN-YINC/2
               WRITE (LIMIT_STRING(32:39),'(f8.3)') YMIN+(NY-1)*YINC
               WRITE (SETBX_STRING(11:17),'(f7.3)') 4.5
               WRITE (SETBX_STRING(19:25),'(f7.3)') 10.
               WRITE (SETBX_STRING(27:33),'(f7.3)') (XMAX-XMIN)*YINC/
     $         ((YMAX-YMIN)*XINC)
               CALL GR_EXEC ('CLEAR PLOT')
               CALL GR_EXEC ('LIM /RG POINTING_MAP')
               CALL GR_EXEC (SETBX_STRING)
               CALL GR_EXEC ('SET BLA -99999.5 1')
               CALL GR_EXEC ('PLOT')
               CALL GR_EXEC ('SET EXP 1')
               CALL GR_EXEC (LIMIT_STRING)
               CALL GR_EXEC ('BOX')
               LABEL_STRING =
     $         'LABEL "'//ITALIC//GREEK//'D Az (arcsec)" /X'
               CALL GR_EXEC (LABEL_STRING)
               LABEL_STRING =
     $         'LABEL "'//ITALIC//GREEK//'D El (arcsec)" /Y'
               CALL GR_EXEC (LABEL_STRING)
               CALL GR_EXEC ('WEDGE')
               IF (GTG_CURS().AND.NITER.EQ.0) THEN
                  CALL GR_CURS(XCOO,YCOO,X4COO,Y4COO,CODEMOUSE)
                  DO I=1,NMAX
                     SUPPX(I) = SUPPX(I)+XCOO+W+AZ_OFFSET/10.
                     SUPPY(I) = SUPPY(I)+YCOO
                  ENDDO
               ENDIF
               CALL GR_SEGM ('plot',ERROR)
               CALL GR4_CONNECT (NMAX,SUPPX,SUPPY,0.,-1.)
               CALL GR_OUT
            ENDIF
            IF (STOP_ITE) GOTO 2000
*
*      Compute and remove plane
            CALL R4TOR4 (MEMORY(IP_SIGNAL),MEMORY(IP_PLANE),
     $      NRECORD*NCHAN)
            CALL REMOVE_PLANE (MEMORY(IP_PLANE),MEMORY(IP_SIGNAL),
     $      MEMORY(IP_SCAN_COORD), MEMORY(IP_SCAN_COORD+NRECORD),
     $      MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_EL),
     $      MEMORY(IP_FLAG),MEMORY(IP_RMS),IREF,A,B,C,SUPPX,SUPPY,
     $      AZ_OFFSET,RMS)
*
*      Fit plane (if niter.ne.0)
            IF (ITER.NE.NITER) THEN
               CALL POINTING_MAP( MEMORY(IP_SIGNAL),
     $         MEMORY(IP),MEMORY(IP+NX*NY),
     $         MEMORY(IP_SCAN_COORD),
     $         MEMORY(IP_SCAN_COORD+NRECORD),
     $         MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_EL),
     $         MEMORY(IP_FLAG),
     $         IFUN,ERROR,NX,NY,XINC,YINC,XMIN,YMIN,CUTOFF,
     $         IREF,ELEV_MEAN,AZ_OFFSET)
               REMOVE=.FALSE.
               CALL BUILD_MDB_IMAGE(MEMORY(IP),NX,NY,XMIN,YMIN,
     $         XINC,YINC,A,B,C,REMOVE)
               IF (ITER.EQ.0)
     $         CALL MESSAGE(2,1,'SOLVE','Fit on doublebeam image')
               CALL EXEC_ALL ('let type mdb')
               CALL NIC_FIT ('fit',ERROR)
               CALL SIC_GET_REAL ('FIT_AZ_POS[1]',OFFX_POS,ERROR)
               CALL SIC_GET_REAL ('FIT_AZ_NEG[1]',OFFX_NEG,ERROR)
               IF (.NOT.ERROR) THEN      ! fit done
                  DX = ABS(OFFX_POS - OFFX_NEG)
                  IF (ABS(DX-WOBBLER_THROW)/WOBBLER_THROW.GT.0.1)
     $            THEN
                     CALL MESSAGE(2,1,'SOLVE','Incorrect fit ')
                     STOP_ITE = .TRUE.
                  ENDIF
               ELSE                      ! no fit done
                  STOP_ITE = .TRUE.
               ENDIF
               CALL R4TOR4(MEMORY(IP_REFSIGNAL),MEMORY(IP_SIGNAL),
     $         NCHAN*NRECORD)
            ENDIF
         ENDDO
*
*
*  Compute and plot plane map
*
2000     DIM(1) = 3
         DIM(2) = 1
         DIM(3) = 1
         DIM(4) = 1
         CALL SIC_DELVARIABLE('PLANE',.FALSE.,ERROR)
         CALL SIC_DEF_REAL ('PLANE',PLANE,1,DIM,.TRUE.,ERROR)
         PLANE(1) = -A/SQRT(A*A+B*B+1)
         PLANE(2) = -B/SQRT(A*A+B*B+1)
         PLANE(3) = 1/SQRT(A*A+B*B+1)
         WRITE (TEXT,'(A,3F6.3)') 'Found plane= ',PLANE
         CALL MESSAGE(2,1,'SOLVE',TEXT)
         CALL POINTING_MAP(MEMORY(IP_PLANE),
     $   MEMORY(IP),MEMORY(IP+NX*NY),
     $   MEMORY(IP_SCAN_COORD), MEMORY(IP_SCAN_COORD+NRECORD),
     $   MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_EL),
     $   MEMORY(IP_FLAG),
     $   IFUN,ERROR,NX,NY,XINC,YINC,XMIN,YMIN,CUTOFF,IREF,ELEV_MEAN,
     $   AZ_OFFSET)
         CALL FREE_VM (2*NRECORD*NCHAN,SIGNAL_PLANE)
*
         CALL SIC_GET_REAL ('RGMIN',RGMIN,ERROR)
         CALL SIC_GET_REAL ('RGMAX',RGMAX,ERROR)
         PLOTS_STRING = 'PLOT /SCAL LIN -0000000.00 +0000000.00'
         LIMIT_STRING = 'LIM -000.000 +000.000 -000.000 +000.000'
         SETBX_STRING = 'SET BOX M 000.000 000.000 000.000'
         WRITE (PLOTS_STRING(16:26),'(f11.2)') RGMIN
         WRITE (PLOTS_STRING(28:38),'(f11.2)') RGMAX
         WRITE (LIMIT_STRING(05:12),'(f8.3)') XMIN
         WRITE (LIMIT_STRING(14:21),'(f8.3)') XMAX
         WRITE (LIMIT_STRING(23:30),'(f8.3)') YMIN
         WRITE (LIMIT_STRING(32:39),'(f8.3)') YMAX
         WRITE (SETBX_STRING(11:17),'(f7.3)') 17.5
         WRITE (SETBX_STRING(19:25),'(f7.3)') 10.
         WRITE (SETBX_STRING(27:33),'(f7.3)') (XMAX-XMIN)*YINC/
     $   ((YMAX-YMIN)*XINC)
         CALL GR_EXEC ('LIM /RG POINTING_MAP')
         CALL GR_EXEC (SETBX_STRING)
         CALL GR_EXEC (PLOTS_STRING)
         CALL GR_EXEC (LIMIT_STRING)
         CALL GR_EXEC ('BOX P N')
         IF (B.NE.0) THEN
            Y8(1) = -(C+X8(1)*A)/B
            Y8(2) = -(C+X8(2)*A)/B
            CALL GR_SEGM ('plot',ERROR)
            CALL GR8_CONNECT (2,X8,Y8,0D0,-1D0)
            CALL GR_OUT
         ENDIF
         CALL GR_EXEC ('ellipse 5 /user 0 0')
         DIM(1) = 4
         CALL SIC_DELVARIABLE('PARAM_POINTING',.FALSE.,ERROR)
         CALL SIC_DEF_REAL ('PARAM_POINTING',PARAM_POINTING,1,
     $   DIM,.TRUE.,ERROR)
         IF (A.NE.0) THEN
            PARAM_POINTING(1)=ATAN(B/A)/PI*180.! angle between vertical
* and line
            PARAM_POINTING(2) = -C/A
            PARAM_POINTING(3) = SQRT(A**2+B**2)
            PARAM_POINTING(4) = RMS
         ENDIF
         CALL DISPLAY_CHANNEL_POSITION(MEMORY(IP_AZOFF),
     $   MEMORY(IP_ELOFF),
     $   NCHAN,MEMORY(IP_EL+NRECORD/2),REFERENCE_RECEIVER,
     $   FIELD_ROTATION)
*
         LABEL_STRING = 'LABEL "'//ITALIC//GREEK//'D Az (arcsec)" /X'
         CALL GR_EXEC (LABEL_STRING)
         CALL GR_EXEC ('WEDGE')
         CALL GR_EXEC ('SET BOX 4 20 4 10')
         CALL GR_EXEC ('LIM 0 2 0 2')
         DRAW_TEXT = 'DRA TEX .5 1. "\\IElevation = 00 deg" 6 /USE'
         I = INDEX(DRAW_TEXT,'00')
         WRITE (DRAW_TEXT(I:I+1),'(I2)') NINT(ELEV_MEAN)
         CALL GR_EXEC (DRAW_TEXT)
         DRAW_TEXT = 'DRA TEX .5 .8 "\\ISource = xxxxxxxxx" 6 /USE'
         I = INDEX(DRAW_TEXT,'xxxxxxxxx')
         WRITE (DRAW_TEXT(I:I+8),'(A)') SOURCE_NAME(1:9)
         CALL GR_EXEC (DRAW_TEXT)
         DRAW_TEXT = 'DRA TEX .5 .6 "\\IScan = xxxx" 6 /USE'
         I = INDEX(DRAW_TEXT,'xxxx')
         WRITE (DRAW_TEXT(I:I+3),'(I4.4)') SCAN_NUMBER
         CALL GR_EXEC (DRAW_TEXT)
         DRAW_TEXT = 'DRA TE 1.4 1 "\\INormal = (xxxxxx,xxxxxx,xxxxxx)"'
     $   //' 6 /USE'
         I = INDEX(DRAW_TEXT,'=')
         WRITE (DRAW_TEXT(I+3:I+8),'(F6.3)') PLANE(1)
         WRITE (DRAW_TEXT(I+10:I+15),'(F6.3)') PLANE(2)
         WRITE (DRAW_TEXT(I+17:I+22),'(F6.3)') PLANE(3)
         CALL GR_EXEC (DRAW_TEXT)
         DRAW_TEXT = 'DRA T 1.4 .8 "\\ISlope =                        "'
     $   //' 6 /USE'
         I = INDEX(DRAW_TEXT,'=')
         WRITE (DRAW_TEXT(I+1:I+23),'(F7.3,A16)') SQRT(A**2+B**2)
     $   ,' (counts/arcsec)'
         CALL GR_EXEC (DRAW_TEXT)
         DRAW_TEXT = 'DRA TE 1.4 .6 "\\IWobbler throw =         "'
     $   //' 6 /USE'
         WRITE (DRAW_TEXT(35:39),'(F5.1)') WOBBLER_THROW
         CALL GR_EXEC (DRAW_TEXT)
*
* Build mdb image (after the removal of plane)
         IF (NITER.EQ.0) THEN
            CALL POINTING_MAP( MEMORY(IP_SIGNAL),
     $      MEMORY(IP),MEMORY(IP+NX*NY),
     $      MEMORY(IP_SCAN_COORD), MEMORY(IP_SCAN_COORD+NRECORD),
     $      MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_EL),
     $      MEMORY(IP_FLAG),
     $      IFUN,ERROR,NX,NY,XINC,YINC,XMIN,YMIN,CUTOFF,IREF,ELEV_MEAN,
     $      AZ_OFFSET)
            REMOVE=.FALSE.
            CALL BUILD_MDB_IMAGE(MEMORY(IP),NX,NY,XMIN,YMIN,XINC,YINC,
     $      A,B,C,REMOVE)
* Fit
            CALL MESSAGE(2,1,'SOLVE','Fit on doublebeam image')
            CALL EXEC_ALL ('let type mdb')
            CALL NIC_FIT ('fit',ERROR)
            CALL SIC_GET_REAL ('FIT_AZ_POS[1]',OFFX_POS,ERROR)
            CALL SIC_GET_REAL ('FIT_AZ_NEG[1]',OFFX_NEG,ERROR)
            DX = ABS(OFFX_POS - OFFX_NEG)
            IF (ABS(DX-WOBBLER_THROW)/WOBBLER_THROW.GT.0.1) THEN
               CALL MESSAGE(2,1,'SOLVE','Incorrect fit ')
            ENDIF
         ENDIF
         IF (ABS(DX-WOBBLER_THROW)/WOBBLER_THROW.LT.0.1) THEN
            CALL SIC_GET_REAL ('FIT_AZ_POS[1]',X,ERROR)
            CALL R4TOR4(X,MEMORY(IP_POINTING),1)
            CALL SIC_GET_REAL ('FIT_EL_POS[1]',X,ERROR)
            CALL R4TOR4(X,MEMORY(IP_POINTING+2),1)
            CALL SIC_GET_REAL ('FIT_FWHM[1]',X,ERROR)
            CALL R4TOR4(X,MEMORY(IP_POINTING+4),1)
            CALL R4TOR4(X,MEMORY(IP_POINTING+6),1)
            CALL SIC_GET_REAL ('FIT_AMP_POS[1]',X,ERROR)
            CALL R4TOR4(X,MEMORY(IP_POINTING+8),1)
            CALL R4TOR4(X,MEMORY(IP_POINTING+10),1)
         ENDIF
      ENDIF
      ERROR=.FALSE.
      RETURN
*
*99    ERROR = .TRUE.
*      RETURN
      END
*
*----------------------------------------------------------------------------
      SUBROUTINE FIT_POINTING (SCAN_AZ,SCAN_EL,SIGNAL,CHAN,
     $YFIT,WEIAZ,WEIEL,AZ,EL,SA,SE,DA,DE,FIT,SAFIT,SEFIT,
     $DAFIT,DEFIT,AZFIT,ELFIT,SMFIT,NSUB,BEAMSIZE,ITRUNC)
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      REAL*8  YFIT(NRECORD), WEIAZ(NRECORD), WEIEL(NRECORD)
      REAL*8  AZ(NRECORD), EL(NRECORD), SMFIT(NRECORD)
      REAL*8  SAFIT(NRECORD),SEFIT(NRECORD),AZFIT(NRECORD)
      REAL*8  DAFIT(NRECORD),DEFIT(NRECORD),ELFIT(NRECORD)
      REAL*8  SA(NRECORD), SE(NRECORD), DA(NRECORD), DE(NRECORD)
      REAL*4  SCAN_AZ(NRECORD), SCAN_EL(NRECORD), SIGNAL(NRECORD)
      REAL*4  FIT(12),BEAMSIZE
      INTEGER CHAN, IAN, IEN, NSUB
*
      REAL*8  OFF_AZ, OFF_EL, WIDTH_AZ, WIDTH_EL, INT_AZ, INT_EL
      REAL*8  CHISQR, OLD_CHISQR, DELTA_AA(5), SIGMA_AA(5), A(5)
      REAL*4  PAR1, PAR3, MINX, MINY, MAXX, MAXY
      REAL*4  PAR2, PAR4, AMPL, MEANSA, MEANDA, MEANSE, MEANDE
      REAL*4  ESTIMATED_COUNTS_AZ, ESTIMATED_COUNTS_EL, VAL1
      INTEGER IA1(6), IA2(6), IE1(6), IE2(6), IA, IE, ID, IK
      INTEGER NTERMS, MODE, LENC, IC, I, J, L, M, IL, ITRUNC
      CHARACTER*120 DRAW_TEXT
      LOGICAL ERROR
      REAL*4  SA_R4,SE_R4,AZ_R4,EL_R4
      CHARACTER GREEK*2, ITALIC*3
*
      ITALIC = CHAR(92)//CHAR(92)//'i'   ! '\\'
      GREEK  = CHAR(92)//'g'
      NTERMS = 5                         ! number of fit parameters
      MODE   = 1                         ! instrumental weighting
      DRAW_TEXT = ' '
*
      DO I=1,NRECORD
         SMFIT(I) = 0
         SAFIT(I) = 0
      ENDDO
      ID = NINT(ABS(SCAN_AZ(1)-SCAN_AZ(NRECORD))/10D0)+5
      VAL1 = SCAN_EL(1)
      IA = 0
      IA2(1) = 1
      DO L=1,NINT(NSUB/4.0)
         IC = 0
         DO WHILE (IC.LE.3)
            IA1(L) = IA2(L)
            IC = 0
            DO WHILE (VAL1.EQ.SCAN_EL(IA2(L)).AND.L.LT.NRECORD)
               IA2(L) = IA2(L)+1
               IC = IC+1
            ENDDO
            IA2(L) = IA2(L)+1
         ENDDO
         IA2(L) = IA2(L)-2
         MEANSA = 0
         MEANSE = 0
         MEANDA = 0
         MEANDE = 0
         MINX = -BLANKING                ! blanking = -99999.5 ?
         MAXX =  BLANKING
         MINY = -BLANKING
         MAXY =  BLANKING
         PAR1 =  BLANKING
         PAR3 =  BLANKING
         DO I=IA1(L),IA2(L)
            IA = IA+1
            AZ(IA) = (SCAN_AZ(I)-SCAN_AZ(NRECORD))/10D0
            SA(IA) = SIGNAL(I)
            IF (ABS(SA(IA)-BLANKING).GT.0.5) THEN
               M = NINT(AZ(IA))+ID+1
               SMFIT(M) = SMFIT(M)+1
               SAFIT(M) = SAFIT(M)+SIGNAL(I)   ! mean
            ENDIF
         ENDDO
         IA2(L+1) = IA2(L)+1
      ENDDO
      IF (ITRUNC.NE.0) THEN
         IK = 1
         DO WHILE (SMFIT(IK).EQ.0)
            IK = IK+1
         ENDDO
         DO L=1,IK+ITRUNC-1
            SMFIT(L) = 0
         ENDDO
         MINX = IK-(ID+1)
         IK = NRECORD
         DO WHILE (SMFIT(IK).EQ.0)
            IK = IK-1
         ENDDO
         DO L=IK-ITRUNC+1,NRECORD
            SMFIT(L) = 0
         ENDDO
         MAXX = IK-(ID+1)
      ENDIF
      IK = 0
      DO L=1,NRECORD
         IF (SMFIT(L).NE.0) THEN
            IK = IK+1
            AZFIT(IK) = L-(ID+1)
            WEIAZ(IK) = SQRT(SMFIT(L))
            SAFIT(IK) = SAFIT(L)/SMFIT(L)
            DAFIT(IK) = AZFIT(IK)*SAFIT(IK)
            MEANSA = MEANSA+SAFIT(IK)
            MEANDA = MEANDA+DAFIT(IK)
            SA_R4 = SAFIT(IK)            ! to be used as argument of amax1 funct
            AZ_R4 = AZFIT(IK)            ! to be used as argument of amax1 funct
            PAR1    = AMAX1 (PAR1,SA_R4)
            MINX = AMIN1 (MINX,AZ_R4)
            MAXX = AMAX1 (MAXX,AZ_R4)
            MINY = AMIN1 (MINY,SA_R4)
            MAXY = AMAX1 (MAXY,SA_R4)
         ENDIF
      ENDDO
*
      DO L=1,NINT(NSUB/4.0)
         IE1(L) = IA2(L)+1
         IF (L+1.LE.NINT(NSUB/4.0)) IE2(L) = IA1(L+1)-1
      ENDDO
      IE2(NINT(NSUB/4.0)) = NRECORD
*
      DO I=1,NRECORD
         SMFIT(I) = 0
         SEFIT(I) = 0
      ENDDO
      ID = NINT(ABS(SCAN_EL(NRECORD))/10D0)+5
      IE = 0
      DO L=1,NINT(NSUB/4.0)
         DO I=IE1(L),IE2(L)
            IE = IE+1
            EL(IE) = SCAN_EL(I)/10D0     ! deciarcsec
            SE(IE) = SIGNAL(I)
            IF (ABS(SE(IE)-BLANKING).GT.0.5) THEN
               M = NINT(EL(IE))+ID+1
               SMFIT(M) = SMFIT(M)+1
               SEFIT(M) = SEFIT(M)+SIGNAL(I)   ! mean
            ENDIF
         ENDDO
      ENDDO
      IF (ITRUNC.NE.0) THEN
         IL = 1
         DO WHILE (SMFIT(IL).EQ.0)
            IL = IL+1
         ENDDO
         DO L=1,IL+ITRUNC-1
            SMFIT(L) = 0
         ENDDO
         MINX = AMIN1(MINX,REAL(IL-(ID+1)))
         IL = NRECORD
         DO WHILE (SMFIT(IL).EQ.0)
            IL = IL-1
         ENDDO
         DO L=IL-ITRUNC+1,NRECORD
            SMFIT(L) = 0
         ENDDO
         MAXX = AMAX1(MAXX,REAL(IL-(ID+1)))
      ENDIF
      IL = 0
      DO L=1,NRECORD
         IF (SMFIT(L).NE.0) THEN
            IL = IL+1
            ELFIT(IL) = L-(ID+1)
            WEIEL(IL) = SQRT(SMFIT(L))
            SEFIT(IL) = SEFIT(L)/SMFIT(L)
            DEFIT(IL) = ELFIT(IL)*SEFIT(IL)
            MEANSE = MEANSE+SEFIT(IL)
            MEANDE = MEANDE+DEFIT(IL)
            SE_R4 = SEFIT(IL)            ! to be used as argument of amax1 funct
            EL_R4 = ELFIT(IL)            ! to be used as argument of amax1 funct
            PAR1    = AMAX1 (PAR1,SE_R4)
            MINX = AMIN1 (MINX,EL_R4)
            MAXX = AMAX1 (MAXX,EL_R4)
            MINY = AMIN1 (MINY,SE_R4)
            MAXY = AMAX1 (MAXY,SE_R4)
         ENDIF
      ENDDO
      PAR2 = MEANDA/MEANSA
      PAR4 = MEANDE/MEANSE
      AMPL = MAXX-MINX
      MINX = MINX-.05*AMPL
      MAXX = MAXX+.05*AMPL
      AMPL = MAXY-MINY
      MINY = MINY-.05*AMPL
      MAXY = MAXY+.05*AMPL
      CALL GR_LIMI (4,MINX,MAXX,MINY,MAXY)
      CALL GR_EXEC1 ('box')
      CALL GR_EXEC1 ('label "Azimuth (arcsec)" /x')
      CALL GR_EXEC1 ('set orie 90')
      IF (JANSKY_FLAG) THEN
         CALL GR_EXEC1 ('dra tex -1.62 0 "Flux (Jy)" 5 /box 4')
      ELSE
         CALL GR_EXEC1 ('dra tex -1.62 0 "Counts" 5 /box 4')
      ENDIF
      CALL GR_EXEC1 ('set orie 0')
      A(2) = PAR2
      IF (BEAMSIZE.NE.0) THEN
         A(3) = 2*BEAMSIZE**2/(8*LOG(2.))
      ELSE
         A(3) = 48                       ! 2*beam^2/(8*log(2)) arcsec
      ENDIF
      A(4) = (SAFIT(1)+SAFIT(2)+SAFIT(IK-1)+SAFIT(IK))/4.
      A(5) = (SAFIT(IK-1)+SAFIT(IK)-(SAFIT(1)+SAFIT(2)))/
     $(AZFIT(IK)-AZFIT(1))/2.
      A(1) = MAXY-A(4)
      DELTA_AA(1) = 10
      DELTA_AA(2) =  1
      IF (BEAMSIZE.NE.0) THEN
         DELTA_AA(3) = 0
      ELSE
         DELTA_AA(3) =  .5
      ENDIF
      DELTA_AA(4) = 10
      DELTA_AA(5) =  1
      CHISQR = BLANKING
100   OLD_CHISQR = CHISQR
      CALL GRIDLS (AZFIT,SAFIT,WEIAZ,IK,NTERMS,MODE,A,DELTA_AA,
     $SIGMA_AA,YFIT,CHISQR)
      IF (ABS(OLD_CHISQR-CHISQR)/CHISQR.GT.1E-3) GOTO 100
      CALL GR_EXEC ('PEN 0')
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_MARKER (IA,AZ,SA,D_BLANKING,10D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 0 /WEI 4')
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_CONNECT (IK,AZFIT,YFIT,0D0,-1D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 0 /DEF')
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_HISTO (IK,AZFIT,SAFIT,0D0,-1D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 1')
      IAN = NINT(IA/2.)
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_MARKER (IAN,AZ,SA,D_BLANKING,10D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 0')
      INT_AZ = YFIT(1)
      J = 1
      DO I=2,IK
         IF (INT_AZ.LT.YFIT(I)) THEN
            INT_AZ = YFIT(I)
            J = I
         ENDIF
      ENDDO
      INT_AZ = INT_AZ-(A(4)-A(5)*AZFIT(J))
      OFF_AZ = A(2)
      WIDTH_AZ = ABS(A(3))
      BEAM_AZIMUTH = 2*SQRT(WIDTH_AZ*LOG(2.))
      FIT(2) =  0
      IF (BEAMSIZE.EQ.0) THEN
         FIT(6) =  0
      ELSE
         FIT(6) =  0
      ENDIF
      FIT(10) = 0
*
      CALL GR_EXEC1 ('set box 15.2 28.2 4 15')
      CALL GR_EXEC1 ('box p n')
      CALL GR_EXEC1 ('label "Elevation (arcsec)" /x')
      A(2) = PAR4
      IF (BEAMSIZE.NE.0) THEN
         A(3) = 2*BEAMSIZE**2/(8*LOG(2.))
      ELSE
         A(3) = 48                       ! 2*beam^2/(8*log(2)) arcsec
      ENDIF
      A(4) = (SEFIT(1)+SEFIT(2)+SEFIT(IK-1)+SEFIT(IK))/4.
      A(5) = (SEFIT(IK-1)+SEFIT(IK)-(SEFIT(1)+SEFIT(2)))/
     $(ELFIT(IK)-ELFIT(1))/2.
      A(1) = MAXY-A(4)
      DELTA_AA(1) = 10
      DELTA_AA(2) =  1
      IF (BEAMSIZE.NE.0) THEN
         DELTA_AA(3) = .0
      ELSE
         DELTA_AA(3) =  .5
      ENDIF
      DELTA_AA(4) = 10
      DELTA_AA(5) =  1
      CHISQR = BLANKING
200   OLD_CHISQR = CHISQR
      CALL GRIDLS (ELFIT,SEFIT,WEIEL,IL,NTERMS,MODE,A,DELTA_AA,SIGMA_AA,
     $YFIT,CHISQR)
      IF (ABS(OLD_CHISQR-CHISQR)/CHISQR.GT.1E-3) GOTO 200
      FIT(4) =  0
      IF (BEAMSIZE.EQ.0) THEN
         FIT(8) =  0
      ELSE
         FIT(8) =  0
      ENDIF
      FIT(12) = 0
      CALL GR_EXEC ('PEN 0')
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_MARKER (IE,EL,SE,D_BLANKING,10D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 0 /WEI 4')
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_CONNECT (IL,ELFIT,YFIT,0D0,-1D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 0 /DEF')
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_HISTO (IL,ELFIT,SEFIT,0D0,-1D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 1')
      IEN = NINT(IE/2.)
      CALL GR_SEGM ('plot',ERROR)
      CALL GR8_MARKER (IEN,EL,SE,D_BLANKING,10D0)
      CALL GR_OUT
      CALL GR_EXEC ('PEN 0')
      INT_EL = YFIT(1)
      J = 1
      DO I=2,IL
         IF (INT_EL.LT.YFIT(I)) THEN
            INT_EL = YFIT(I)
            J = I
         ENDIF
      ENDDO
      INT_EL = INT_EL-(A(4)-A(5)*ELFIT(J))
      OFF_EL = A(2)
      WIDTH_EL = ABS(A(3))
      BEAM_ELEVATION = 2*SQRT(WIDTH_EL*LOG(2.))
*
*      IF (WIDTH_EL.GT.3.AND.ABS(OFF_EL).LT.20)
*     $ESTIMATED_COUNTS_AZ = INT_AZ/EXP(-OFF_EL**2/WIDTH_EL)
      ESTIMATED_COUNTS_AZ = INT_AZ
*      IF (WIDTH_AZ.GT.3.AND.ABS(OFF_AZ).LT.20)
*     $ESTIMATED_COUNTS_EL = INT_EL/EXP(-OFF_AZ**2/WIDTH_AZ)
      ESTIMATED_COUNTS_EL = INT_EL
*
      CALL GR_EXEC1 ('set box 1.8 28.2 4 15')
      CALL GR_EXEC1 ('lim -1 1 0 1')
      CALL GR_EXEC1 ('set exp 1.')
      DRAW_TEXT(1:25) = 'dra tex 0 1.2 "'//ITALIC//'SCAN : '
      DRAW_TEXT(30:44) = '  -  CHANNEL : '
      DRAW_TEXT(47:) = '  -  SOURCE : '//
     $SOURCE_NAME(1:LENC(SOURCE_NAME))//'" 5 /user'
      WRITE (DRAW_TEXT(26:29),'(i4)') SCAN_NUMBER
      WRITE (DRAW_TEXT(45:46),'(i2)') CHAN
      CALL GR_EXEC1 (DRAW_TEXT)
      DRAW_TEXT = '     '
      CALL GR_EXEC1 ('set exp .9')
      DRAW_TEXT(1:) = 'dra tex -1e+0 1.10 "'//ITALIC//GREEK//
     $'DAZ = -00.00" 6 /user'
      WRITE (DRAW_TEXT(32:37),'(f6.2)') OFF_AZ
      FIT(1) = OFF_AZ
      CALL GR_EXEC1 (DRAW_TEXT)
      DRAW_TEXT(1:) = 'dra tex -5e-1 1.10 "'//ITALIC//
     $'INTENSITY = 0000000.0" 6 /user'
      IF (ESTIMATED_COUNTS_AZ.LT.100) THEN
         WRITE (DRAW_TEXT(36:44),'(f9.2)') ESTIMATED_COUNTS_AZ
      ELSE
         WRITE (DRAW_TEXT(36:44),'(f9.1)') ESTIMATED_COUNTS_AZ
      ENDIF
      FIT(9) = ESTIMATED_COUNTS_AZ
      CALL GR_EXEC1 (DRAW_TEXT)
      DRAW_TEXT(1:) = 'dra tex -1e+0 1.05 "'//ITALIC//
     $'FWHM = -00.00" 6 /user'
      WRITE (DRAW_TEXT(31:36),'(f6.2)') BEAM_AZIMUTH
      FIT(5) = BEAM_AZIMUTH
      CALL GR_EXEC1 (DRAW_TEXT)
      DRAW_TEXT(1:) = 'dra tex  1e-2 1.10 "'//ITALIC//GREEK//
     $'DEL = -00.00" 6 /user'
      WRITE (DRAW_TEXT(32:37),'(f6.2)') OFF_EL
      FIT(3) = OFF_EL
      CALL GR_EXEC1 (DRAW_TEXT)
      DRAW_TEXT(1:) = 'dra tex +5e-1 1.10 "'//ITALIC//
     $'INTENSITY = 0000000.0" 6 /user'
      IF (ESTIMATED_COUNTS_EL.LT.100) THEN
         WRITE (DRAW_TEXT(36:44),'(f9.2)') ESTIMATED_COUNTS_EL
      ELSE
         WRITE (DRAW_TEXT(36:44),'(f9.1)') ESTIMATED_COUNTS_EL
      ENDIF
      FIT(11) = ESTIMATED_COUNTS_EL
      CALL GR_EXEC1 (DRAW_TEXT)
      DRAW_TEXT(1:) = 'dra tex  1e-2 1.05 "'//ITALIC//
     $'FWHM = -00.00" 6 /user'
      WRITE (DRAW_TEXT(31:36),'(f6.2)') BEAM_ELEVATION
      FIT(7) = BEAM_ELEVATION
      CALL GR_EXEC1 (DRAW_TEXT)
      IF (DATE_OF_OBSERVATION.GE.19981207) THEN
         DRAW_TEXT(1:) = 'dra tex  -5e-1 1.05 "'//ITALIC//
     $   'COL* = -000.00 (FITTED)" 6 /user'
         WRITE (DRAW_TEXT(32:38),'(f7.2)') COLSTAR+OFF_AZ
         CALL GR_EXEC1 (DRAW_TEXT)
         DRAW_TEXT(1:) = 'dra tex   5e-1 1.05 "'//ITALIC//
     $   'NULE = -000.00 (FITTED)" 6 /user'
         WRITE (DRAW_TEXT(32:38),'(f7.2)') PNULE+OFF_EL
         CALL GR_EXEC1 (DRAW_TEXT)
      ENDIF
      CALL GR_LIMI (4,MINX,MAXX,MINY,MAXY)
      RETURN
      END
*----------------------------------------------------------------------------
* Adapted from Bevington
*
      SUBROUTINE GRIDLS (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA,SIGMAA,
     $YFIT,CHISQR)
      INTEGER NPTS,MODE,NTERMS,I,J,NFREE,LOOP_END
      REAL*8 X(NPTS),Y(NPTS),SIGMAY(NPTS),CHISQ3,VAL
      REAL*8 A(NTERMS),DELTAA(NTERMS),SIGMAA(NTERMS),YFIT(NPTS)
      REAL*8 FCHISQ,FUNCTN,FN,DELTA,SAVE,CHISQR,CHISQ1,CHISQ2
      LOGICAL EXIT_LOOP
*
      NFREE = NPTS-NTERMS
      CHISQR = 0
*
      IF (NFREE.GT.0) THEN
         DO J = 1,NTERMS
            DO I = 1,NPTS
               YFIT(I) = FUNCTN (X,I,A)
            ENDDO
            CHISQ1 = FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
            CHISQ2 = CHISQ1
            FN = 0
            DELTA = DELTAA(J)
            LOOP_END = 1
            DO WHILE (CHISQ1.EQ.CHISQ2.AND.LOOP_END.NE.200)
               A(J) = A(J)+DELTA
               LOOP_END = LOOP_END+1
               DO I=1,NPTS
                  YFIT(I) = FUNCTN (X,I,A)
               ENDDO
               CHISQ2 = FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
            ENDDO
            IF (CHISQ1.LT.CHISQ2) THEN
               DELTA = -DELTA
               A(J) = A(J)+DELTA
               DO I=1,NPTS
                  YFIT(I) = FUNCTN (X,I,A)
               ENDDO
               SAVE = CHISQ1
               CHISQ1 = CHISQ2
               CHISQ2 = SAVE
            ENDIF
            EXIT_LOOP = .TRUE.
            LOOP_END = 1
            DO WHILE (EXIT_LOOP)
               FN = FN+1
               A(J) = A(J)+DELTA
               DO I=1,NPTS
                  YFIT(I) = FUNCTN (X,I,A)
               ENDDO
               CHISQ3 = FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
               EXIT_LOOP = .FALSE.
               IF (CHISQ3.LT.CHISQ2) THEN
                  CHISQ1 = CHISQ2
                  CHISQ2 = CHISQ3
                  LOOP_END = LOOP_END+1
                  IF (LOOP_END.NE.200) EXIT_LOOP = .TRUE.
               ENDIF
            ENDDO
            IF (CHISQ3.NE.CHISQ2) THEN
               DELTA = DELTA *(1./(1+(CHISQ1-CHISQ2)/
     $         (CHISQ3-CHISQ2))+0.5)
               A(J) = A(J)-DELTA
            ELSE
               DELTA = 0
            ENDIF
            VAL = CHISQ3-2*CHISQ2+CHISQ1
            IF (VAL.GT.0) THEN
               SIGMAA(J) = DELTAA(J)*SQRT(2./(NFREE*VAL))
            ELSE
               SIGMAA(J) = 1E+7          ! fix a reasonable maximum
            ENDIF
            DELTAA(J) = DELTAA(J)*FN/3
         ENDDO
         DO I=1,NPTS
            YFIT(I) = FUNCTN(X,I,A)
         ENDDO
         CHISQR = FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
      ENDIF
      RETURN
      END
*
      REAL*8 FUNCTION FCHISQ(Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
c ---------------------------------------------------------------------------
c
c	purpose
c		evaluate reduced chi square for fit to data
c		fchisq = sum ((y-yfit)**2 / sigma**2) / nfree
c
c	usage
c		result = fchisq (y, sigmay, npts, nfree, mode, yfit)
c
c	description of parameters
c		y	- array of data points
c		sigmay	- array of standard deviations for data points
c		npts	- number of data points
c		mode	- determines method of weighting least-squares fit
c			  +1 (instrumental) weight(i) = 1./sigma(i)**2
c			   0 (no weighting) weight(i) = 1.
c			  -1 (statistical)  weight(i) = 1./y(i)
c		yfit	- array of calculated values of fitted function
c
c	subroutines and functions required
c		none
c
c -----------------------------------------------------------------------------
      INTEGER NPTS,MODE,I,NFREE
      REAL*8  Y(NPTS),SIGMAY(NPTS),YFIT(NPTS)
      REAL*8 CHISQ,WEIGHT
*
c	if (nfree.le.0) stop 'fchisq: nfree <= 0'
*
c -----	accumulate chi square
*
      CHISQ = 0
      DO I = 1,NPTS
         WEIGHT = 1
         IF (MODE.LT.0.AND.Y(I).NE.0) WEIGHT = ABS(1/Y(I))
         IF (MODE.GT.0) WEIGHT = 1/SIGMAY(I)**2
         CHISQ = CHISQ+WEIGHT*(Y(I)-YFIT(I))**2
      ENDDO
*
c -----	divide by number of degrees of freedom.
*
      FCHISQ = CHISQ/NFREE
      RETURN
      END
*
      REAL*8 FUNCTION FUNCTN (X,I,A)
      REAL*8 X(1),A(5),E,XI
      INTEGER I
*
* ----- gaussian + linear fitting function
*
      XI = X(I)
      E = (A(2)-XI)**2/(ABS(A(3))+1E-9)  ! to prevent floating exceptions
      IF (E.LT.500) FUNCTN = A(1)*EXP(-E)+A(4)+A(5)*XI
      RETURN
      END
*
*
*----------------------------------------------------------------------------
*
      SUBROUTINE POINTING_LIMITS(XCOORD,YCOORD,OA,OE,EL,IREF,
     $XMIN,XMAX,YMIN,YMAX,AZ_OFFSET)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      REAL    XCOORD(1),YCOORD(1)
      REAL    OA(1),OE(1),EL(1), AZ_OFFSET
      INTEGER IREF                       ! reference channel
      REAL    XMIN,XMAX,YMIN,YMAX
*
      INTEGER IREC,ICHAN
      REAL    COS_EL,SIN_EL
      REAL    XCOO,YCOO
*
      XMIN = 1.D32
      XMAX = -1.D32
      YMIN = 1.D32
      YMAX = -1.D32
      DO IREC=1,NRECORD
         COS_EL = COS(EL(IREC)*DEG_TO_RAD)
         SIN_EL = SIN(EL(IREC)*DEG_TO_RAD)
         DO ICHAN = 1,NCHAN
            XCOO = (XCOORD(IREC)-AZ_OFFSET)/10.
     $      +(OA(IREF)-OA(ICHAN))*COS_EL
     $      +(OE(IREF)-OE(ICHAN))*SIN_EL
            YCOO = YCOORD(IREC)/10.
     $      -(OA(IREF)-OA(ICHAN))*SIN_EL
     $      +(OE(IREF)-OE(ICHAN))*COS_EL
            IF (XMIN.GT.XCOO) XMIN = XCOO
            IF (XMAX.LT.XCOO) XMAX = XCOO
            IF (YMIN.GT.YCOO) YMIN = YCOO
            IF (YMAX.LT.YCOO) YMAX = YCOO
         ENDDO
      ENDDO
!      XMIN = XMIN
!      XMAX = XMAX
!      YMIN = YMIN
!      YMAX = YMAX
!      WRITE(6,*)'xmin,xmax,ymin,ymax',XMIN,XMAX,YMIN,YMAX
      END
*----------------------------------------------------------------------------
*
      SUBROUTINE  POINTING_MAP(SIGNAL,MAP_DATA,MAP_WEIGHT,
     $XCOORD,YCOORD,OA,OE,EL,CHAN_FLAG,IFUN,ERROR,NX,NY,
     $XINC,YINC,XMIN,YMIN,CUTOFF,IREF,ELEV_MEAN,AZ_OFFSET)
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'const.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      REAL    SIGNAL(NRECORD,NCHAN)
      REAL    MAP_DATA(1)                ! output
      REAL    MAP_WEIGHT(1)              ! output
      REAL    XCOORD(1),YCOORD(1),OA(1),OE(1),EL(1)
      LOGICAL ERROR
      INTEGER NX,NY,IFUN,IREF
      REAL    XINC,YINC,XMIN,YMIN
      REAL    CUTOFF, ELEV_MEAN, AZ_OFFSET
      LOGICAL CHAN_FLAG(1)
*
* local data
      REAL    MAGI
      REAL    XINV,YINV,FUNWX,FUNWY
      REAL    XINVFX,XINVFY
      REAL    TOTWX,TOTWY,XLIMM,YLIMM
      REAL    ZZ,Z
      REAL    VAL
      INTEGER IREC,ICHAN,I
      INTEGER IFIRST_X,ILAST_X,IFIRST_Y,ILAST_Y
      INTEGER JJ,ICONO,ICO,KK
      REAL    XPTS,YPTS
      REAL    SOFFY,SOFFX,X
      REAL*8  W_MIN,W_MAX
      REAL*8  ALIM,ALOW
      REAL    COS_EL,SIN_EL
* Initializations
      XINV=1./XINC
      YINV=1./YINC
*
      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
*
      FUNWX = VLIGHT/(FREQUENCY*1.E9)/(MAGI*TELESCOPE_DIAMETER)
     $/(XINC*SEC_TO_RAD)                 ! lambda/(1.67*diam_telescope)/xinc
      FUNWY = VLIGHT/(FREQUENCY*1.E9)/(MAGI*TELESCOPE_DIAMETER)
     $/(YINC*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-1.D0+TOTWX
      YLIMM = NY-1.D0+TOTWY
*
*
      DO I=1, NX*NY
         MAP_DATA(I) = 0
         MAP_WEIGHT(I) = 0
      ENDDO
*
* Add scan data to map_data and map_weight
*
      DO ICHAN =1,NCHAN
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            ELEV_MEAN = 0
            DO IREC = 1, NRECORD
               ELEV_MEAN = ELEV_MEAN+EL(IREC)
               VAL = SIGNAL(IREC,ICHAN)
               IF (ABS(VAL-BLANKING).GT.0.5) THEN
                  COS_EL = COS(FIELD_ROTATION*EL(IREC)*DEG_TO_RAD)
                  SIN_EL = SIN(FIELD_ROTATION*EL(IREC)*DEG_TO_RAD)
                  XPTS = (XCOORD(IREC)-AZ_OFFSET)/10
     $            +(OA(IREF)-OA(ICHAN))*COS_EL
     $            +(OE(IREF)-OE(ICHAN))*SIN_EL
                  YPTS = YCOORD(IREC)/10
     $            -(OA(IREF)-OA(ICHAN))*SIN_EL
     $            +(OE(IREF)-OE(ICHAN))*COS_EL
                  XPTS = (XPTS - XMIN) * XINV
                  YPTS = (YPTS - YMIN) * YINV
                  IF(XPTS.LE.-TOTWX.OR.XPTS.GE.XLIMM) GO TO 8
                  IF(YPTS.LE.-TOTWY.OR.YPTS.GE.YLIMM) GO TO 8
                  SOFFY = YPTS + 1D0
*
                  IFIRST_Y = SOFFY - TOTWY +1D0
                  ILAST_Y  = SOFFY + TOTWY
                  IF (IFIRST_Y.LT.1)    IFIRST_Y = 1
                  IF (ILAST_Y .GT.NY) ILAST_Y  = NY
*
                  DO JJ = IFIRST_Y, ILAST_Y
                     X = (DFLOAT(JJ) - SOFFY) * XINVFY
                     X = ABS(X)
                     ICONO = INT(5D1*X) + 1
                     Z = FUNC(ICONO)
                     ICO = (JJ-1) * NX   ! index of beginning of line
                     SOFFX = XPTS + 1.D0
                     IFIRST_X = SOFFX - TOTWX + 1.0D0
                     ILAST_X  = SOFFX + TOTWX
                     IF (IFIRST_X.LT.1)  IFIRST_X = 1
                     IF (ILAST_X .GT.NX) ILAST_X  = NX
                     DO  KK = IFIRST_X, ILAST_X
                        X = (DFLOAT(KK)-SOFFX) * XINVFX
                        X = ABS(X)
                        ICONO = INT(5D1*X) + 1
                        ZZ = Z * FUNC(ICONO)
*     ZZ=y_weight* x_weight*point_weight
*
                        MAP_DATA(ICO+KK) = MAP_DATA(ICO+KK) +
     $                  VAL*ZZ
                        MAP_WEIGHT(ICO+KK) = MAP_WEIGHT(ICO+KK)
     $                  + ZZ
                     ENDDO
                  ENDDO
               ENDIF
8              CONTINUE
            ENDDO
         ENDIF
      ENDDO
      ELEV_MEAN = ELEV_MEAN/NRECORD
*
*
* --- Compute maximum and minimum weight
*
      CALL NIC_WMAX(MAP_DATA,MAP_WEIGHT,W_MAX,W_MIN,BLANKING,
     $NX*NY)
      IF (W_MAX.LT.1.D-06) THEN
         WRITE(6,*) ' '
         WRITE(6,*) ' *** No point in the selected field. ',
     $   'Mapping 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 NIC_MNORM(MAP_DATA,MAP_WEIGHT,ALIM,ALOW,BLANKING,
     $NX*NY,1.)
      CALL NIC_NORM_WEIGHT(MAP_WEIGHT,NX*NY,IFUN)
*
      RETURN
      END
*
*----------------------------------------------------------------------
*
*    WMAX ()
*
c   this subroutine computes the max and min of a nod2 array, amap,   *
c   containing real*4 numbers in its data section.                    *
*   ----------------------------------------------------------------------
*
      SUBROUTINE NIC_WMAX(MAPD,MAPW,AMAX,AMIN,BLANKING,NPOINT)
*
      REAL*4  MAPD(1),MAPW(1)
      REAL*4  BLANKING
      REAL*8  AMAX,AMIN
      INTEGER NPOINT
      INTEGER I
*
      AMAX = -1.D35
      AMIN =  1.D35
      DO I = 1,NPOINT
         IF(ABS(MAPD(I) - BLANKING).GE. 0.5) THEN
            IF(MAPW(I).GT.AMAX) AMAX = MAPW(I)
            IF(MAPW(I).LT.AMIN) AMIN = MAPW(I)
         ENDIF
      ENDDO
      RETURN
      END
*
* ----------------------------------------------------------------------
      SUBROUTINE NIC_MNORM(MAPD,MAPW,ALIMIT,ALOWLT,BLANKING,NPOINT,CPJY)
*
      REAL*4   MAPD(1),MAPW(1)
      REAL*8   ALIMIT,ALOWLT
      REAL*4   BLANKING, CPJY
      INTEGER  NPOINT,I
*
      DO I=1,NPOINT
         IF (ABS(MAPD(I)-BLANKING) .GE. 0.5) THEN
            IF (MAPW(I).GE.ALIMIT .AND. MAPW(I).GE.ALOWLT) THEN
               MAPD(I) = MAPD(I)/MAPW(I)/CPJY
            ELSE
               MAPD(I) = BLANKING
               MAPW(I) = 0
            ENDIF
         ENDIF
      ENDDO
      RETURN
      END
*
* ------------------------------------------------------------------------
      SUBROUTINE NIC_NORM_WEIGHT(MAPW,NPOINT,TAPER)
*
      REAL*4   MAPW(1)
      INTEGER  NPOINT
      INTEGER  TAPER
      INTEGER  I
*
      IF (TAPER.EQ.0) THEN               ! parabolic
         DO I=1,NPOINT
            MAPW(I) = MAPW(I)*0.15
         ENDDO
      ELSEIF (TAPER.EQ.1) THEN           ! gaussian
         DO I=1,NPOINT
            MAPW(I) = MAPW(I)*0.15
         ENDDO
      ELSEIF (TAPER.EQ.2) THEN           ! linear
         DO I=1,NPOINT
            MAPW(I) = MAPW(I)*0.17
         ENDDO
      ELSEIF (TAPER.EQ.3) THEN           ! no taper
         DO I=1,NPOINT
            MAPW(I) = MAPW(I)*0.20
         ENDDO
      ENDIF
      RETURN
      END
*
*----------------------------------------------------------------------------
      SUBROUTINE REMOVE_PLANE (PLANE,SIGNAL,XCOORD,YCOORD,OA,OE,EL,
     $CHAN_FLAG,CHAN_RMS,IREF,A,B,C,SUPPX,SUPPY,AZ_OFFSET,RMS)
*
      INCLUDE 'const.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      INTEGER NMAX
      PARAMETER (NMAX=2*16+1)
* parameters
      REAL    PLANE(NRECORD*NCHAN), SIGNAL(NRECORD*NCHAN)
      REAL    XCOORD(NRECORD),YCOORD(NRECORD)
      REAL    OA(NCHAN), OE(NCHAN),EL(NRECORD)
      LOGICAL CHAN_FLAG(NCHAN)
      REAL    CHAN_RMS(NCHAN)
      INTEGER IREF
      REAL    A,B,C,AZ_OFFSET
      REAL    SUPPX(NMAX),SUPPY(NMAX)
      REAL    RMS
* local variable
      INTEGER*4 IP_FIT_X, IP_FIT_Y, IP_FIT_Z, IP_FIT_W
      INTEGER*4 POINTER, IP, IADDR_PLANE
      REAL*8  GONS(NMAX-1,4),BOUND(5), XR8, YR8, SUM2
      INTEGER ICHAN, IREC, SIC_GETVM, IER, COUNTER
      REAL    MEAN,VALN
      REAL    VAL,COS_EL,SIN_EL,XPT,YPT,WEI
      INTEGER COUNTER_BLANK,COUNTER_RMS,I
      LOGICAL GR8_IN
*
* Compute parameters for GR8_IN
      DO I=1,NMAX-1
         GONS(I,1) = SUPPX(I)
         GONS(I,2) = SUPPY(I)
         GONS(I,3) = SUPPX(I+1)-SUPPX(I)
         GONS(I,4) = SUPPY(I+1)-SUPPY(I)
      ENDDO
      BOUND(2)= GONS(1,1)
      BOUND(3)= GONS(1,1)
      BOUND(4)= GONS(1,2)
      BOUND(5)= GONS(1,2)
      DO I=2,NMAX-1
         BOUND(2) = MIN(BOUND(2),GONS(I,1))    ! xmin
         BOUND(3) = MAX(BOUND(3),GONS(I,1))    ! xmax
         BOUND(4) = MIN(BOUND(4),GONS(I,2))    ! ymin
         BOUND(5) = MAX(BOUND(5),GONS(I,2))    ! ymax
      ENDDO
*
* Initialize arrays with positions, signal and weight
      IER = SIC_GETVM (NRECORD*NCHAN*4,IADDR_PLANE)
      IF (IER.NE.1) THEN
         CALL MESSAGE (2,1,'SOLVE','Memory allocation failure')
         RETURN
      ENDIF
      IP_FIT_X  = POINTER(IADDR_PLANE,MEMORY)
      IP_FIT_Y  = IP_FIT_X+NRECORD*NCHAN
      IP_FIT_Z  = IP_FIT_Y+NRECORD*NCHAN
      IP_FIT_W  = IP_FIT_Z+NRECORD*NCHAN
      COUNTER = 0
      DO ICHAN =1,NCHAN
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            IF (CHAN_RMS(ICHAN).NE.0) THEN
               WEI = 1/CHAN_RMS(ICHAN)**2
            ELSE
               WEI = 1
            ENDIF
            DO IREC = 1, NRECORD
               VAL = PLANE(IREC+(ICHAN-1)*NRECORD)
               IF (ABS(VAL-BLANKING).GT.0.5) THEN
                  COS_EL = COS(FIELD_ROTATION*EL(IREC)*DEG_TO_RAD)
                  SIN_EL = SIN(FIELD_ROTATION*EL(IREC)*DEG_TO_RAD)
                  XPT = (XCOORD(IREC)-AZ_OFFSET)/10
     $            +(OA(IREF)-OA(ICHAN))*COS_EL
     $            +(OE(IREF)-OE(ICHAN))*SIN_EL
                  YPT = YCOORD(IREC)/10
     $            -(OA(IREF)-OA(ICHAN))*SIN_EL
     $            +(OE(IREF)-OE(ICHAN))*COS_EL
                  IP = IP_FIT_X+COUNTER
                  CALL R4TOR4 (XPT,MEMORY(IP),1)
                  IP = IP_FIT_Y+COUNTER
                  CALL R4TOR4 (YPT,MEMORY(IP),1)
                  IP = IP_FIT_Z+COUNTER
                  CALL R4TOR4 (VAL,MEMORY(IP),1)
                  IP = IP_FIT_W+COUNTER
                  CALL R4TOR4 (WEI,MEMORY(IP),1)
                  COUNTER = COUNTER + 1
               ENDIF
            ENDDO
         ENDIF
      ENDDO
* Fit plane
      IF (COUNTER.NE.0) THEN
         CALL FIT_PLANE (MEMORY(IP_FIT_X),MEMORY(IP_FIT_Y),
     $   MEMORY(IP_FIT_Z),MEMORY(IP_FIT_W),COUNTER,A,B,C,
     $   SUPPX,SUPPY)
         COUNTER = 0
         COUNTER_BLANK=0
         COUNTER_RMS = 0
         SUM2 = 0
         DO ICHAN =1,NCHAN
            IF (.NOT.CHAN_FLAG(ICHAN)) THEN
               MEAN = 0
               VALN = 0
               DO IREC = 1, NRECORD
                  VAL = PLANE(IREC+(ICHAN-1)*NRECORD)
                  IF (ABS(VAL-BLANKING).GT.0.5) THEN
                     IP = IP_FIT_Z+COUNTER
                     CALL R4TOR4 (MEMORY(IP),VAL,1)
                     PLANE(IREC+(ICHAN-1)*NRECORD) = VAL
                     MEAN = MEAN+VAL
                     VALN = VALN+1
                     COUNTER = COUNTER + 1
                  ENDIF
               ENDDO
               DO IREC = 1, NRECORD
                  VAL = SIGNAL(IREC+(ICHAN-1)*NRECORD)
                  IF (ABS(VAL-BLANKING).GT.0.5) THEN
                     SIGNAL(IREC+(ICHAN-1)*NRECORD) =
     $               SIGNAL(IREC+(ICHAN-1)*NRECORD) - MEAN/VALN
                  ENDIF
               ENDDO
               DO IREC = 1, NRECORD
                  VAL = SIGNAL(IREC+(ICHAN-1)*NRECORD)
                  IF (ABS(VAL-BLANKING).GT.0.5) THEN
                     IP = IP_FIT_X+COUNTER_BLANK
                     CALL R4TOR8(MEMORY(IP),XR8,1)
                     IP = IP_FIT_Y+COUNTER_BLANK
                     CALL R4TOR8(MEMORY(IP),YR8,1)
                     IF (.NOT.GR8_IN(XR8,YR8,NMAX-1,NMAX-1,GONS,BOUND))
     $               THEN
                        SUM2 = VAL**2
                        COUNTER_RMS=COUNTER_RMS+1
                     ENDIF
                     COUNTER_BLANK = COUNTER_BLANK  + 1
                  ENDIF
               ENDDO
            ENDIF
         ENDDO
      ENDIF
      IF (COUNTER_RMS.NE.0) THEN
         RMS = SQRT(SUM2/COUNTER_RMS)
      ELSE
         RMS=0
      ENDIF
      CALL FREE_VM (NRECORD*NCHAN*4,IADDR_PLANE)
      RETURN
      END
*
*----------------------------------------------------------------------------
      SUBROUTINE FIT_PLANE (X,Y,Z,W,N,A,B,C,SUPPX,SUPPY)
*
      INTEGER NMAX
      PARAMETER (NMAX=2*16+1)
      INTEGER N
      REAL    A,B,C
      REAL    X(N),Y(N),Z(N),W(N)
      REAL    SUPPX(NMAX), SUPPY(NMAX)
*
* local variables
      REAL*8  G,GX,GY,GZ,GXX,GYY,GXY,GZX,GZY
      REAL*8  GONS(NMAX-1,4),BOUND(5), XR8, YR8, ZR8, WR8
      REAL    D
      INTEGER I
      LOGICAL GR8_IN
*
      DO I=1,NMAX-1
         GONS(I,1) = SUPPX(I)
         GONS(I,2) = SUPPY(I)
         GONS(I,3) = SUPPX(I+1)-SUPPX(I)
         GONS(I,4) = SUPPY(I+1)-SUPPY(I)
      ENDDO
      BOUND(2)= GONS(1,1)
      BOUND(3)= GONS(1,1)
      BOUND(4)= GONS(1,2)
      BOUND(5)= GONS(1,2)
      DO I=2,NMAX-1
         BOUND(2) = MIN(BOUND(2),GONS(I,1))    ! xmin
         BOUND(3) = MAX(BOUND(3),GONS(I,1))    ! xmax
         BOUND(4) = MIN(BOUND(4),GONS(I,2))    ! ymin
         BOUND(5) = MAX(BOUND(5),GONS(I,2))    ! ymax
      ENDDO
      G   = 0
      GX  = 0
      GY  = 0
      GZ  = 0
      GXX = 0
      GYY = 0
      GXY = 0
      GZX = 0
      GZY = 0
      DO I=1,N
         XR8 = X(I)
         YR8 = Y(I)
         ZR8 = Z(I)
         WR8 = W(I)
         IF (.NOT.GR8_IN(XR8,YR8,NMAX-1,NMAX-1,GONS,BOUND)) THEN
            G   = G   + WR8
            GX  = GX  + WR8*XR8
            GY  = GY  + WR8*YR8
            GZ  = GZ  + WR8*ZR8
            GXX = GXX + WR8*XR8*XR8
            GYY = GYY + WR8*YR8*YR8
            GXY = GXY + WR8*XR8*YR8
            GZX = GZX + WR8*ZR8*XR8
            GZY = GZY + WR8*ZR8*YR8
         ENDIF
      ENDDO
      D = G*(GXX*GYY-GXY*GXY)+GX*(GY*GXY-GX*GYY)+GY*(GX*GXY-GY*GXX)
      IF (D.NE.0) THEN
         A = GZ*(GY*GXY-GX*GYY)+GZX*(G*GYY-GY*GY)+GZY*(GX*GY-G*GXY)
         B = GZ*(GX*GXY-GY*GXX)+GZX*(GX*GY-G*GXY)+GZY*(G*GXX-GX*GX)
         C = GZ*(GXX*GYY-GXY*GXY)+GZX*(GY*GXY-GX*GYY)+
     $   GZY*(GX*GXY-GY*GXX)
         A = A/D
         B = B/D
         C = C/D
         DO I=1,N
            XR8 = X(I)
            YR8 = Y(I)
            Z(I) = (A*X(I)+B*Y(I)+C)
         ENDDO
      ENDIF
      RETURN
      END
*----------------------------------------------------------------------------
      SUBROUTINE COMPUTE_OFFSET (SCAN_AZ, AZ_OFFSET)
      INCLUDE 'nic.inc'
      REAL*4  SCAN_AZ(NRECORD), AZ_OFFSET
      AZ_OFFSET = SCAN_AZ(NRECORD)
      RETURN
      END
*----------------------------------------------------------------------------
      SUBROUTINE BUILD_MDB_IMAGE(MAP_D,NX,NY,XMIN,YMIN,XINC,YINC,
     $A,B,C,REMOVE)
*
*      make GILDAS image file
*----------------------------------------------------------------------------
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:format.inc'
*
      REAL    MAP_D(1),XMIN,YMIN,XINC,YINC
      INTEGER NX,NY
      REAL    A,B,C
      LOGICAL REMOVE
*  local variables
      INTEGER*4 IP,POINTER
      INTEGER I,J,II
      REAL    MAXMIN(6),VAL
      LOGICAL ERROR
      CHARACTER*12 SNAM
      INTEGER   BLC(4), TRC(4)
      DATA   BLC /0,0,0,0/, TRC /0,0,0,0/
*
* Remove plane (if required)
      IF (REMOVE) THEN
         DO I=1,NX
            DO J=1,NY
               II = I + (J-1)*NX
               IF (ABS(MAP_D(II)-BLANKING).GT.0.5) THEN
                  MAP_D(II) = MAP_D(II)-
     $            A*(XMIN+(I-1)*XINC)-B*(YMIN+(J-1)*YINC)-C
               ENDIF
            ENDDO
         ENDDO
      ENDIF
* Find min max
      MAXMIN(1) =  1.D31                 ! MAXMIN(1) = minimum
      MAXMIN(2) = -1.D31                 ! MAXMIN(2) = maximum
      DO J=1,NY
         DO I=1,NX
            II = I + (J-1)*NX
            VAL = MAP_D(II)
            IF (ABS(VAL-BLANKING).GT.0.5)
     $      MAXMIN(1) = MIN (MAXMIN(1),VAL)
            IF (MAXMIN(1).EQ.VAL) THEN
               MAXMIN(3) = I
               MAXMIN(4) = J
            ENDIF
            MAXMIN(2) = MAX (MAXMIN(2),VAL)
            IF (MAXMIN(2).EQ.VAL) THEN
               MAXMIN(5) = I
               MAXMIN(6) = J
            ENDIF
         ENDDO
      ENDDO
* Prepare header
      X_NDIM   = 2
      IF (.NOT.JANSKY_FLAG) THEN
         X_UNIT   = ' Counts'
      ELSE
         X_UNIT   = ' Jy'
      ENDIF
      X_DIM(1) = NX
      X_DIM(2) = NY
      X_DIM(3) = 1
      X_DIM(4) = 1
      X_SIZE   = X_DIM(1)*X_DIM(2)*X_DIM(3)*X_DIM(4)
*
      X_REF1 = 1
      X_VAL1 = -ABS(XMIN)*SEC_TO_RAD
      X_INC1 = XINC*SEC_TO_RAD
*
      X_REF2 = 1
      X_VAL2 = -ABS(YMIN)*SEC_TO_RAD
      X_INC2 = YINC*SEC_TO_RAD
*
      X_REF3 = 0
      X_VAL3 = 0
      X_INC3 = 1
*
      X_REF4 = 0
      X_VAL4 = 0
      X_INC4 = 1
*
      X_TYPE = 'GILDAS_IMAGE'
      X_NAME = SOURCE_NAME
      X_FORM = FMT_R4
      X_GENE = 29
      X_SYST = 'AZIMUTHAL'
      X_BLAN = 2
      X_BVAL = BLANKING
      X_EVAL = 1
      X_EXTR = 10
      X_RMIN = MAXMIN(1)
      X_RMAX = MAXMIN(2)
      X_MIN1 = MAXMIN(3)
      X_MIN2 = MAXMIN(4)
      X_MIN3 = 1
      X_MAX1 = MAXMIN(5)
      X_MAX2 = MAXMIN(6)
      X_MAX3 = 1
      X_MIN4 = 1
      X_MAX4 = 1
*
      X_DESC = 18
      X_POSI = 12
      X_RA   = SLAM*DEG_TO_RAD
      X_DEC  = SBET*DEG_TO_RAD
      X_PROJ = 9
      X_SPEC = 12
      X_RESO = 3
      X_FREQ = FREQUENCY*1.D3            ! frequency in Mhz
*
      SNAM = ' '
      WRITE (SNAM(1:4),'(i4.4)') SCAN_NUMBER
      CALL SIC_LOWER (SNAM)
      CALL SIC_PARSEF(SNAM,X_FILE,' ','.mdb')
      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(MAP_D(1),MEMORY(IP),NX*NY)
      CALL GDF_FRIS (X_ISLO,ERROR)
*
      RETURN
      END
*
      SUBROUTINE  DISPLAY_CHANNEL_POSITION(OA,OE,N,ELEV,IREF,
     $FIELD_ROTATION)
*
      INCLUDE 'const.inc'
*
      INTEGER      N,IREF,FIELD_ROTATION
      REAL         OA(N),OE(N)
      REAL         ELEV
* local variables
      INTEGER      ICHAN
      REAL         COO(2)
      REAL         COS_EL,SIN_EL
      CHARACTER*80 STRING
*
      WRITE (STRING,*)'ellipse /user'
      COS_EL = COS(FIELD_ROTATION*ELEV/180.*PI)
      SIN_EL = SIN(FIELD_ROTATION*ELEV/180.*PI)
*
      CALL GR_EXEC('set exp 0.5')
      DO ICHAN=1,N
         COO(1) =  (OA(IREF)-OA(ICHAN))*COS_EL
     $   +(OE(IREF)-OE(ICHAN))*SIN_EL
         COO(2) = -(OA(IREF)-OA(ICHAN))*SIN_EL
     $   +(OE(IREF)-OE(ICHAN))*COS_EL
         WRITE (STRING,'(a,2f8.3)')'ellipse 5. /user', COO(1),COO(2)
         CALL GR_EXEC(STRING)
         IF (ICHAN.LE.9) THEN
            WRITE (STRING,'(a,2f8.3,a,i1,a)')'draw text ',
     $      COO(1),COO(2),' "',ICHAN,'" 5 /user'
         ELSE
            WRITE (STRING,'(a,2f8.3,a,i2,a)')'draw text ',
     $      COO(1),COO(2),' "',ICHAN,'" 5 /user'
         ENDIF
         CALL GR_EXEC(STRING)
      ENDDO
      CALL GR_EXEC('set exp 1.')
      END