* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
*
*  Computes pointing maps from rasterdips
*
      SUBROUTINE SOLVE_RASTERDIP (LINE,ERROR)
*
*----------------------------------------------------------------------
*       Support routine for command
*       SOLVE
* Arguments :
*       LINE    C*(*)   Command line                    Input
*       ERROR   L       Logical error flag              Output
*
*----------------------------------------------------------------------
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      LOGICAL ERROR,SUBSCAN_OK
      CHARACTER*(*) LINE
      INTEGER*4 POINTER, IP_SUBFLAG
      INTEGER*4 SIGNAL_PLANE, IP_PLANE
      INTEGER*4 IP, IP_REFSIGNAL, SIC_GETVM
      INTEGER  IER, I, IEFF, DIM(4), NRECS
      REAL     XMIN,XMAX,YMIN,YMAX,PLANE(3),X4(2),Y4(2),AZ,EL
      REAL     SUBELEV,XINC,YINC,CUTOFF,ELEV_MEAN,EXTVAL(2),A,B,C
      REAL     SINA, COSA, TANB
      REAL     PARAM_FIT(60)             ! 15 x 4
* [1] : line orientation (degrees)
* [2] : line offset (arcseconds)
* [3] : slope (counts/arcseconds)
* [4] : elevation (degrees)
      INTEGER  NX,NY,IM,IN,IFUN,IREF
      CHARACTER LIMIT_STRING*40, SETBX_STRING*40, LABEL_STRING*40
      CHARACTER GREEK*2, ITALIC*3, TEXT*60
      DATA X4/ -100.,100./
      SAVE PLANE
      SAVE PARAM_FIT
*
      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
*
      IF (NSUBSCAN.LE.1) THEN
         CALL MESSAGE (2,1,'SOLVE','Bad number of subscans ')
         RETURN
      ENDIF
*
      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)
      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 GR_EXEC ('set plot port')
      IREF = 1
      CALL COMPUTE_MAXMIN(MEMORY(IP_SIGNAL),NCHAN*NRECORD,EXTVAL)
      IP_SUBFLAG  = POINTER(IADDR_SUBFLAG,MEMORY)
      IEFF = 0
      DO I=1,NSUBSCAN
         IF (SUBSCAN_OK (MEMORY(IP_SUBFLAG),I)) IEFF = IEFF+1
      ENDDO
      DO I=1,NSUBSCAN
         IF (SUBSCAN_OK (MEMORY(IP_SUBFLAG),I))THEN
            CALL SMEAR (MEMORY(IP_SCAN_COORD),
     $      MEMORY(IP_SCAN_COORD+NRECORD),MEMORY(IPSUBREC),I)
            CALL COMPUTE_LIMITS (MEMORY(IP_SCAN_COORD),
     $      MEMORY(IP_SCAN_COORD+NRECORD),
     $      MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_EL),IREF,
     $      XMIN,XMAX,YMIN,YMAX,MEMORY(IPSUBREC),I,NRECS,SUBELEV,
     $      AZ,MEMORY(IP_AZ))
*
            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*NRECS*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*NRECS
            CALL R4TOR4(MEMORY(IP_SIGNAL),MEMORY(IP_REFSIGNAL),
     $      NCHAN*NRECS)
*
            CALL NIC_FUNGEN(FUNC,IFUN)
            LIMIT_STRING = 'LIM -000.000 +000.000 -000.000 +000.000'
            WRITE (LIMIT_STRING(05:12),'(f8.3)') XMIN
            WRITE (LIMIT_STRING(14:21),'(f8.3)') XMIN+NX*XINC
            WRITE (LIMIT_STRING(23:30),'(f8.3)') YMIN
            WRITE (LIMIT_STRING(32:39),'(f8.3)') YMIN+NY*YINC
            IM = MOD(I-1,4)              ! number of columns
            IN = INT((I-1)/4)            ! number of lines
            SETBX_STRING = 'SET BOX 000.000 000.000 000.000 000.000'
            WRITE (SETBX_STRING(09:15),'(f7.3)')  2.5+4.5*IM
            WRITE (SETBX_STRING(17:23),'(f7.3)')  2.5+4.5*(IM+1)
            WRITE (SETBX_STRING(25:31),'(f7.3)')  27.-4.5*(IN+1)
            WRITE (SETBX_STRING(33:39),'(f7.3)')  27.-4.5*IN
            CALL COMPUTE_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,MEMORY(IPSUBREC),I,MEMORY(IP_RMS),
     $      A,B,C)
            CALL GR_EXEC ('LIM /RG POINTING_MAP')
            CALL GR_EXEC (SETBX_STRING)
            CALL GR_EXEC ('SET BLA -99999.5 1')
            WRITE (TEXT,'(a,2f9.1)') 'plot /scal lin ',
     $      EXTVAL(1)/1.3,EXTVAL(2)/1.3
            CALL GR_EXEC (TEXT)
            CALL GR_EXEC ('SET EXP 1')
            CALL GR_EXEC (LIMIT_STRING)
            IF (I.EQ.1) THEN
               CALL SIC_DELVARIABLE('PARAM_POINTING',.FALSE.,ERROR)
               DIM(1) = IEFF
               DIM(2) = 4
               CALL SIC_DEF_REAL ('PARAM_POINTING',PARAM_FIT,2,
     $         DIM,.TRUE.,ERROR)
               CALL GR_EXEC ('set exp 1.')
               WRITE (TEXT,'(a9,2f8.3,a11,i5,a9)') 'draw tex ',
     $         XMIN,YMAX*1.05,' "\\iSCAN: ',SCAN_NUMBER,'" 9 /user'
               CALL GR_EXEC (TEXT)
               WRITE (TEXT,'(a9,2f8.3,a14,i4,a9)') 'draw tex ',
     $         XMAX,YMAX*1.05,' "\\iAZIMUTH: ',INT(AZ),
     $         '" 9 /user'
               CALL GR_EXEC (TEXT)
            ENDIF
            IF (A.NE.0) THEN
               PARAM_FIT(I)        = ATAN(B/A)/PI*180.     ! angle between vertical and line
               PARAM_FIT(I+IEFF)   = -C/A
               PARAM_FIT(I+IEFF*2) = SQRT(A**2+B**2)
               PARAM_FIT(I+IEFF*3) = SUBELEV
            ELSE
               PARAM_FIT(I) = 0
               PARAM_FIT(I+IEFF) = 0
               PARAM_FIT(I+IEFF*2) = 0
               PARAM_FIT(I+IEFF*3) = 0
            ENDIF
            IF (B.NE.0) THEN
               Y4(1) =  -(C+X4(1)*A)/B
               Y4(2) =  -(C+X4(2)*A)/B
               CALL GR_SEGM ('plot',ERROR)
               CALL GR4_CONNECT (2,X4,Y4,0.,-1.)
               CALL GR_OUT
            ENDIF
            CALL GR_EXEC ('set exp .8')
            WRITE (TEXT,'(a,2f10.3,a,i2,a)') 'draw text ',
     $      XMAX-(XMAX-XMIN)/4,YMIN+(YMAX-YMIN)/9,
     $      ' "\\iEL=',NINT(SUBELEV),'" 5 /user'
            CALL GR_EXEC (TEXT)
            CALL GR_EXEC ('ellipse 5 /user 0 0')
            CALL GR_EXEC ('BOX N N')
            CALL SIC_DELVARIABLE('POINTING_MAP',.FALSE.,ERROR)
            CALL R4TOR4 (MEMORY(IP_SIGNAL),MEMORY(IP_PLANE),
     $      NRECS*NCHAN)
         ENDIF
      ENDDO
      WRITE (TEXT,'(a,2f9.1)') 'wedge top .3 /scal lin ',
     $EXTVAL(1),EXTVAL(2)
      CALL GR_EXEC ('set box 2.5 20.5 28 28.2')
      CALL GR_EXEC (TEXT)
      CALL GR_EXEC ('set exp 1.')
!
! plot results from fit
!
      CALL GR_EXEC ('set box 9 14 3 8')
      CALL GR_EXEC ('box')
      DO I=1,NCHAN
         CALL GETPOS (MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),I,AZ,EL)
         WRITE (TEXT,'(a,2f10.2)') 'ellipse 5.5 /user ',AZ,EL
         CALL GR_EXEC (TEXT)
      ENDDO
      CALL GR_EXEC ('pen 1')
      CALL GR_SEGM ('plot',ERROR)
      DO I=1,IEFF
         COSA = COS((PARAM_FIT(I+IEFF*3)+180)*DEG_TO_RAD)
         SINA = SIN((PARAM_FIT(I+IEFF*3)+180)*DEG_TO_RAD)
         TANB = TAN((PARAM_FIT(I+IEFF*3)+PARAM_FIT(I)+270)*DEG_TO_RAD)
         Y4(1) = X4(1)*TANB+PARAM_FIT(I+IEFF)*(SINA-TANB*COSA)
         Y4(2) = X4(2)*TANB+PARAM_FIT(I+IEFF)*(SINA-TANB*COSA)
         CALL GR4_CONNECT(2,X4,Y4,0.,-1.)
      ENDDO
      CALL GR_OUT
      CALL GR_EXEC ('pen 0')
      LABEL_STRING =
     $'LABEL "'//ITALIC//'Nasmyth '//GREEK//'DVe (arcsec)" /Y'
      CALL GR_EXEC (LABEL_STRING)
      CALL GR_EXEC ('set box 9 14 2.5 8')
      LABEL_STRING =
     $'LABEL "'//ITALIC//'Nasmyth '//GREEK//'DHo (arcsec)" /X'
      CALL GR_EXEC (LABEL_STRING)
      RETURN
      END
      SUBROUTINE SMEAR (AZ_OUTPUT,EL_OUTPUT,NSUBREC,NSUB)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      INTEGER*4 NSUBREC(NSUBSCAN), NSUB, I, L
      REAL*4    AZ_OUTPUT(NRECORD), EL_OUTPUT(NRECORD)
*
* Mimick pointing scans /length 60
*
      L = NSUBREC(NSUB)
      DO I=1,L/2
         AZ_OUTPUT(I) = 0                ! 300*(2*(i-1.)/(l/2-1)-1)
         EL_OUTPUT(I) = 0
      ENDDO
      DO I=L/2+1,L
         AZ_OUTPUT(I) = 0
         EL_OUTPUT(I) = 0                ! 300*(2*(i-l/2-1.)/(l-l/2-1)-1)
      ENDDO
*
      RETURN
      END
      SUBROUTINE COMPUTE_LIMITS(XCOORD,YCOORD,OA,OE,EL,IREF,
     $XMIN,XMAX,YMIN,YMAX,NSUBREC,I,NRECS,SUBELEV,AZI,AZ)
*
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      REAL    XCOORD(1),YCOORD(1)
      REAL    OA(1),OE(1),AZ(1),EL(1),FC,AZI
      INTEGER IREF                       ! reference channel
      REAL    XMIN,XMAX,YMIN,YMAX
*
      INTEGER IREC,ICHAN,NSUBREC(NSUBSCAN),I,IOLD,INEW,NRECS
      REAL    COS_EL,SIN_EL,SUBELEV
      REAL    XCOO,YCOO
*
      NRECS = NSUBREC(I)
*
      XMIN = 1.D32
      XMAX = -1.D32
      YMIN = 1.D32
      YMAX = -1.D32
      FC = 1.3
*
      INEW = 0
      DO IREC=1,I
         IOLD = INEW
         INEW = IOLD+NSUBREC(IREC)
      ENDDO
      IOLD = IOLD+1
*
      DO IREC=IOLD,INEW
         COS_EL = COS(FIELD_ROTATION*EL(IREC)*DEG_TO_RAD)
         SIN_EL = SIN(FIELD_ROTATION*EL(IREC)*DEG_TO_RAD)
         DO ICHAN = 1,NCHAN
            XCOO = XCOORD(IREC-IOLD+1)/10.
     $      +(OA(IREF)-OA(ICHAN))*COS_EL
     $      +(OE(IREF)-OE(ICHAN))*SIN_EL
            YCOO = YCOORD(IREC-IOLD+1)/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
      IF (-XMIN.GT.XMAX) THEN
         XMAX = -XMIN
      ELSE
         XMIN = -XMAX
      ENDIF
      IF (-YMIN.GT.YMAX) THEN
         YMAX = -YMIN
      ELSE
         YMIN = -YMAX
      ENDIF
      IF (YMAX.GT.XMAX) THEN
         XMAX = YMAX
         XMIN = -YMAX
      ELSE
         YMAX = XMAX
         YMIN = -XMAX
      ENDIF
      XMIN = XMIN*FC
      XMAX = XMAX*FC
      YMIN = YMIN*FC
      YMAX = YMAX*FC
      SUBELEV = NINT(EL(IOLD))
      AZI = AZ(1)
      END
      SUBROUTINE  COMPUTE_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,NSUBREC,ISUB,
     $CHAN_RMS,A,B,C)
*
      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*8  W_MIN,W_MAX
      REAL*8  ALIM,ALOW
      INTEGER*4 IP, IP_FIT_X, IP_FIT_Y, IP_FIT_Z, IP_FIT_W
      INTEGER*4 IADDR_PLANE, SIC_GETVM, POINTER
      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)
      REAL    A,B,C
      INTEGER NX,NY,IFUN,IREF,NSUBREC(NSUBSCAN),IER
      REAL    XINC,YINC,XMIN,YMIN
      REAL    CUTOFF, ELEV_MEAN,CHAN_RMS(NCHAN)
      LOGICAL ERROR,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,WEI
      INTEGER IREC,ICHAN,I,COUNTER
      INTEGER IFIRST_X,ILAST_X,IFIRST_Y,ILAST_Y
      INTEGER JJ,ICONO,ICO,KK,ISUB,IOLD,INEW
      REAL    XPTS,YPTS
      REAL    SOFFY,SOFFX,X
      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
*
* Initialize arrays with positions, signal and weight
*
      IER = SIC_GETVM (NSUBREC(ISUB)*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+NSUBREC(ISUB)*NCHAN
      IP_FIT_Z  = IP_FIT_Y+NSUBREC(ISUB)*NCHAN
      IP_FIT_W  = IP_FIT_Z+NSUBREC(ISUB)*NCHAN
      COUNTER = 0
*
* Add scan data to map_data and map_weight
*
      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
*
            ELEV_MEAN = 0
*
            INEW = 0
            DO IREC=1,ISUB
               IOLD = INEW
               INEW = IOLD+NSUBREC(IREC)
            ENDDO
            IOLD = IOLD+1
            DO IREC = IOLD,INEW
               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-IOLD+1)/10
     $            +(OA(IREF)-OA(ICHAN))*COS_EL
     $            +(OE(IREF)-OE(ICHAN))*SIN_EL
                  YPTS = YCOORD(IREC-IOLD+1)/10
     $            -(OA(IREF)-OA(ICHAN))*SIN_EL
     $            +(OE(IREF)-OE(ICHAN))*COS_EL
*
                  IP = IP_FIT_X+COUNTER
                  CALL R4TOR4 (XPTS,MEMORY(IP),1)
                  IP = IP_FIT_Y+COUNTER
                  CALL R4TOR4 (YPTS,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
*
                  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
      CALL COMPUTE_FIT (MEMORY(IP_FIT_X),MEMORY(IP_FIT_Y),
     $MEMORY(IP_FIT_Z),MEMORY(IP_FIT_W),COUNTER,A,B,C)
      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)
      CALL FREE_VM (NSUBREC(ISUB)*NCHAN*4,IADDR_PLANE)
*
      RETURN
      END
*
*
      SUBROUTINE COMPUTE_FIT (X,Y,Z,W,N,A,B,C)
*
      INTEGER N
      REAL    A,B,C
      REAL    X(N),Y(N),Z(N),W(N)
*
* local variables
      REAL*8  G,GX,GY,GZ,GXX,GYY,GXY,GZX,GZY
      REAL*8  XR8, YR8, ZR8, WR8
      REAL    D
      INTEGER I
*
      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)
         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
      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
      ENDIF
      RETURN
      END
*
      SUBROUTINE COMPUTE_MAXMIN (SIGNAL,NUM,VAL)
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      INTEGER I, NUM
      REAL SIGNAL(NUM), VALUE, VAL(2)
*
      VAL(1) =  1E34
      VAL(2) = -1E34
      DO I=1,NUM
         VALUE = SIGNAL(I)
         IF (ABS(VALUE-BLANKING).GT.0.5) THEN
            IF (VAL(1).GT.VALUE) VAL(1) = VALUE
            IF (VAL(2).LT.VALUE) VAL(2) = VALUE
         ENDIF
      ENDDO
      END
*
      FUNCTION SUBSCAN_OK (SUBSCAN_FLAG,NSUB)
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      INTEGER*4 SUBSCAN_FLAG(1)
      LOGICAL SUBSCAN_OK
      INTEGER I,K,NSUB
*
      SUBSCAN_OK = .TRUE.
      DO I=1,NCHAN
         K = (I-1)*NSUBSCAN+1
         DO WHILE (SUBSCAN_FLAG(K).EQ.0 .AND. K.LE.I*NSUBSCAN)
            K = K+1
         ENDDO
         IF (K.LE.I*NSUBSCAN) THEN       ! flags found
            DO K=(I-1)*NSUBSCAN+1,I*NSUBSCAN
               IF (SUBSCAN_FLAG(K).NE.0.AND.K-(I-1)*NSUBSCAN.EQ.NSUB)
     $         SUBSCAN_OK = .FALSE.
            ENDDO
         ENDIF
      ENDDO
      RETURN
      END
!
      SUBROUTINE GETPOS (INP_AZ,INP_EL,CHAN,AZ,EL)
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
!
      INTEGER CHAN
      REAL INP_AZ(NCHAN), INP_EL(NCHAN), AZ,EL
!
      AZ = INP_AZ(CHAN)
      EL = INP_EL(CHAN)
      RETURN
      END
