      SUBROUTINE LSQRV(Q,TQ,NI,NF,NFREQ,FREQ,AMP,PHI,SIGMAA,SIGMAP,A0,  
     *                 SIGA0,MD,YMEAN,YSTDEV,RMEAN,RSTDEV)  
C     
C     ROUTINE TO FIT A FUNCTION OF THE FORM     
C             A0 + AMP(1)*COS(2*PI*FREQ(1)*TQ(I) + PHI(1))  
C     TO THE DATA IN Q AND TQ BY LEAST SQUARES.  THE ROUTINE WAS  
C     MOTIVATED BY A NEED TO FIT RADIAL VELOCITY ORBITS     
C     
      double precision work1(3,3)
      dimension iwork2(3)
      DIMENSION Q(NF),FREQ(1),AMP(1),PHI(1),SIGMAA(1),SIGMAP(1)    
      DOUBLE PRECISION TQ(NF)
      DOUBLE PRECISION SQQ,SUMRES,SSQRES,ARG,TWOPI,TI,OMEGA 
      DOUBLE PRECISION V(2),JUNK(3),COV(3,3),COVINV(3,3),UNIT(3,3)
      DOUBLE PRECISION Y(3),CM(3,3),DERIV(3),ANSWER(3)
      TWOPI=2.0D0*3.141592653D0     
      N=NF-NI+1   
      OMEGA=TWOPI*FREQ(1)     
      IDEBUG=0    
      IF(MD.LT.0)IDEBUG=1     
      DO 20 I=1,3 
        Y(I)=0.0  
        ANSWER(I)=0.0   
        DO 20 J=1,3     
          COV(I,J)=0.0  
          CM(I,J)=0.0   
   20 CONTINUE    
C     
C     CALCULATE THE NORMAL EQUATION MATRIX AND Y VECTOR     
C     
      SQQ=0.0     
      DO 50 I=NI,NF     
        TI=TQ(I)  
        ARG=DMOD(OMEGA*TI,TWOPI)    
        DERIV(1)=DCOS(ARG)    
        DERIV(2)=DSIN(ARG)    
        DERIV(3)=1.0    
        DO 40 J=1,3     
          DO 35 K=1,3   
            COV(J,K)=COV(J,K)+DERIV(J)*DERIV(K) 
   35     CONTINUE
          Y(J)=Y(J)+Q(I)*DERIV(J)   
   40   CONTINUE  
        SQQ=SQQ+Q(I)*Q(I)     
   50 CONTINUE    
      YMEAN=Y(3)/N
      YSTDEV=DSQRT((SQQ-Y(3)*Y(3)/N)/(N-1))     
C     
C     INVERT THE NORMAL MATRIX 
C     
      DO 60 J=1,3 
      DO 60 K=1,3 
   60 COVINV(J,K)=COV(J,K)    
      V(1)=1.0    
      num = 3

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 j = 1,num
            covinv(i,j) = work1(i,j)
         enddo
      enddo
C     
C     FROM HERE DOWN TO ABOUT LABEL 135 ARE DEBUG CALCULATIONS & PRINTS 
C     
   65 IF(IDEBUG.EQ.1)THEN     
        PRINT 70  
   70   FORMAT(//' COVARIANCE MATRIX IS : '/)   
        DO 90 I=1,3     
          PRINT 80,(COV(I,J),J=1,3) 
   80     FORMAT(1X,1P10E13.4)
          PRINT 80
   90   CONTINUE  
        PRINT 95  
   95   FORMAT(///' INVERSE IS : '/)
        DO 100 I=1,3    
          PRINT 80,(COVINV(I,J),J=1,3)    
          PRINT 80
  100   CONTINUE  
        PRINT 105 
  105   FORMAT(///' THE Y VECTOR IS : '/) 
        PRINT 80,(Y(J),J=1,3) 
        PRINT 80  
C     
C       TEST FOR GOOD MATRIX INVERSION BY PRODUCING A UNIT MATRIX 
C       AND PRINTING IT AND THEN CALCULATE THE CONTRIBUTION TO THE
C       SOLUTIONS BY CALCULATING A MATRIX GIVEN BY COVINV(I,J)*Y(J)     
C     
        PRINT 120 
        DO 125 I=1,3    
          DO 115 K=1,3  
            UNIT(I,K)=0.0     
            DO 110 J=1,3
              UNIT(I,K)=UNIT(I,K)+COVINV(I,J)*COV(J,K)
  110       CONTINUE    
  115     CONTINUE
  120     FORMAT(///' PRODUCT OF COVARIANCE AND ITS INVERSE IS : '/)    
          PRINT 80,(UNIT(I,K),K=1,3)
          PRINT 80
  125   CONTINUE  
        PRINT 128 
  128   FORMAT(///' PRODUCT OF COVARIANCE INVERSE AND Y VECTOR :'/)     
        DO 135 I=1,3    
          DO 130 J=1,3  
            UNIT(I,J)=COVINV(I,J)*Y(J)    
  130     CONTINUE
          PRINT 80,(UNIT(I,J),J=1,3)
          PRINT 80
  135   CONTINUE  
      END IF
C     
C     CALCULATE THE ANSWERS   
C     
      DO 140 I=1,3
        DO 140 J=1,3    
          ANSWER(I)=ANSWER(I)+COVINV(I,J)*Y(J)  
  140 CONTINUE    
      A0=ANSWER(3)
      AMP(1)=DSQRT(ANSWER(1)*ANSWER(1)+ANSWER(2)*ANSWER(2)) 
      PHI(1)=DATAN2(-ANSWER(2),ANSWER(1)) 
C     
C     CALCULATE SIGMAE, THE VARIANCE OF THE RESIDUALS AND REPLACE 
C     Q BY THE RESIDUALS IF MD = 1  
C     
      SUMRES=0.0  
      SSQRES=0.0  
      DO 150 I=NI,NF    
        RESID=Q(I)
        TI=TQ(I)  
        ARG=DMOD(OMEGA*TI+PHI(1),TWOPI)   
        RESID=RESID-A0-AMP(1)*DCOS(ARG)   
        SUMRES=SUMRES+RESID   
        SSQRES=SSQRES+RESID**2
        IF(MD.EQ.1)Q(I)=RESID 
  150 CONTINUE    
      SIGMAE=(SSQRES-(SUMRES**2/N))/(N-1) 
      RMEAN=SUMRES/N    
      RSTDEV=SQRT(SIGMAE)     
      SIG1=DSQRT(SIGMAE*COVINV(1,1))
      SIG2=DSQRT(SIGMAE*COVINV(2,2))
      SIGA0=DSQRT(SIGMAE*COVINV(3,3))     
      SIGMAA(1)=DSQRT((ANSWER(1)*SIG1)**2+(ANSWER(2)*SIG2)**2)/AMP(1)   
      SIGMAP(1)=DSQRT((ANSWER(1)*SIG2)**2+(ANSWER(2)*SIG1)**2)/AMP(1)**2
      RETURN
      END   
