module linear_least_squares
       use least_squares

       ! functions for fitting y=a(2)*x+a(1)
       interface lsq_linear
          module procedure mp_lsq_linear_sd
          module procedure mp_lsq_linear_nosd
          module procedure mp_lsq_linear_noa
       end interface

       ! functions for fitting y=a(1)*x
       interface lsq_linear_no_b
          module procedure mp_lsq_linear_no_b_sd
          module procedure mp_lsq_linear_no_b_nosd
          module procedure mp_lsq_linear_no_b_noa
       end interface

contains
       
       ! functions for fitting y=a(2)*x+a(1)
       
       subroutine mp_lsq_linear_sd(x,y,y_sd,yfit,a,a_sd)
              implicit none
              real, dimension(:), intent(in)  :: x, y     ! x values, y values
              real, dimension(:), intent(in)  :: y_sd     ! standard deviation in y values
              real, dimension(:), intent(out) :: yfit     ! values of lsq curve at xi
              real, dimension(2), intent(out) :: a, a_sd  ! coefficient vector, standard deviation

              call lsq_engine(x,y,y_sd,yfit,a,a_sd,mp_fj_linear)
       end subroutine mp_lsq_linear_sd

       subroutine mp_lsq_linear_nosd(x,y,yfit,a)
              implicit none
              real, dimension(:), intent(in)  :: x, y     ! x values, y values
              real, dimension(:), intent(out) :: yfit     ! values of lsq curve at xi
              real, dimension(2), intent(out) :: a        ! coefficient vector
              real, dimension(size(y))        :: y_sd     ! standard deviation in y values
              real, dimension(size(a))        :: a_sd     ! standard deviation in coefficients
              integer i

              do i=1,size(y)
                 y_sd(i) = 1
              end do

              call lsq_engine(x,y,y_sd,yfit,a,a_sd,mp_fj_linear)
       end subroutine mp_lsq_linear_nosd

       subroutine mp_lsq_linear_noa(x,y,yfit)
              implicit none
              real, dimension(:), intent(in)  :: x, y     ! x values, y values
              real, dimension(:), intent(out) :: yfit     ! values of lsq curve at xi
              real, dimension(2)              :: a        ! coefficient vector
              real, dimension(size(y))        :: y_sd     ! standard deviation in y values
              real, dimension(size(a))        :: a_sd     ! standard deviation in coefficients
              integer :: i

              y_sd=(/(1.,i=1,size(y))/)

              call lsq_engine(x,y,y_sd,yfit,a,a_sd,mp_fj_linear)
       end subroutine mp_lsq_linear_noa

       subroutine mp_fj_linear(xi,f)
              ! linear weighting function for use in lsq_engine, y=a(2)*x+a(1)
              implicit none
              real, intent(out), dimension(:,:) :: F      ! weight values (nterms,npts)
              real, intent(in), dimension(:)    :: XI     ! x values (npts)
              integer :: i

              DO I=1,size(F,2)
                 F(1,i)=1.
              end do

              DO I=1,size(F,2)
                 F(2,i)=xi(i)
              end do
       end subroutine mp_fj_linear
       
       ! functions for fitting y=a(1)*x
       
       subroutine mp_lsq_linear_no_b_sd(x,y,y_sd,yfit,a,a_sd)
              implicit none
              real, dimension(:), intent(in)  :: x, y     ! x values, y values
              real, dimension(:), intent(in)  :: y_sd     ! standard deviation in y values
              real, dimension(:), intent(out) :: yfit     ! values of lsq curve at xi
              real, dimension(1), intent(out) :: a, a_sd  ! coefficient vector, standard deviation

              call lsq_engine(x,y,y_sd,yfit,a,a_sd,mp_fj_linear_no_b)
       end subroutine mp_lsq_linear_no_b_sd

       subroutine mp_lsq_linear_no_b_nosd(x,y,yfit,a)
              implicit none
              real, dimension(:), intent(in)  :: x, y     ! x values, y values
              real, dimension(:), intent(out) :: yfit     ! values of lsq curve at xi
              real, dimension(1), intent(out) :: a        ! coefficient vector
              real, dimension(size(y))        :: y_sd     ! standard deviation in y values
              real, dimension(size(a))        :: a_sd     ! standard deviation in coefficients
              integer i

              do i=1,size(y)
                 y_sd(i) = 1
              end do

              call lsq_engine(x,y,y_sd,yfit,a,a_sd,mp_fj_linear_no_b)
       end subroutine mp_lsq_linear_no_b_nosd

       subroutine mp_lsq_linear_no_b_noa(x,y,yfit)
              implicit none
              real, dimension(:), intent(in)  :: x, y     ! x values, y values
              real, dimension(:), intent(out) :: yfit     ! values of lsq curve at xi
              real, dimension(1)              :: a        ! coefficient vector
              real, dimension(size(y))        :: y_sd     ! standard deviation in y values
              real, dimension(size(a))        :: a_sd     ! standard deviation in coefficients
              integer :: i

              y_sd=(/(1.,i=1,size(y))/)

              call lsq_engine(x,y,y_sd,yfit,a,a_sd,mp_fj_linear_no_b)
       end subroutine mp_lsq_linear_no_b_noa

       subroutine mp_fj_linear_no_b(xi,f)
              ! linear weighting function for use in lsq_engine, y=a(1)*x
              implicit none
              real, intent(out), dimension(:,:) :: F      ! weight values (nterms,npts)
              real, intent(in), dimension(:)    :: XI     ! x values (npts)
              integer :: i
              
              DO i=1,size(f,2)
                 f(1,i)=xi(i)
              end do
       end subroutine mp_fj_linear_no_b
end module linear_least_squares