PROGRAM FIT

      IMPLICIT NONE

      DOUBLE PRECISION PAR,YMODEL,OMEGA,D2,ST1,EPS,MINSTEP,MIN1A,C,PI,
     &     YM
      INTEGER NPAR,I,J,NDATA,NFIT,WINDOW,RES,PLT,AMP,PER
      LOGICAL CONV
      CHARACTER PFILE*6,RFILE*6,DATAFILE*15,ANSWER*1,AFILE*6,TFILE*6
      PARAMETER (NPAR=6,PFILE='fitplt',RFILE='fitout',AFILE='fitamp',
     &     TFILE='fitper',PI=3.14159265358979323846,RES=11,PLT=12,
     &     AMP=13,PER=14)
      DIMENSION PAR(NPAR),C(3)

      COMMON /SUNE/T(200000),Y(200000),DATA(2)
      DOUBLE PRECISION T,Y
      INTEGER DATA
      
      EXTERNAL CHISQ

      NDATA=200000
      EPS=1.0E-6
      MINSTEP=1.0E-16

      WRITE(*,*) 'Size of moving window: '
      READ(*,*) WINDOW
      CALL LAES(T,Y,NDATA,DATAFILE)
      
 5    OPEN(UNIT=RES,FILE=RFILE,STATUS='UNKNOWN')
      OPEN(UNIT=PLT,FILE=PFILE,STATUS='UNKNOWN')
      OPEN(UNIT=AMP,FILE=AFILE,STATUS='UNKNOWN')
      OPEN(UNIT=PER,FILE=TFILE,STATUS='UNKNOWN')

      WRITE(RES,'(3A)') '# Fit of f(t) = a*sin(omega*t + phi) + m',
     &     ' to the data in ',DATAFILE
      WRITE(RES,'(A1,A12,8(XXA))') '#',' time       ','amplitude  ',
     &     'omega      ','phi        ','empty      ','empty      ',
     &     'mean value ','period     ','amplitude^2'
      WRITE(PLT,'(3A)') '# Fit of f(t) = a*sin(omega*t + phi) + m',
     &     ' to the data in ',DATAFILE
      WRITE(PLT,'(A1,A12,2(XXA11))') '#',' time       ',
     & 'data       ','fit        '

      NFIT=(NDATA-MOD(NDATA,WINDOW))/WINDOW
      WRITE(*,*) 'Now doing',NFIT,' fits.'
      DO 20,I=0,NFIT
         DATA(1)=I*WINDOW
         DATA(2)=DATA(1)+WINDOW

         OMEGA=0.18
         ST1=0.001
         CALL BISEC(CHISQ,OMEGA,MIN1A,ST1,EPS,MINSTEP,D2)
         CALL GETC(OMEGA,C)
         
         IF (C(1).EQ.0.0) THEN
            PAR(1)=C(2)
            PAR(3)=ASIN(C(2)/PAR(1))
         ELSE
            PAR(3)=ATAN2(C(2),C(1))
            PAR(1)=C(2)/SIN(PAR(3))
         END IF
         PAR(2)=OMEGA
         PAR(6)=C(3)
*         WRITE(*,'(A,F9.2,A,F9.2,A,F8.4,A,F8.4)') '[',T(DATA(1)),':',
*     &        T(DATA(2)),']  T = ',2.0*PI/OMEGA,' amplitude = ',
*     &        PAR(1)*1000.0
         WRITE(RES,'(9(XG12.5))') T(DATA(1)+WINDOW/2),PAR,2.0*PI/OMEGA
     &        ,PAR(1)**2.0
         WRITE(AMP,*) PAR(1)
         WRITE(PER,*) 2.0*PI/OMEGA
         DO 10 J=DATA(1),DATA(2)
            YM=YMODEL(T(J),PAR)
            WRITE(PLT,'(4(XG14.5))') T(J),Y(J),YM,ABS(YM-Y(J))
*            WRITE(PLT,*) T(J),Y(J),YM,(YM-Y(J))*10.0+PAR(6),
*     &           ABS((YM-Y(J))*10.0)+PAR(6)
 10      CONTINUE
*         CALL GNU(T(DATA(1)),T(DATA(2)))
 20   CONTINUE
      CLOSE(RES)
      CLOSE(PLT)
      CLOSE(AMP)
      CLOSE(PER)
      WRITE(*,'(3A)') 'Done. You can now [q]uit or [c]ontinue',
     &     ' to work with ',DATAFILE
      READ(*,*) ANSWER
      IF ((ANSWER.EQ.'c').OR.(ANSWER.EQ.'C')) THEN
         WRITE(*,*) 'Size of moving window: '
         READ(*,*) WINDOW
         GOTO 5
      END IF
      END

*************************************************************************
      
      DOUBLE PRECISION FUNCTION YMODEL(T,P)

      IMPLICIT NONE

      DOUBLE PRECISION T,P
      DIMENSION P(6)
      YMODEL=P(1)*SIN(P(2)*T+P(3))+P(6)
*      YMODEL=P(1)*SIN(P(2)*T+P(3))+P(4)*SIN(2.0*P(2)*T+P(5)) + P(6)
      RETURN
      END

*************************************************************************

      SUBROUTINE GNU(T1,T2)

      IMPLICIT NONE

      DOUBLE PRECISION T1,T2

      INTEGER IUNIT
      PARAMETER (IUNIT=10)
      OPEN(UNIT=IUNIT,FILE='gnufile',STATUS='UNKNOWN',FORM='FORMATTED')
      WRITE(IUNIT,'(A)') 'set term X11'
      WRITE(IUNIT,'(A,F8.1,A,F8.1,A)') 'set xrange[',T1,':',T2,']'
      WRITE(IUNIT,'(2A)')  "plot 'fitplt' title 'data' w l,",
     &     "'fitplt' u 1:3 title 'fit' w l" 
      WRITE(IUNIT,'(A)')  'pause -1 "Press <return> to continue."'
      CLOSE(IUNIT)
      CALL SYSTEM("gnuplot gnufile")
      OPEN(UNIT=IUNIT,FILE='gnufile',STATUS='OLD',FORM='FORMATTED')
      CLOSE(IUNIT,STATUS='DELETE')
      RETURN
      END