* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      SUBROUTINE NIC_FIT (LINE,ERROR)
*----------------------------------------------------------------------
*       Support routine for command
*       FIT /INPUT filename /FUNCT [1|2] [param1 param2 param3 param4]
*           /FILE  filename [NEW]
* Arguments :
*       LINE    C*(*)   Command line                    Input
*       ERROR   L       Logical error flag              Output
*
*----------------------------------------------------------------------
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
      CHARACTER*(*) LINE
      LOGICAL  ERROR
*
      INTEGER*4 IPMAP,IADDR,IPDATA
      INTEGER*4 POINTER
*
      REAL*8    FSUMSQ,XPARAM(8),XPARAM0(8),AMPLITUDE,WIDTH,ANGLE
      REAL*8    XPARAM_MIN(8),FSUMSQ_MIN,WIDTH_DEF
      REAL*8    TEMP(6)
      REAL*8    REF1, INC1, VAL1, REF2, INC2, VAL2
      REAL*8    Y, DY(8)
      REAL*8    TOL, DENORM
      INTEGER*4 IP_AZ_OFF, IP_EL_OFF, IP_AMP_POS, IP_AMP_NEG
      INTEGER*4 IP_AZ_POS, IP_WIDTH, IP_WIDTH1, IP_WIDTH2
      INTEGER*4 IP_IW, IP_IWA, IPFVEC
*
      INTEGER   LWORK,MPAR,NDATA,IMAX,IMIN,IFUNC
      INTEGER*4 NX, NY,ICHAN,FIRST_CHAN,LAST_CHAN,NC,I,J
      INTEGER*4 IER, SIC_GETVM
      INTEGER*4 SIC_GETLUN,SIC_FRELUN, LUN(2), NLUN, SIC_OPEN
      INTEGER*4 NITE
      REAL      LEVEL, MEAN_ELE, UT, LST, RES(NREC), DUMMY(NREC)
      LOGICAL   SIC_PRESENT,PARAM_OPT,INPUT_OPT,LEVEL_OPT
      LOGICAL   OUTPUT_OPT
      LOGICAL   GAUSS_2D
      LOGICAL   EKH_SAA_TYPE
      LOGICAL   SIC_CTRLC
      INTEGER   BLC(4), TRC(4), LENC
      INTEGER   LWA,IOPT,NPRINT,INFO,INFO_SOL
      CHARACTER SNAM*40,STRING*40,TYPE*40,FILENAME*40, MONTH(12)*3
      CHARACTER INPUT_LINE*300, STRING_100*80, STRING_110*80
      CHARACTER STRING_101*80, STRING_111*80
      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/
      COMMON/PARAM_MIN/XPARAM_MIN,FSUMSQ_MIN,INFO_SOL
      DATA BLC /0,0,0,0/, TRC /0,0,0,0/
      DATA MONTH /'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG',
     $'SEP','OCT','NOV','DEC'/
      DATA WIDTH_DEF /11.D0/             ! default value for width
      EXTERNAL FCN
*
      CALL XSETF(0)
      FSUMSQ_MIN=1.E32
*
* Read file xxxx.mdb, or xxxx.ekh or xxxx.saa or xxxx.cnv according to
* variable TYPE
* If option /INPUT filename, open filename
*
*
      INPUT_OPT = SIC_PRESENT(2,0)
      IF (INPUT_OPT) THEN
         CALL SIC_CH (LINE,2,1,X_FILE,NC,.TRUE.,ERROR)
         IF (ERROR) RETURN
         CALL SIC_LOWER(X_FILE)
         I = INDEX(X_FILE,'.')
         TYPE = X_FILE(I+1:I+3)
      ELSE
         SNAM = ' '
         WRITE (SNAM(1:4),'(i4.4)') SCAN_NUMBER
         CALL SIC_GET_CHAR ('TYPE',TYPE,NC,ERROR)
         CALL SIC_LOWER(TYPE)
         CALL SIC_PARSEF(SNAM,X_FILE,' ','.'//TYPE(1:3))
      ENDIF
      CALL GDF_GEIS  (X_ISLO,ERROR)      ! reserve slot
      CALL GDF_REIS  (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
      IF (ERROR) THEN
         IF (INPUT_OPT) THEN
            WRITE (STRING,'(a,1x,a)') 'Error on file ',X_FILE(1:NC)
         ELSE
            WRITE (STRING,'(a,1x,a)') 'Error on file ',
     $      SNAM(1:4)//'.'//TYPE(1:3)
         ENDIF
         CALL MESSAGE (2,1,'FIT',STRING)
         GOTO 99
      ENDIF
      CALL GDF_READX  (X_ISLO,ERROR)
      CALL GDF_GEMS  (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
      IPMAP = POINTER (X_ADDR,MEMORY)
*
* Time estimator
*
      CALL R4TOR4 (MEMORY(IPSUBSCAN+2*NSUBSCAN),LST,1)
      UT = LST
      CALL STUT (UT,1,IDAY,IMONTH,IYEAR,SITE_LONGITUDE/DEG_TO_RAD)
*
      DOUBLEBEAM_FIT = .FALSE.
      CNV_MAP    = .FALSE.
      EKH_SAA_TYPE   = .FALSE.
      IF (TYPE(1:3).EQ.'mdb') DOUBLEBEAM_FIT = .TRUE.
      IF (TYPE(1:3).EQ.'cnv') CNV_MAP    = .TRUE.
      IF (TYPE(1:3).EQ.'ekh'.OR.TYPE(1:3).EQ.'saa')
     $EKH_SAA_TYPE   = .TRUE.
      NX = X_DIM(1)
      NY = X_DIM(2)
      NCHAN_FIT = X_DIM(3)
      IF (CNV_MAP) NCHAN_FIT = 1
      REF1 = X_REF1
      REF2 = X_REF2
      IF (X_INC1.NE.1) THEN
         VAL1 = X_VAL1/SEC_TO_RAD
         INC1 = X_INC1/SEC_TO_RAD
         VAL2 = X_VAL2/SEC_TO_RAD
         INC2 = X_INC2/SEC_TO_RAD
      ELSE
         VAL1 = X_VAL1
         INC1 = X_INC1
         VAL2 = X_VAL2
         INC2 = X_INC2
      ENDIF
*   Test if computing of offsets is possible (elevation must be available)
      COMPUTE_OFF_FIT = .FALSE.
      IF (DOUBLEBEAM_FIT.AND..NOT.INPUT_OPT) COMPUTE_OFF_FIT = .TRUE.
      IF (EKH_SAA_TYPE  .AND..NOT.INPUT_OPT) COMPUTE_OFF_FIT = .TRUE.
      IF ((DOUBLEBEAM_FIT.OR.EKH_SAA_TYPE).AND.INPUT_OPT
     $.AND.IADDRRECORD.NE.0) THEN
         WRITE(SNAM,'(i4.4)') SCAN_NUMBER
         IF (SNAM(1:4).EQ.X_FILE(1:4)) COMPUTE_OFF_FIT=.TRUE.
      ENDIF
*
      STRING_100 = '! Chan  Az(arcsec) El(arcsec)  A(counts)'//
     $'   FWHM(arcsec) El(degree) sigma info'
      STRING_110 = '! Plane RA(arcsec) Dec(arcsec) A(counts)'//
     $'   FWHM(arcsec) El(degree) sigma info'
      STRING_101 = '! Chan  Az(arcsec) El(arcsec)  A(counts)'//
     $'   FWHM(arcsec)   PAng(deg) sigma info'
      STRING_111 = '! Plane RA(arcsec) Dec(arcsec) A(counts)'//
     $'   FWHM(arcsec)   PAng(deg) sigma info'
*
      STRING_100(34:39) = X_UNIT(2:7)
      STRING_101(34:39) = X_UNIT(2:7)
      STRING_110(34:39) = X_UNIT(2:7)
      STRING_111(34:39) = X_UNIT(2:7)
*
*
* Option /CHANNEL
*
      FIRST_CHAN = 1
      LAST_CHAN  = NCHAN_FIT
      PARAM_OPT = .FALSE.
      IF (.NOT.CNV_MAP.AND. .NOT.INPUT_OPT) THEN
         DO WHILE (RECEIVER_FLAG(FIRST_CHAN).AND.
     $      FIRST_CHAN.LE.NCHAN)
            FIRST_CHAN = FIRST_CHAN+1
         ENDDO
         IF (FIRST_CHAN.GT.NCHAN) THEN
            CALL MESSAGE (2,1,'FIT','No channel available')
            GOTO 99
         ENDIF
      ENDIF
      IF (SIC_PRESENT(1,0)) THEN
         LAST_CHAN = 1                   ! chan 1 is default for option /CHAN
         IF (SIC_PRESENT(1,1)) THEN
            CALL SIC_I4 (LINE,1,1,ICHAN,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            IF (ICHAN.GT.NCHAN_FIT.OR.ICHAN.LT.1) THEN
               CALL MESSAGE (2,1,'FIT','Bad channel number ')
               GOTO 99
            ENDIF
            FIRST_CHAN = ICHAN
            LAST_CHAN  = ICHAN
         ENDIF
      ENDIF
*
* Option /FUNCTION
*
      GAUSS_2D = .FALSE.
      IFUNC    = 1
      IF (SIC_PRESENT(3,1)) THEN
         CALL SIC_I4 (LINE,3,1,IFUNC,.TRUE.,ERROR)
         IF (IFUNC.EQ.2) GAUSS_2D = .TRUE.
         IF (IFUNC.EQ.2 .AND. DOUBLEBEAM_FIT)
     $   CALL MESSAGE(2,1,'FIT',
     $   'Option /FUNCT 2 not allowed with doublebeam files')
* Initialize starting parameters (if they are present)
         IF (SIC_PRESENT(3,2)) THEN
            PARAM_OPT = .TRUE.
            CALL SIC_R8 (LINE,3,2,XPARAM0(1),.TRUE.,ERROR)
            CALL SIC_R8 (LINE,3,3,XPARAM0(2),.TRUE.,ERROR)
            CALL SIC_R8 (LINE,3,4,AMPLITUDE, .TRUE.,ERROR)
            CALL SIC_R8 (LINE,3,5,XPARAM0(4),.TRUE.,ERROR)
            XPARAM0(3) = AMPLITUDE
         ENDIF
      ENDIF
*
* Option /FILE
*
      NLUN   = 1
      LUN(1) = 6                         ! to write results on screen
      OUTPUT_OPT = .FALSE.
      IF (SIC_PRESENT(4,0)) THEN
         IF (SIC_PRESENT(4,1)) THEN
            CALL SIC_CH (LINE,4,1,FILENAME,NC,.TRUE.,ERROR)
            IF (ERROR) GOTO 99
            CALL SIC_LOWER(FILENAME)
            I = INDEX(FILENAME,'.')
            TYPE = FILENAME(I+1:I+3)
            IF (TYPE(1:3).EQ.'mdb'.OR.TYPE(1:3).EQ.'ekh'.OR.
     $      TYPE(1:3).EQ.'saa'.OR.TYPE(1:3).EQ.'cnv') THEN
               CALL MESSAGE (2,1,'FIT','Reserved outputfile extensions')
               GOTO 99
            ENDIF
            SNAM = FILENAME
         ENDIF
*
*         Open file for results
*
         IER = SIC_GETLUN(LUN(2))
         CALL SIC_PARSEF(SNAM,FILENAME,' ','.fit')
         NC = LENC(FILENAME)
         IF (SIC_PRESENT(4,2)) THEN
            IER = SIC_OPEN(LUN(2),FILENAME(1:NC),'NEW',.FALSE.)
         ELSE
            IER = SIC_OPEN(LUN(2),FILENAME(1:NC),'OLD',.FALSE.)
            IF (IER.NE.0) THEN
               IER = SIC_OPEN(LUN(2),FILENAME(1:NC),'NEW',.FALSE.)
            ELSE
600            READ (LUN(2),'(A)',END=700,ERR=700) INPUT_LINE
               GOTO 600
            ENDIF
         ENDIF
700      OUTPUT_OPT = .TRUE.
         NLUN = 2
*         ELSE
*            CALL MESSAGE (2,1,'FIT','Missing FILE name')
*            GOTO 99
      ENDIF
*      ENDIF
*
*     Option /LEVEL
*
      LEVEL_OPT = .FALSE.
      IF (SIC_PRESENT(5,0)) THEN
         IF (SIC_PRESENT(5,1)) THEN
            LEVEL_OPT = .TRUE.
            CALL SIC_R4 (LINE,5,1,LEVEL,.TRUE.,ERROR)
            IF (LEVEL.LT.0 .OR. LEVEL.GT.100) THEN
               CALL MESSAGE (2,1,'FIT','Invalid blanking level')
               GOTO 99
            ENDIF
            LEVEL = LEVEL/100.
         ELSE
            CALL MESSAGE (2,1,'FIT','Missing parameter')
            GOTO 99
         ENDIF
      ENDIF
*
* Allocate memory
*
*
* Get some memory and compute pointers
      MPAR = 8                           ! max number of parameters (must be even)
      NDATA = NX*NY
!      if (mod(ndata,2).ne.0) ndata=ndata+1   !even for 64 bits machines
      LWA = MPAR*(NDATA+5)+NDATA
      LWORK = MPAR + 2*LWA + 3*NDATA     ! mpar*i4 + lwa*r8+ ndata*r8 +ndata*r4
      IER = SIC_GETVM(LWORK , IADDR)
      IF (IER.NE.1) THEN
         CALL MESSAGE (2,1,'FIT','Memory allocation failure')
         GOTO 99
      ENDIF
      IP_IW  = POINTER(IADDR, MEMORY)
      IP_IWA = IP_IW  + MPAR
      IPFVEC = IP_IWA + 2*LWA
      IPDATA = IPFVEC + 2*NDATA
*
* Remove old FIT variables, create new variables
*
      CALL REMOVE_FIT_VARIABLE
      CALL CREATE_FIT_VARIABLE
*
* Loop on channels
*
*
*   Loop
      DO ICHAN = FIRST_CHAN, LAST_CHAN
         IF (SIC_CTRLC()) GOTO 98
         IF (.NOT.RECEIVER_FLAG(ICHAN).OR.INPUT_OPT.OR.CNV_MAP) THEN
            CALL SCALE_DATA(MEMORY(IPMAP+(ICHAN-1)*NX*NY),
     $      MEMORY(IPDATA),NX*NY,AMPLITUDE,IMAX,IMIN,DOUBLEBEAM_FIT,
     $      CNV_MAP)
            IF (NDATA.NE.0 .AND. ABS(AMPLITUDE).GT.1.D-32) THEN
               IF (PARAM_OPT) THEN
                  WRITE(6,*)(XPARAM0(I),I=1,4)
                  XPARAM(1) = XPARAM0(1)
                  XPARAM(2) = XPARAM0(2)
                  XPARAM(3) = XPARAM0(3)/AMPLITUDE
                  XPARAM(4) = XPARAM0(4)
               ELSE
                  I = IMAX - ((IMAX-1)/NX)*NX
                  J = IMAX/NX + 1
                  XPARAM(1) = VAL1 + INC1*(I-REF1)
                  XPARAM(2) = VAL2 + INC2*(J-REF2)
                  IF (ICHAN.EQ.FIRST_CHAN) THEN
                     XPARAM(3) = 1
                     XPARAM(4) = WIDTH_DEF
                  ENDIF
               ENDIF
               WIDTH = WIDTH_DEF
               IF (TELESCOPE(1:3).EQ.'HHT'.OR.
     $         TELESCOPE(1:3).EQ.'CSO') THEN
                  WIDTH=WIDTH_DEF*30./TELESCOPE_DIAMETER*250./FREQUENCY
               ENDIF
*
* One beam
* --------
               IF (.NOT.DOUBLEBEAM_FIT) THEN
*     Search 2 parameters : coordinates of the gaussian: x,y
                  MPAR =2
                  IOPT = 2
*                 tol = sqrt(r1mach(4))
                  TOL = 1.D-5
                  NPRINT = 0
                  CALL DNLS1E (FCN, IOPT, NDATA, MPAR, XPARAM,
     $            MEMORY(IPFVEC), TOL, NPRINT, INFO,
     $            MEMORY(IP_IW), MEMORY(IP_IWA), LWA)
                  IF (INFO.EQ.0.OR.INFO.GE.4) THEN
                     CALL PRINT_SLATEC_ERROR(INFO)
                  ENDIF
                  FSUMSQ  = DENORM(NDATA, MEMORY(IPFVEC))
                  CALL SAVE_FIT_RESULT_SLATEC(FSUMSQ,INFO,
     $            XPARAM(1),XPARAM(2),1.D0,WIDTH,WIDTH,0.D0,0.D0)
*     Search 4 parameters : x,y,intensity ,width
                  IF (PARAM_OPT) THEN
                     XPARAM(1) = XPARAM0(1)
                     XPARAM(2) = XPARAM0(2)
                  ENDIF
                  XPARAM(4) = WIDTH_DEF
                  MPAR = 4
                  CALL DNLS1E (FCN, IOPT, NDATA, MPAR, XPARAM,
     $            MEMORY(IPFVEC), TOL, NPRINT, INFO,
     $            MEMORY(IP_IW), MEMORY(IP_IWA), LWA)
                  FSUMSQ  = DENORM(NDATA, MEMORY(IPFVEC))
                  IF (INFO.EQ.0.OR.INFO.GE.4) THEN
                     CALL PRINT_SLATEC_ERROR(INFO)
                     FSUMSQ = FSUMSQ_MIN+1     ! force FSUMSQ>FSUMSQ_MIN
                  ENDIF
                  IF (FSUMSQ.LT.FSUMSQ_MIN)
     $            CALL SAVE_FIT_RESULT_SLATEC(FSUMSQ,INFO,
     $            XPARAM(1),XPARAM(2),XPARAM(3),XPARAM(4),
     $            XPARAM(4),0.D0,0.D0)
                  IF (.NOT.GAUSS_2D) THEN
                     MEAN_ELE = 0
                     IF (.NOT.INPUT_OPT)
     $               CALL GET_ELEVATION (XPARAM_MIN(1), XPARAM_MIN(2),
     $               MEAN_ELE,MEMORY(IP_SCAN_COORD),
     $               MEMORY(IP_SCAN_COORD+NRECORD),MEMORY(IP_EL))
                     DO I=1,NLUN
                        IF (ICHAN.EQ.FIRST_CHAN.AND..NOT.CNV_MAP) THEN
                           IF (.NOT.INPUT_OPT) THEN
                              WRITE(LUN(I),112) SCAN_NUMBER,IDAY,
     $                        MONTH(IMONTH),IYEAR,UT,LST
                           ELSE
                              WRITE(LUN(I),116)X_FILE(1:NC)
                           ENDIF
                           WRITE(LUN(I),'(A)') STRING_100
                        ENDIF
                        IF (ICHAN.EQ.FIRST_CHAN.AND. CNV_MAP) THEN
                           IF (.NOT.INPUT_OPT) THEN
                              WRITE(LUN(I),112) SCAN_NUMBER,IDAY,
     $                        MONTH(IMONTH),IYEAR,UT,LST
                           ELSE
                              WRITE(LUN(I),116)X_FILE(1:NC)
                           ENDIF
                           WRITE(LUN(I),'(A)') STRING_110
                        ENDIF
                        IF (X_RMAX.LT.10000.) THEN
                           WRITE(LUN(I),1)
     $                     ICHAN,XPARAM_MIN(1),XPARAM_MIN(2),
     $                     XPARAM_MIN(3)*AMPLITUDE,
     $                     XPARAM_MIN(4),MEAN_ELE,
     $                     FSUMSQ_MIN*AMPLITUDE/SQRT(NDATA*1.D0),
     $                     INFO_SOL
1                          FORMAT(I4,2X,2F11.2,2X,F8.2,4X,F9.2,
     $                     3X,F9.2,F10.0,I3)
                        ELSE
                           WRITE(LUN(I),2)
     $                     ICHAN,XPARAM_MIN(1),XPARAM_MIN(2),
     $                     XPARAM_MIN(3)*AMPLITUDE,
     $                     XPARAM_MIN(4),MEAN_ELE,
     $                     FSUMSQ_MIN*AMPLITUDE/SQRT(NDATA*1.D0),
     $                     INFO_SOL
2                          FORMAT(I4,2X,2F11.2,F10.0,4X,F9.2,
     $                     3X,F9.2,F10.0,I3)
                        ENDIF
                     ENDDO
                  ENDIF
*
*     Search 6 parameters (if gauss_2d)
                  IF (GAUSS_2D) THEN
                     IF (PARAM_OPT) THEN
                        XPARAM(1) = XPARAM0(1)
                        XPARAM(2) = XPARAM0(2)
                     ENDIF
                     XPARAM(4) = 1./XPARAM(4)
                     XPARAM(5) = XPARAM(4)*1.1
                     DO I=1,5
                        TEMP(I) = XPARAM(I)    ! save initial value
                     ENDDO
                     MPAR  = 6
                     ANGLE = 0
                     TOL=1.D-3
                     DO WHILE (ANGLE.LT.90)
                        XPARAM(6) = ANGLE*DEG_TO_RAD
                        DO I=1,5
                           XPARAM(I) = TEMP(I)
                        ENDDO
                        CALL DNLS1E (FCN, IOPT, NDATA, MPAR, XPARAM,
     $                  MEMORY(IPFVEC), TOL, NPRINT, INFO,
     $                  MEMORY(IP_IW), MEMORY(IP_IWA), LWA)
                        IF (INFO.EQ.0.OR.INFO.GE.4) THEN
                           CALL PRINT_SLATEC_ERROR(INFO)
                           FSUMSQ = FSUMSQ_MIN+1     ! force FSUMSQ>FSUMSQ_MIN
                        ENDIF
                        FSUMSQ  = DENORM(NDATA, MEMORY(IPFVEC))
                        IF (FSUMSQ.LT.FSUMSQ_MIN.OR.INFO.LT.4)
     $                  CALL SAVE_FIT_RESULT_SLATEC(FSUMSQ,INFO,
     $                  XPARAM(1),XPARAM(2),XPARAM(3),
     $                  XPARAM(4),XPARAM(5),XPARAM(6),0.D0)
                        ANGLE = ANGLE + 20.    ! step = 20 degrees
                     ENDDO
                     XPARAM_MIN(4) = 1./XPARAM_MIN(4)
                     XPARAM_MIN(5) = 1./XPARAM_MIN(5)
                     XPARAM_MIN(6) = XPARAM_MIN(6)/DEG_TO_RAD
                     IF (XPARAM_MIN(5).GE.XPARAM_MIN(4)) THEN
                        TEMP(1) = XPARAM_MIN(5)
                        XPARAM_MIN(5) = XPARAM_MIN(4)
                        XPARAM_MIN(4) = TEMP(1)
                        XPARAM_MIN(6) = 90  + XPARAM_MIN(6)
                     ENDIF
                     DO I=1,NLUN
                        IF (ICHAN.EQ.FIRST_CHAN.AND..NOT.CNV_MAP) THEN
                           IF (.NOT.INPUT_OPT) THEN
                              WRITE(LUN(I),112) SCAN_NUMBER,IDAY,
     $                        MONTH(IMONTH),IYEAR,UT,LST
                           ELSE
                              WRITE(LUN(I),116)X_FILE(1:NC)
                           ENDIF
                           WRITE(LUN(I),'(A)') STRING_101
                        ENDIF
                        IF (ICHAN.EQ.FIRST_CHAN.AND.CNV_MAP) THEN
                           IF (.NOT.INPUT_OPT) THEN
                              WRITE(LUN(I),112) SCAN_NUMBER,IDAY,
     $                        MONTH(IMONTH),IYEAR,UT,LST
                           ELSE
                              WRITE(LUN(I),116)X_FILE(1:NC)
                           ENDIF
                           WRITE(LUN(I),'(A)') STRING_111
                        ENDIF
                        XPARAM_MIN(6) = MOD(XPARAM_MIN(6),180.D0)
                        IF (XPARAM_MIN(6).LT.0)
     $                  XPARAM_MIN(6) = XPARAM_MIN(6)+180
                        IF (XPARAM_MIN(6).GT.90)
     $                  XPARAM_MIN(6) = XPARAM_MIN(6)-180
                        WRITE(LUN(I),
     $                  '(i4,2x,2f11.2,f10.0,2f9.2,f9.2,f9.0,i3)')
     $                  ICHAN,XPARAM_MIN(1),XPARAM_MIN(2),
     $                  XPARAM_MIN(3)*AMPLITUDE,
     $                  XPARAM_MIN(4),XPARAM_MIN(5),
     $                  XPARAM_MIN(6),
     $                  FSUMSQ_MIN*AMPLITUDE/SQRT(NDATA*1.D0),
     $                  INFO_SOL
                     ENDDO
                  ENDIF
                  XPARAM_MIN(3)=XPARAM_MIN(3)*AMPLITUDE
               ELSE
*
*     Double beam
*
*     Find position of gaussians
                  I = IMIN - ((IMIN-1)/NX)*NX
                  J = IMIN/NX + 1
                  XPARAM(3) = VAL1 + INC1*(I-REF1)
                  XPARAM(4) = VAL2 + INC2*(J-REF2)
* 4 parameters x_pos, y_pos, x_neg, y_neg
                  MPAR = 4
                  IOPT = 2
*                 tol = sqrt(r1mach(4))
                  TOL = 1.D-5
                  NPRINT = 0
                  CALL DNLS1E (FCN, IOPT, NDATA, MPAR, XPARAM,
     $            MEMORY(IPFVEC), TOL, NPRINT, INFO,
     $            MEMORY(IP_IW), MEMORY(IP_IWA), LWA)
                  IF (INFO.EQ.0.OR.INFO.GE.4) THEN
                     CALL PRINT_SLATEC_ERROR(INFO)
                  ENDIF
                  FSUMSQ  = DENORM(NDATA, MEMORY(IPFVEC))
                  CALL SAVE_FIT_RESULT_SLATEC(FSUMSQ,INFO,XPARAM(1),
     $            XPARAM(2),XPARAM(3),XPARAM(4),1.D0,-1.D0,WIDTH)
* 5 parameters x_pos, y_pos, x_neg, y_neg, width
                  MPAR = 5
                  XPARAM(5) = WIDTH
                  CALL DNLS1E (FCN, IOPT, NDATA, MPAR, XPARAM,
     $            MEMORY(IPFVEC), TOL, NPRINT, INFO,
     $            MEMORY(IP_IW), MEMORY(IP_IWA), LWA)
                  IF (INFO.EQ.0.OR.INFO.GE.4) THEN
                     CALL PRINT_SLATEC_ERROR(INFO)
                     FSUMSQ = FSUMSQ_MIN+1     ! force FSUMSQ>FSUMSQ_MIN
                  ENDIF
                  FSUMSQ  = DENORM(NDATA, MEMORY(IPFVEC))
                  IF (FSUMSQ.LT.FSUMSQ_MIN)
     $            CALL SAVE_FIT_RESULT_SLATEC(FSUMSQ,INFO,XPARAM(1),
     $            XPARAM(2),XPARAM(3),XPARAM(4),1.D0,-1.D0,
     $            XPARAM(5))
* 6 parameters x_pos, y_pos, x_neg, y_neg, a_pos, a_neg
                  MPAR = 6
                  WIDTH = XPARAM(5)
                  XPARAM(5) = 1
                  XPARAM(6) = -1
                  CALL DNLS1E (FCN, IOPT, NDATA, MPAR, XPARAM,
     $            MEMORY(IPFVEC), TOL, NPRINT, INFO,
     $            MEMORY(IP_IW), MEMORY(IP_IWA), LWA)
                  IF (INFO.EQ.0.OR.INFO.GE.4) THEN
                     CALL PRINT_SLATEC_ERROR(INFO)
                     FSUMSQ = FSUMSQ_MIN+1     ! force FSUMSQ>FSUMSQ_MIN
                  ENDIF
                  FSUMSQ  = DENORM(NDATA, MEMORY(IPFVEC))
                  IF (FSUMSQ.LT.FSUMSQ_MIN)
     $            CALL SAVE_FIT_RESULT_SLATEC(FSUMSQ,INFO,XPARAM(1),
     $            XPARAM(2),XPARAM(3),XPARAM(4),XPARAM(5),
     $            XPARAM(6),WIDTH)
                  DO I=1,NLUN
                     IF (ICHAN.EQ.FIRST_CHAN) THEN
                        IF (SCAN_TYPE.NE.4) THEN
                           IF (.NOT.INPUT_OPT) THEN
                              WRITE(LUN(I),112) SCAN_NUMBER,IDAY,
     $                        MONTH(IMONTH),IYEAR,UT,LST
                           ELSE
                              WRITE(LUN(I),116) X_FILE(1:NC)
                           ENDIF
                           WRITE(LUN(I),200)
                        ELSE
                           WRITE(LUN(I),201)
                        ENDIF
                     ENDIF
                     WRITE(LUN(I),
     $               '(i3,4f8.2,2f10.0,f9.2,f8.0,i3)')
     $               ICHAN,XPARAM_MIN(1),XPARAM_MIN(2),XPARAM_MIN(3),
     $               XPARAM_MIN(4), XPARAM_MIN(5)*AMPLITUDE,
     $               XPARAM_MIN(6)*AMPLITUDE,XPARAM_MIN(7),
     $               FSUMSQ_MIN*AMPLITUDE/SQRT(NDATA*1.D0),
     $               INFO_SOL
                  ENDDO
*
                  XPARAM_MIN(5)=XPARAM_MIN(5)*AMPLITUDE
                  XPARAM_MIN(6)=XPARAM_MIN(6)*AMPLITUDE
               ENDIF
               CALL UPDATE_FIT_VARIABLE (XPARAM_MIN,ICHAN)
            ELSE
               WRITE(6,'(a,i3,a)')
     $         '! ----- No data available for channel: ',
     $         ICHAN,' -----'
            ENDIF
         ENDIF
      ENDDO
      IF (DOUBLEBEAM_FIT.AND.SCAN_TYPE.NE.4)
     $CALL UPDATE_BEAMSEP(FIRST_CHAN,LAST_CHAN)
      IF (SECTION(7).EQ.1) CALL FIT_VARIABLES(.FALSE.,.FALSE.)
      IF (SECTION(7).EQ.2) CALL FIT_VARIABLES(.FALSE.,.TRUE.)
      IF (OUTPUT_OPT.AND.FIRST_CHAN.EQ.1.AND.
     $LAST_CHAN.EQ.NCHAN) THEN
         IP_AZ_POS  = POINTER (IADDR_FIT,MEMORY)
         IF (DOUBLEBEAM_FIT) THEN
            IP_AMP_POS = IP_AZ_POS  + NCHAN*5
            IP_AMP_NEG = IP_AMP_POS + NCHAN
            IP_AZ_OFF  = IP_AMP_NEG + NCHAN
            IP_EL_OFF  = IP_AZ_OFF  + NCHAN
            IP_WIDTH   = IP_AZ_POS  + NCHAN*4
         ELSE
            IP_AMP_POS = IP_AZ_POS  + NCHAN*2
            IP_AZ_OFF  = IP_AMP_POS + NCHAN*4
            IP_EL_OFF  = IP_AZ_OFF  + NCHAN
            IP_WIDTH1  = IP_AMP_POS + NCHAN
            IP_WIDTH2  = IP_WIDTH1  + NCHAN
         ENDIF
         CALL R4TOR4(MEMORY(IP_AZ_OFF),RES,NCHAN)
         DO I=1,NCHAN
            WRITE (INPUT_LINE((I-1)*8+1:I*8),'(f8.2)') RES(I)
         ENDDO
         WRITE (LUN(2),113) INPUT_LINE(1:NCHAN*8)
         CALL R4TOR4(MEMORY(IP_EL_OFF),RES,NCHAN)
         DO I=1,NCHAN
            WRITE (INPUT_LINE((I-1)*8+1:I*8),'(f8.2)') RES(I)
         ENDDO
         WRITE (LUN(2),114) INPUT_LINE(1:NCHAN*8)
         CALL R4TOR4(MEMORY(IP_AMP_POS),RES,NCHAN)
         IF (DOUBLEBEAM_FIT) THEN
            CALL R4TOR4(MEMORY(IP_AMP_NEG),DUMMY,NCHAN)
            DO I=1,NCHAN
               RES(I) = (RES(I)+ABS(DUMMY(I)))
            ENDDO
            CALL R4TOR4(MEMORY(IP_WIDTH),DUMMY,NCHAN)
            DO I=1,NCHAN
               RES(I) = RES(I)*DUMMY(I)
            ENDDO
         ELSE
            CALL R4TOR4(MEMORY(IP_WIDTH1),DUMMY,NCHAN)
            DO I=1,NCHAN
               RES(I) = RES(I)*DUMMY(I)
            ENDDO
            CALL R4TOR4(MEMORY(IP_WIDTH2),DUMMY,NCHAN)
            DO I=1,NCHAN
               RES(I) = RES(I)*DUMMY(I)
            ENDDO
         ENDIF
         DUMMY(1) = RES(1)
         DO I=1,NCHAN
            IF (ABS(RES(I)).GT.1.D-05)
     $      WRITE (INPUT_LINE((I-1)*8+1:I*8),'(f8.2)')
     $      DUMMY(1)/RES(I)
         ENDDO
         WRITE (LUN(2),115) INPUT_LINE(1:NCHAN*8)
      ENDIF
*
98    CALL GDF_FRIS (X_ISLO,ERROR)
      CALL FREE_VM(LWORK,IADDR)
      IF (OUTPUT_OPT) THEN
         CLOSE(UNIT=LUN(2))
         IER = SIC_FRELUN(LUN(2))
      ENDIF
      RETURN
*
99    CALL GDF_FRIS (X_ISLO,ERROR)
      ERROR = .TRUE.
      RETURN
*
*
990   CALL MESSAGE (2,1,'FIT','FIT needs NAG library ')
      CALL GDF_FRIS (X_ISLO,ERROR)
      CALL FREE_VM(LWORK,IADDR)
      IF (OUTPUT_OPT) THEN
         CLOSE(UNIT=LUN(2))
         IER = SIC_FRELUN(LUN(2))
      ENDIF
      RETURN
*
* Format
112   FORMAT('! Scan: ',I4.4,'  Date: ',I2.2,'-',A3,'-',
     $I4,'  Time: ',F7.1,' (UT) / ',F7.1,' (LST)')
113   FORMAT('!',/'! Az: ',A)
114   FORMAT('!',/'! El: ',A)
115   FORMAT('!',/'! Ga: ',A)
116   FORMAT('! Input: ',A)
*
200   FORMAT('! Chan  X_pos   Y_pos   X_neg   Y_neg ',
     $'  I_pos    I_neg     FWHM    Sigma  Info')
201   FORMAT('! Plane X_pos   Y_pos   X_neg   Y_neg ',
     $'  I_pos    I_neg     FWHM    Sigma  Info')
      END
*
      SUBROUTINE FCN(IFLAG,M,NVPAR,XPAR,F,FJAC,LJC)
C-----------------------------------------------------------------------
C FCN is called by DNLS1E
C
C     IFLAG   Integer    Code for print out                Input
C     M       Integer    Number of data points             Input
C     NVPAR   Integer    Number of variable parameters     Input
C     X       Real*8     Variable Parameters               Input
C     F       Real*8     Function value at each point      Output
C     FJAC    Real*8     Partial derivatives of F          Output
C     LJC     Integer    First dimension of FJAC           Input
c
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
C Global:
      INCLUDE 'nic.inc'
C Dummy:
      INTEGER M, NVPAR, LJC,IFLAG
      REAL*8 XPAR(NVPAR), F(M), FJAC(LJC,NVPAR)
C Local:
      INTEGER K, I, L, J,II
      REAL*8 XX,YY
      REAL*8 Y,DY(8)
      REAL*8  REF1,INC1,VAL1,REF2,VAL2,INC2,WIDTH
      LOGICAL TEST_BLANKING
*
      INTEGER*4 IPDATA
      INTEGER 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/
*
      IF (IFLAG.EQ.0) THEN
         WRITE(6,*) 'iflag=0...'
         RETURN
      ENDIF
*
* Compute F and its derivatives
*
      K = 1
      DO II=1,NX*NY
         IF (.NOT.TEST_BLANKING(MEMORY(IPDATA+II-1))) THEN
            I = II - ((II-1)/NX)*NX
            J = II/NX + 1
            XX = VAL1 + INC1*(I-REF1)
            YY = VAL2 + INC2*(J-REF2)
*    compute y and dy(nvpar)
            CALL COMP_GAUSS(XX,YY,XPAR,DOUBLEBEAM_FIT)
            IF (IFLAG.EQ.1) THEN
*              compute  f(k) = y - memory(ipmap+ii-1)
               CALL DIFF(Y,MEMORY(IPDATA+II-1),F(K))
            ELSEIF (IFLAG.EQ.2) THEN
               DO L=1, NVPAR
                  FJAC(K,L) = DY(L)
               ENDDO
            ENDIF
            K = K+1
         ENDIF
      ENDDO
      RETURN
      END
*
*-----------------------------------------------------------------------
      SUBROUTINE SAVE_FIT_RESULT_SLATEC(FSUMSQ,INFO,
     $XPAR1,XPAR2,XPAR3,XPAR4,XPAR5,XPAR6,XPAR7)
*-----------------------------------------------------------------------
*
*  Save result of a fit
*
*  Input:  fsumsq,xpar1,xpar2,xpar3,xpar4,xpar5,xpar6,xpar7
*  Output: fsumsq_min,xpar(6)
*-----------------------------------------------------------------------
      REAL*8   FSUMSQ
      REAL*8   XPAR1,XPAR2,XPAR3,XPAR4,XPAR5,XPAR6,XPAR7
      INTEGER  INFO
*
      REAL*8   XPARAM_MIN(8),FSUMSQ_MIN
      INTEGER*4 INFO_SOL
      COMMON/PARAM_MIN/XPARAM_MIN,FSUMSQ_MIN,INFO_SOL
*
      FSUMSQ_MIN = FSUMSQ
      INFO_SOL  = INFO
      XPARAM_MIN(1) = XPAR1
      XPARAM_MIN(2) = XPAR2
      XPARAM_MIN(3) = XPAR3
      XPARAM_MIN(4) = XPAR4
      XPARAM_MIN(5) = XPAR5
      XPARAM_MIN(6) = XPAR6
      XPARAM_MIN(7) = XPAR7
      RETURN
      END
*-----------------------------------------------------------------------
      SUBROUTINE PRINT_SLATEC_ERROR(INFO)
      INTEGER INFO
*
      IF (INFO.EQ.0) THEN
         PRINT *,'F-DNLS1E, improper input parameters, info= ',INFO
      ELSEIF (INFO.EQ.6 .OR. INFO.EQ.7) THEN
         PRINT *,'F-DNLS1E, TOL too small, info= ',INFO
      ELSEIF (INFO.EQ.5) THEN
         PRINT *,'F-DNLS1E, not converging, info= ',INFO
      ELSEIF (INFO.EQ.4) THEN
         PRINT *,'W-DNLS1E, FVEC orthog. to Jacob. col., info= ',
     $   INFO
      ENDIF
      END
