      SUBROUTINE FIT(AT,AF,AA,SIGAA,AJDH,IPART,TITLE)
      COMMON /PLOTFL/PLNAME
      double precision work1(3,3)
      DIMENSION AT(1),AF(2,256),AA(1),AMP(256),ALINE(60),BIN(10,50)     
      DIMENSION COV(3,3),COVINV(3,3),JUNK(3),V(2),Y(3),DERIV(3)   
      DIMENSION SIGAA(1),SIGBIN(50),AFIT(53),PH(53),iwork2(3)
      CHARACTER*1 ALINE,TITLE*80,PLNAME*12
      DOUBLE PRECISION AJDH   
C     
C     ROUTINE TO LEAST SQUARES FIT THE AMPLITUDES WITH A SINUSOID 
C     
      ISQ=0 
      TWOPI=2.0*3.141592653   
      READ 10,STARTP,ENDP,DP,AMAX   
   10 FORMAT()    
      IF(DP.LT.0.0)THEN 
        DP=-DP    
        ISQ=1     
        PRINT 11  
   11   FORMAT(12X,'FIT TO THE SQUARES OF THE AMPLITUDES')  
      END IF
      FMEAN=0.0   
      IF(ENDP.LT.0.0)THEN     
        ENDP=-ENDP
        READ *,FMEAN    
      END IF
      PRINT 20,AJDH,STARTP,ENDP,DP  
   20 FORMAT(12X,22HTIME STARTS AT HJD OF ,F18.7/     
     *       12X,23HSTARTING TEST PERIOD = ,F9.1,5H SECS/   
     *       12X,23HENDING TEST PERIOD   = ,F9.1,5H SECS/   
     *       12X,23HPERIOD INCREMENT     = ,F9.1,5H SECS/)  
      IF(FMEAN.NE.0.0)PRINT 25,FMEAN
   25 FORMAT(12X,'TRYING TO FORCE MEAN = ',F9.1/)     
      NP=(ENDP-STARTP)/DP+1   
      SUMMIN=1.E30
      NUM=3 
      DO 30 J=1,IPART   
        IF(ISQ.EQ.0)AMP(J)=AA(J)-FMEAN    
        IF(ISQ.EQ.1)AMP(J)=AA(J)**2-FMEAN 
        IF(J.EQ.1)PRINT 28    
   28   FORMAT(12X,'DATA TO BE FITTED ARE : '/) 
        PRINT 29,AT(J),AMP(J) 
   29   FORMAT(12X,'TIME = ',F10.1,'     A = ',1PE15.5)     
   30 CONTINUE    
      DO 100 J=1,NP     
        AMEAN=0.0 
        A=0.0     
        B=0.0     
        DO 32 L=1,NUM   
          Y(L)=0.0
          DO 32 M=1,NUM 
            COV(L,M)=0.0
   32   CONTINUE  
        P=STARTP+(J-1)*DP     
        F=1.0/P   
        W=TWOPI*F 
        DO 35 K=1,IPART 
          T=AT(K) 
          ARG=AMOD(W*T,TWOPI) 
          S=SIN(ARG)    
          C=COS(ARG)    
          DERIV(1)=1    
          DERIV(2)=C    
          DERIV(3)=S    
          DO 34 L=1,NUM 
            DO 33 M=1,NUM     
              COV(L,M)=COV(L,M)+DERIV(L)*DERIV(M)     
              COVINV(L,M)=COV(L,M)  
   33       CONTINUE    
            Y(L)=Y(L)+DERIV(L)*AMP(K)     
   34     CONTINUE
   35   CONTINUE  
        V(1)=1.0  

c   Invert covinv; inverse appears in work1

        do 800 i=1,num
        do 801 jj=1,num
 801    work1(i,jj)=0.
 800    work1(i,i)=1.
        call ludcmp(covinv,num,3,iwork2,work2)
        do 802 i=1,num
        call lubksb(covinv,num,3,iwork2,work1(1,i))
 802    continue

c       Put the inverse matrix into COVINV
        do i = 1,num
           do jj = 1,num
              covinv(i,jj) = work1(i,jj)
           enddo
        enddo

   36   DO 38 M=1,NUM   
          AMEAN=AMEAN+COVINV(1,M)*Y(M)    
          A=A+COVINV(2,M)*Y(M)
          B=B+COVINV(3,M)*Y(M)
   38   CONTINUE  
        SUMSQ=0.0 
        DO 45 K=1,IPART 
          T=AT(K) 
          TERM=AMP(K)-B*SIN(W*T)-A*COS(W*T)-AMEAN     
          SUMSQ=SUMSQ+TERM*TERM     
   45   CONTINUE  
        IF(SUMSQ.GT.SUMMIN)GO TO 50 
        IND=J     
        SUMMIN=SUMSQ    
        AMNMIN=AMEAN    
        SIGM=SQRT(COVINV(1,1)*SUMSQ/(IPART-NUM))
        AMIN=A    
        SIGAMN=SQRT(COVINV(2,2)*SUMSQ/(IPART-NUM))    
        BMIN=B    
        SIGBMN=SQRT(COVINV(3,3)*SUMSQ/(IPART-NUM))    
        PMIN=P    
   50   CONTINUE  
        IF(J.EQ.1)PRINT 60    
   60   FORMAT(//15X,6HPERIOD,7X,'AMP',8X,'PHI',9X,'SUMSQ'/12X,46(1H=)) 
        R=SQRT(A**2+B**2)     
        PHI=ATAN2(-B,A) 
        PRINT 70,P,R,PHI,SUMSQ
   70   FORMAT(12X,F9.1,3X,1PE10.3,1X,E10.3,3X,E10.4) 
        IF(J.EQ.NP)PRINT 80   
   80   FORMAT(12X,46(1H=)/)  
  100 CONTINUE    
      AMEAN=AMNMIN+FMEAN
      R=SQRT(AMIN**2+BMIN**2) 
      PHI=ATAN2(-BMIN,AMIN)   
      SIGR=SQRT((AMIN*SIGAMN)**2+(BMIN*SIGBMN)**2)/R  
      SIGPHI=SQRT((AMIN*SIGBMN)**2+(BMIN*SIGAMN)**2)/(R**2) 
      PRINT 105,PMIN,AMIN,BMIN,R,SIGR,PHI,SIGPHI,AMEAN,SIGM 
  105 FORMAT(12X,'BEST FIT PARAMETERS ARE :'/   
     *       14X,'BEST PERIOD = ',T40,F9.1,' SECS'/   
     *       14X,'BEST A = ',T40,F9.1/14X,'BEST B = ',T40,F9.1/   
     *       14X,'BEST AMPLITUDE = ',T40,F9.1/  
     *       14X,'WITH SIGMA = ',T40,F9.1/
     *       14X,'BEST PHI = ',T40,F9.3/14X,'WITH SIGMA = ',T40,F9.3/   
     *       14X,'FITTED MEAN = ',T40,F9.1/14X,'WITH SIGMA = ',T40,F9.1)
      PAV=0.0     
      ISKIP=0     
      IF(AMAX.LT.0.0)THEN     
        AMAX=-AMAX
        ISKIP=1   
      END IF
      DO 110 J=1,IPART  
        PAV=PAV+AF(2,J)/IPART 
        IF(ISKIP.EQ.1)GO TO 110     
        IF(ISQ.EQ.0.AND.AA(J).GT.AMAX)AMAX=AA(J)
        IF(ISQ.EQ.1.AND.AA(J)**2.GT.AMAX)AMAX=AA(J)**2
        IF(AMAX.LT.R+AMEAN)AMAX=R+AMEAN   
  110 CONTINUE    
      NB=50 
      DO 111 J=1,NB     
        DO 111 K=1,10   
          BIN(K,J)=0.0  
  111 CONTINUE    
      DPHASE=1./NB
      DO 114 J=1,IPART  
        IPHASE=AT(J)/PMIN     
        PHASE=AT(J)/PMIN-IPHASE     
        NOBIN=PHASE/DPHASE+1  
        NOBIN=MIN0(MAX0(1,NOBIN),NB)
        DO 112 K=1,10   
          IF(BIN(K,NOBIN).GT.0.0)GO TO 112
          IF(ISQ.EQ.0)BIN(K,NOBIN)=AA(J)  
          IF(ISQ.EQ.1)BIN(K,NOBIN)=AA(J)**2     
          IF(ISQ.EQ.0.AND.K.EQ.1)SIGBIN(NOBIN)=SIGAA(J)     
          IF(ISQ.EQ.1.AND.K.EQ.1)SIGBIN(NOBIN)=SIGAA(J)*AA(J)     
          GO TO 114     
  112   CONTINUE  
  114 CONTINUE    
      STEPA=60.0/AMAX   
      PRINT 120,PAV,PMIN
  120 FORMAT(1H1,11X,107(1H-)/12X,1HI,5X,19HBEST FIT OF MAXIMUM,  
     *       20H AMPLITUDES FOR THE ,F5.1,17H S PERIOD WITH A ,   
     *       F9.1,11H S SINUSOID,19X,1HI/12X,107(1H-)/12X,1HI,    
     *       6X,5HPHASE,72X,3HAMP,10X,'SIGMA',4X,1HI/12X,107(1H-))
      DO 125 J=1,60     
  125 ALINE(J)=' '
      DO 150 J=1,NB     
        PHASE=FLOAT(J-1)/NB   
        PH(J)=PHASE     
        ARG=TWOPI*PHASE+PHI   
        AFIT(J)=(AMEAN+R*COS(ARG))  
        JJ=AFIT(J)*STEPA
        JJ=MIN0(MAX0(1,JJ),60)
        ALINE(JJ)='+'   
        DO 126 K=1,10   
          IF(BIN(K,J).EQ.0.0)GO TO 126    
          JJ=BIN(K,J)*STEPA   
          JJ=MIN0(MAX0(1,JJ),60)    
          ALINE(JJ)='*' 
  126   CONTINUE  
        IF(ALINE(JJ).EQ.'*')THEN    
          PRINT 130,PHASE,ALINE,BIN(1,J),SIGBIN(J)    
        ELSE
          PRINT 140,PHASE,ALINE     
        END IF    
  130   FORMAT(12X,'I',F11.3,1X,'I',60A1,1X,'I',F13.5,2X,F13.5,2X,'I')  
  140   FORMAT(12X,'I',F11.3,1X,'I',60A1,1X,'I',30X,'I')    
        DO 145 K=1,60   
  145   ALINE(K)=' '    
  150 CONTINUE    
      PRINT 170   
  170 FORMAT(12X,107(1H-)/)   
      XL=8.0     
      YL=5.0
C
C     INITIALISE THE PLOTTING
C
c     CALL SETCC
c     CALL PLFILE(PLNAME)
c     CALL PLOTS(IERR,IDUMMY,28)
c     IF(IERR.NE.0) STOP 'ERROR ON OPENING SPOOLTEK CALL'
c     CALL FACTOR(2.54) 
      QBQB=XL+2.0
      PBPB=YL+2.0
c     CALL FRAME(QBQB,PBPB)  
c     CALL PLOT(1.0,1.0,-3)   
      DPHDX=2.0/XL
      DADY=AMAX/YL
c     CALL AXIS(0.0,0.0,'PHASE',-5,XL,0.0,0.0,DPHDX)  
c     CALL AXIS(0.0,0.0,'AMPLITUDE',9,YL,90.0,0.0,DADY)     
c     CALL AXIS(0.0,YL,' ',1,XL,0.0,0.0,DPHDX)  
c     CALL AXIS(XL,0.0,' ',-1,YL,90.0,0.0,DADY) 
c     CALL SYMBOL(0.0,YL+0.6,0.2,TITLE,0.0,78)
      PH(51)=1.0  
      ARG=PHI     
      AFIT(51)=(AMEAN+R*COS(ARG))   
      IC=1  
      AFIT(53)=DADY     
      PH(53)=DPHDX
  180 CONTINUE
c     CALL LINE(PH,AFIT,51,1,0,0)   
      DO 200 J=1,IPART  
        XX=(AT(J)/PMIN-IFIX(AT(J)/PMIN))/DPHDX  
        IF(ISQ.EQ.0)THEN
          YY=AA(J)/DADY 
          DY=SIGAA(J)   
        ELSE
          YY=AA(J)**2/DADY    
          DY=SIGAA(J)*AA(J)/DADY    
        END IF    
c       CALL SYMBOL(XX,YY,0.08,11,0.0,-1) 
c       CALL PLOT(XX,AMIN1(YL,YY+DY),3)   
c       CALL PLOT(XX,AMAX1(0.0,YY-DY),2)  
  200 CONTINUE    
      IF(IC.EQ.1)THEN   
        IC=IC+1   
c       CALL PLOT(XL/2.0,0.0,-3)    
        GO TO 180 
      END IF
c     CALL PLOT(0.0,0.0,999)
      RETURN
      END   
