!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! BstFit
!
! This program provides the user with three options:
! (1) Compute least-squares lines using vertical, horizontal, and
!     perpendicular distances
! (2) Show the distances between points and a given line
! (3) Compute linear, exponential, and power least-squares curves
! 
! Points are read from an ASCII text file, where each line contains a x-value
! followed by a space and a y-value.  Values are expected to be real numbers.
! The number of points is not specified; the program reads data from each line
! until the end of the file (or 1,000 points, whichever comes first).
!
! Output is displayed to the screen.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module Mdataset

  ! Declare all variables
  implicit none

  ! Declare a Tdataset type to hold data points, count, sums, etc.
  type :: Tdataset
    real             :: datapt(4,1000)    ! (x,y) data points from input file
                                          !    stored as (x, y, ln x, ln y)
    integer          :: npts              ! Number of points in input file
    real             :: valsum(4)         ! Sums of values of
                                          !    x, y, ln x, and ln y
    real             :: sqsum(4)          ! Sums of values of x**2, y**2,
                                          !    (ln x)**2, and (ln y)**2
    real             :: xysum(3)          ! Sums of products
                                          !    (xy, x ln y, ln x ln y)
    logical          :: has_np(2)         ! Are any x- or y-values nonpositive?
  end type Tdataset

end module Mdataset

program BstFitProj

  use Mdataset

  ! Declare all variables
  implicit none

  ! Type                Name(s)             Description
  integer			 :: if
  integer            :: ichoice           ! The user's choice from the menu
  type(Tdataset)     :: ds                ! Points and computed values
  real               :: aval, bval        ! Values of a and b (in y=ax+b)
  real               :: slopem, intrcpt   ! Slope and y-intercept of a line
                                          !   supplied by the user
  real               :: ssvd, sshd, sspd  ! Sum of squares of horiz, vert, and
                                          !   perpendicular distances
  logical            :: suces             ! Exponential or power fit success?
  integer            :: outu              ! Output file unit number

  ! Introduce the program
  write (*, *) 'This program performs least-squares fitting on ', &
                            'a set of n (x,y) pairs'
  write (*, *) ! Blank line

  ! Ask the user for the input file, and read the data from it  
  call LoadDataSet(ds)
  write (*, *) ! Blank line
  
  ! Show the list of points
  write (*, *) 'Data points are:'
  call ShowDataPoints(ds)
  write (*, *) ! Blank line
  
  ! Ask the user for the name of the output file, and open it
  outu = OpenOutput()
  write (*, *) ! Blank line
  
  ! The user hasn't chosen anything from the menu yet
  ichoice = 0
  
  ! Let the user continue choosing from the menu until s/he chooses to exit
  do while (ichoice .ne. 4)
  
    ! Display a list of options, and allow the user to select one
    write (*, *) ! Blank line
    write (*, *) '===================================================='
    write (*, *) 'MENU'
    write (*, *) '===================================================='
    write (*, *) '1. Compute three linear least squares lines'
    write (*, *) '2. Compute distances between points and a given line'
    write (*, *) '3. Find linear, exponential, and power least squares'
    write (*, *) '4. Exit'
    write (*, *) '===================================================='
    write (*, *) ! Blank line
    write (*, *) 'Choose an option: '
    read  (*, *) ichoice
    
    if (if .eq. 5) then
      print *, 'Why did you name a variable "if" anyway?'
    end if

    select case (ichoice)
      case (1) ! Option 1 - Three linear least squares lines

        write (outu, *) 'THREE LEAST SQUARES LINES:'

        ! Linear least squares, vertical distances
        call LinearLeastSquares(ds, 1, 2, aval, bval)
        write (outu, *) ! Blank line
        write (outu, *) 'Using vertical distances:'
        write (outu, *) 'y = ', trim(real2str(aval)), &
                                     'x + ', trim(real2str(bval))

        ! Linear least squares, horizontal distances
        call LinearLeastSquares(ds, 2, 1, aval, bval)
        write (outu, *) ! Blank line
        write (outu, *) 'Using horizontal distances:'
        write (outu, *) 'y = ', trim(real2str(aval)), &
                                     'x + ', trim(real2str(bval))

        ! Linear least squares, perpendicular distances
        call PerpLeastSquares(ds, 0, aval, bval)
        write (outu, *) ! Blank line
        write (outu, *) 'Using perpendicular distances:'
        write (outu, *) 'y = ', trim(real2str(aval)), &
                                     'x + ', trim(real2str(bval))
        
        write (outu, *) '--OR--'

        call PerpLeastSquares(ds, 1, aval, bval)
        write (outu, *) 'y = ', trim(real2str(aval)), &
                                     'x + ', trim(real2str(bval))

        write (outu, *) ! Blank line

      case (2) ! Option 2 - Distances between points and a given line

        write (outu, *) 'SUM OF SQUARES OF DISTANCES BETWEEN ', &
                                     'POINTS AND A GIVEN LINE:'
      
        write (*, *) ! Blank line
        write (*, *) 'Enter the line''s slope and y-intercept: '
        read (*, *) slopem, intrcpt
        
        write (outu, *) ! Blank line
        write (outu, *) 'Using f(x) = ', trim(real2str(slopem)), &
                                     'x + ', trim(real2str(intrcpt))

        write (outu, *) ! Blank line
        call SquareDistances(ds, slopem, intrcpt, ssvd, sshd, sspd)

        write (outu, *) 'Sum of squares of vertical distances:      ', &
                                      ssvd
        write (outu, *) ! Blank line
        write (outu, *) 'Sum of squares of horizontal distances:    ', &
                                      sshd
        write (outu, *) ! Blank line
        write (outu, *) 'Sum of squares of perpendicular distances: ', &
                                      sspd
        write (outu, *) ! Blank line

      case (3) ! Option 3 - Linear, exponential, and power least squares

        write (outu, *) 'LINEAR, EXPONENTIAL, AND POWER LEAST SQUARES:'
      
        call LinearLeastSquares(ds, 1, 2, aval, bval)
        write (outu, *) ! Blank line
        write (outu, *) 'Linear fit:       y = ', trim(real2str(aval)), &
                                     'x + ', trim(real2str(bval))

        call ExpLeastSquares(ds, aval, bval, suces)
        write (outu, *) ! Blank line
        if (suces .eqv. .true.) then
          write (outu, *) 'Exponential fit:  y = ', trim(real2str(bval)), &
                                       ' * e^(', trim(real2str(aval)), 'x)'
        else
          write (outu, *) 'Unable to perform exponential fit: ', &
                                       'there is at least one nonpositive y-value'
        end if

        call PowerLeastSquares(ds, aval, bval, suces)
        write (outu, *) ! Blank line
        if (suces .eqv. .true.) then
          write (outu, *) 'Power fit:        y = ', trim(real2str(bval)), &
                                       ' * x^', trim(real2str(aval))
        else
          write (outu, *) 'Unable to perform power fit: ', &
                                       'there is at least one nonpositive x- or y-value'
        end if

        write (outu, *) ! Blank line
        
      case (4) ! Option 4 - Exit
        ! Allow loop to stop

      case default
        write (*, *) 'The option you chose is not valid.'
    end select
    
  end do
  
  stop
  
contains
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! LoadDataSet - Ask the user for an input filename, and read data points from
  !               that file into the ds%datapt array, storing the number
  !               of points in ds%npts and computing sums
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine LoadDataSet(ds)
  
    ! Declare all variables
    implicit none
    ! Type                         Name(s)             Description
    type(Tdataset), intent(out) :: ds                ! Data points and sums
    character (len=64)          :: infile            ! Input filename
    integer                     :: ieof              ! File status
    integer                     :: indx              ! Loop counter
    integer                     :: inu               ! Input file unit number
    integer                     :: imaxpairs         ! Maximum number of pairs
    
    ! Default maximum number of pairs to array size
    imaxpairs = ubound(ds%datapt,2)
    
    ! Ask the user for the name of the input file
    write (*, *) 'Enter the input filename or - for keyboard input: '
    read (*, *) infile
    
    ! Determine the input unit number
    if (infile .eq. '-') then
    
      inu = 5     ! Read from keyboard (standard input)
      
      write (*, *) 'How many pairs? '
      read (*, *) imaxpairs
      
      ! Make sure this number is legal
      if (imaxpairs .lt. 1 .or. imaxpairs .gt. ubound(ds%datapt,2)) then
        write (*, *) 'Number of pairs must be between 1 and ', &
          trim(int2str(ubound(ds%datapt,2)))
        stop
      end if
      
      write (*, *) 'Enter pairs, one per line'
      write (*, *) ! End line before input
      
    else
    
      inu = 10    ! Read from file
    
      ! Open the file
      open (10, file = infile, status = 'old', iostat = ieof)
      if (ieof .ne. 0) then
        write (*, *) 'The specified file could not be opened.'
        stop
      end if
      
    end if

    ! Initialize the members of the dataset structure
    ds%npts      = 0
    ds%valsum(1) = 0.0
    ds%valsum(2) = 0.0
    ds%valsum(3) = 0.0
    ds%valsum(4) = 0.0
    ds%sqsum(1)  = 0.0
    ds%sqsum(2)  = 0.0
    ds%sqsum(3)  = 0.0
    ds%sqsum(4)  = 0.0
    ds%xysum(1)  = 0.0
    ds%xysum(2)  = 0.0
    ds%xysum(3)  = 0.0
    ds%has_np(1) = .false.
    ds%has_np(2) = .false.

    ! Read data points until EOF or imaxpairs reached
    do while (ds%npts < imaxpairs)
    
      ! Load the points into the datapt array, 
      read (inu, *, end = 110) &
        ds%datapt(1, ds%npts + 1), &
        ds%datapt(2, ds%npts + 1)
      
      ! Calculate ln x, or set has_np if the value is nonpositive
      if (ds%datapt(1, ds%npts + 1) .le. 0) then
        ds%datapt(3, ds%npts + 1) = 0
        ds%has_np(1) = .true.
      else
        ds%datapt(3, ds%npts + 1) = log(ds%datapt(1, ds%npts + 1))
      end if
      
      ! Calculate ln y, or set has_np if the value is nonpositive
      if (ds%datapt(2, ds%npts + 1) .le. 0) then
        ds%datapt(4, ds%npts + 1) = 0
        ds%has_np(2) = .true.
      else
        ds%datapt(4, ds%npts + 1) = log(ds%datapt(2, ds%npts + 1))
      end if
        
      ! Update the number of points
      ds%npts      = ds%npts      + 1

      ! Update the value sums and square sums
      do indx = 1,4
        ds%valsum(indx) = ds%valsum(indx) + ds%datapt(indx, ds%npts)
        ds%sqsum(indx)  = ds%sqsum(indx)  + ds%datapt(indx, ds%npts)**2
      end do
      
      ! Update the xy sums (xy, x ln y, ln x ln y)
      ds%xysum(1)  = ds%xysum(1) + &
        ds%datapt(1, ds%npts) * ds%datapt(2, ds%npts)
      ds%xysum(2)  = ds%xysum(2) + &
        ds%datapt(1, ds%npts) * ds%datapt(4, ds%npts)
      ds%xysum(3)  = ds%xysum(3) + &
        ds%datapt(3, ds%npts) * ds%datapt(4, ds%npts)
      
    end do
    
    ! If we get here and we're reading from a file, there are more than
    ! ubound(ds%datapt,2) points in the file, so display a warning
    if (inu .eq. 10) then
      write (*, *) 'WARNING: Only the first ', &
        trim(int2str(ubound(ds%datapt,2))), ' points will be used'
    end if

    ! We have read all points from the file; close it
110 if (inu .eq. 10) then
      close (10)
    end if
    
    return

  end subroutine LoadDataSet
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! ShowDataPoints - Display the list of data points
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine ShowDataPoints(ds)
  
    ! Declare all variables
    implicit none
    ! Type                         Name(s)      Description
    type(Tdataset), intent(in)  :: ds         ! Data points and sums
    integer                     :: indx       ! Loop counter
  
    ! Display a header
    write (*, '(a5,a2,2a10)') '=====', '==', '==========', &
                                                          '=========='
    write (*, '(a5,2x,2a10)') 'i', 'x_i', 'y_i'
    write (*, '(a5,a2,2a10)') '=====', '==', '==========', &
                                                          '=========='

    ! Display the list of points
    do indx = 1,ds%npts
      write (*, '(i5,a,2f10.4)') &
        indx, ': ', ds%datapt(1,indx), ds%datapt(2,indx)
    end do
    
    ! Display a footer
    write (*, '(a5,a2,2a10)') '=====', '==', '==========', &
                                                          '=========='
    
    return
    
  end subroutine ShowDataPoints
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! SquareDistances - Calculate the sums of squares of distances of a set of
  !                   points from a given line f(x) = mx + b
  !
  ! finv denotes the inverse of the function f
  !
  ! Sum of squares of vertical distances is returned in ssvd
  ! Sum of squares of horizontal distances is returned in sshd
  ! Sum of squares of perpendicular distances is returned in sspd
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine SquareDistances(ds, slopem, intrcpt, ssvd, sshd, sspd)
  
    ! Declare all variables
    implicit none
    ! Type                         Name(s)      Description
    type(Tdataset), intent(in)  :: ds         ! Data points and sums
    real, intent(in)            :: slopem     ! Slope of the line
    real, intent(in)            :: intrcpt    ! Line's y-intercept
    real                        :: fx         ! f(x) at a given x
    real                        :: finvy      ! finv(y) at a given y
    real                        :: perpx      ! x-coord where a perpendicular
                                              !   line through (x_i,y_i)
                                              !   intersects
    real                        :: perpy      ! y-coord of intersection
    real :: vdist, hdist, pdist               ! Distances of a particular point
    real, intent(out)           :: ssvd       ! Sum of squares of vert dist's
    real, intent(out)           :: sshd       ! Sum of squares of horiz dist's
    real, intent(out)           :: sspd       ! Sum of squares of perp dist's
    integer                     :: indx       ! Loop counter
    
    ! Initialize sums to 0.0
    ssvd = 0.0
    sshd = 0.0
    sspd = 0.0
  
    ! Add the squares of the distances
    do indx = 1,ds%npts
      ! Calculate f(x_i) = m x_i + b
      fx = slopem * ds%datapt(1,indx) + intrcpt
      
      ! Calculate finv(y_i) = (y_i - b) / m
      finvy = (ds%datapt(2,indx) - intrcpt) / slopem
      
      ! Vertical distance = |f(x_i) - y_i|
      vdist = abs(fx - ds%datapt(2,indx))
      
      ! Horizontal distance = |finv(y_i) - x_i|
      hdist = abs(finvy - ds%datapt(1,indx))
      
      ! Perpendicular distance is more complicated
      
      ! The perpendicular line through (x_i, y_i) will intersect
      ! at x-coordinate (y_i + x_i/m - b)/(m + 1/m)
      ! and y-coordinate f(perpx)
      perpx = (ds%datapt(2,indx) + ds%datapt(1,indx)/slopem - intrcpt) &
              / (slopem + 1.0/slopem)
      perpy = slopem * perpx + intrcpt
      
      ! Perpendicular distance is now given by the distance formula
      !    ________________________________
      !   /              2                2
      ! \/  (y_i - perpy)  + (x_i - perpx)
      !
      pdist = sqrt((ds%datapt(2,indx) - perpy)**2 + (ds%datapt(1,indx) - perpx)**2)
      
      ! Now add these distances squared to the sums
      ssvd = ssvd + vdist**2
      sshd = sshd + hdist**2
      sspd = sspd + pdist**2
    end do
    
    return
    
  end subroutine SquareDistances
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! OpenOutput - Open an output file and return the unit number
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  integer function OpenOutput() result(outu)

    ! Declare all variables
    implicit none
    ! Type                  Name(s)             Description
    character (len=64) :: outfile           ! Output filename
    integer            :: ieof              ! File status

    ! Ask the user for the name of the output file
    write (*, *) 'Enter the output filename or - for screen display: '
    read (*, *) outfile
    
    if (outfile .eq. '-') then
      outu = 6  ! Standard output
    else
    
      ! Open the file
      open (10, file = outfile, status = 'replace', iostat = ieof)
      if (ieof .ne. 0) then
        write (*, *) 'The specified file could not be opened.'
        stop
      end if
      
      outu = 10 ! File output
    end if
    
    return
    
  end function OpenOutput
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! LinearLeastSquares - Compute a linear least squares line y = ax + b
  !                      using vertical or horizontal distances
  !
  ! For vertical distances, set xcol = 1 and ycol = 2
  ! For horizontal distances, set xcol = 2 and ycol = 1
  ! Other values of xcol and ycol are used by the power and exponential
  !   least squares functions (see below)
  ! Values of a and b will be stored in aval and bval
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine LinearLeastSquares(ds, xcol, ycol, aval, bval)
  
    ! Declare all variables
    implicit none
    ! Type                         Name(s)      Description
    type(Tdataset), intent(in)  :: ds         ! Data points and sums
    integer, intent(in)         :: xcol       ! Which column contains x-values
                                              !   (independent data)?
    integer, intent(in)         :: ycol       ! Which column contains y-values
                                              !   (dependent data)?
    integer                     :: xycol      ! Which xyval to use?
    real, intent(out)           :: aval       ! Value of a
    real, intent(out)           :: bval       ! Value of b
    real                        :: dval       ! Value of delta (denominator)
    
    ! Determine which xyval to use (xy, x ln y, or ln x ln y)
    if (xcol .eq. 1 .and. ycol .eq. 2 .or. xcol .eq. 2 .and. ycol .eq. 1) then
      xycol = 1
    else if (xcol .eq. 1 .and. ycol .eq. 4 .or. xcol .eq. 4 .and. ycol .eq. 1) then
      xycol = 2
    else if (xcol .eq. 3 .and. ycol .eq. 4 .or. xcol .eq. 4 .and. ycol .eq. 3) then
      xycol = 3
    else
      write (*, *) 'ERROR: Invalid call to LinearLeastSquares'
      write (*, *) 'Unable to continue'
      stop
    end if
    
    ! Compute the denominator:
    !
    !             n                n        2
    !            ___    2      (  ___     )
    ! delta = n  \    x    -  (   \    x   )
    !            /__   i      (   /__   i  )
    !                          (          )
    !            i=1              i=1
    
    dval = ds%npts * ds%sqsum(xcol) - ds%valsum(xcol)**2
    
    ! Compute a:
    !
    !          n            n        n
    !         ___          ___      ___
    !      n  \    x y  -  \    x   \    y
    !         /__   i i    /__   i  /__   i
    !
    !         i=1          i=1      i=1
    ! a = ----------------------------------
    !                   delta
    
    aval = ( ds%npts * ds%xysum(xycol) -         &
             ds%valsum(xcol) * ds%valsum(ycol) ) &
           / dval
    
    
    ! Compute b:
    !
    !       n        n            n        n
    !      ___      ___    2     ___      ___
    !      \    y   \    x    -  \    x   \    x y
    !      /__   i  /__   i      /__   i  /__   i i
    !
    !      i=1      i=1          i=1      i=1
    ! b = ------------------------------------------
    !                   delta
    
    bval = ( ds%valsum(ycol) * ds%sqsum(xcol) -   &
             ds%valsum(xcol) * ds%xysum(xycol)  ) &
           / dval
           
    ! If we are calculating using horizontal distances, we need to
    ! solve x=ay+b for y (y = 1/a * x - b/a) and change a and b
    ! accordingly
    if (xcol .gt. ycol) then
      aval = 1.0 / aval
      bval = -bval * aval
    end if
    
    return
    
  end subroutine LinearLeastSquares
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! PerpLeastSquares - Compute a linear least squares line y = ax + b
  !                    using perpendicular distances
  !
  ! To use a = -p + sqrt(...), set subtr to an even value
  ! To use a = -p - sqrt(...), set subtr to an odd value
  ! Values of a and b will be stored in aval and bval
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine PerpLeastSquares(ds, subtr, aval, bval)
  
    ! Declare all variables
    implicit none
    ! Type                         Name(s)      Description
    type(Tdataset), intent(in)  :: ds         ! Data points and sums
    integer, intent(in)         :: subtr      ! Determines which of two values
                                              !   of a will be chosen
    real, intent(out)           :: aval       ! Value of a
    real, intent(out)           :: bval       ! Value of b
    real                        :: pval       ! Value of p
    
    ! Compute p:
    !
    !         (   n           n       )    [  (  n      )  2   (  n      )  2 ]
    !        (   ___    2    ___    2  )   [ (  ___      )    (  ___      )   ]
    !      n (   \    x   -  \    y    ) - [ (  \    x   )  - (  \    y   )   ]
    !        (   /__   k     /__   k   )   [ (  /__   k  )    (  /__   k  )   ]
    !        (                         )   [  (         )      (         )    ]
    !         (  k=1         k=1      )    [    k=1              k=1          ]  
    ! p = ----------------------------------------------------------------------
    !         (     n           n       n     )
    !        (     ___         ___     ___     )
    !      2 (  n  \   x y  -  \   x   \   y   )
    !        (     /__  k k    /__  k  /__  k  )
    !        (                                 )
    !         (    k=1         k=1     k=1    )
    
    pval = (ds%npts*(ds%sqsum(1)-ds%sqsum(2)) - (ds%valsum(1)**2 - ds%valsum(2)**2)) &
           / (2 * (ds%npts*ds%xysum(1) - ds%valsum(1)*ds%valsum(2)))
    
    !                        _______
    !                       / 2
    ! Compute a = -p +/-   / p  + 1
    !                    \/
    
    aval = -pval + (-1)**subtr * sqrt(pval**2 + 1)

    !               n           n
    !              ___         ___
    !              \    y  - a \    x
    !              /__   k     /__   k
    !
    !              k=1         k=1
    ! Compute b = ---------------------
    !                       n
    
    bval = ( ds%valsum(2) - aval * ds%valsum(1) ) / ds%npts
    
    return
    
  end subroutine PerpLeastSquares
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                     ax
  ! ExpLeastSquares - Compute an exponential least squares curve y = b e
  !                   using vertical distances
  !
  ! Values of a and b will be stored in aval and bval
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine ExpLeastSquares(ds, aval, bval, suces)
  
    ! Declare all variables
    implicit none
    ! Type                         Name(s)      Description
    type(Tdataset), intent(in)  :: ds         ! Data points and sums
    real, intent(out)           :: aval       ! Value of a
    real, intent(out)           :: bval       ! Value of b
    logical, intent(out)        :: suces      ! Success?

    ! Make sure all y-values are positive
    if (ds%has_np(2) .eqv. .true.) then
      suces = .false.
      return
    end if
    
    ! Exponential least squares is done by reducing the problem to a linear
    ! one, which is more easily solved.
    
    ! We want    y  = b e^(ax)
    ! ==>     ln(y) = ln(b e^(ax)
    ! ==>     ln(y) = ln(b) + ax
    ! ==>     ln(y) = ax + ln(b)
    
    ! Use linear least squares on ln(y) = ax + ln(b)
    call LinearLeastSquares(ds, 1, 4, aval, bval)
    ! The variable bval now contains ln(b); make it contain b instead
    bval = exp(bval)
    
    ! Now     ln(y) = ax + ln(b)
    ! So         y  = e^(ax + ln(b))
    ! ==>        y  = e^(ax) e^(ln(b))
    ! ==>        y  = e^(ax) b
    ! ==>        y  = b e^(ax)
    
    suces = .true.
    return
  
  end subroutine ExpLeastSquares
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                           a
  ! PowerLeastSquares - Compute an least squares curve y = b x
  !                     using vertical distances
  !
  ! Values of a and b will be stored in aval and bval
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine PowerLeastSquares(ds, aval, bval, suces)
  
    ! Declare all variables
    implicit none
    ! Type                         Name(s)      Description
    type(Tdataset), intent(in)  :: ds         ! Data points and sums
    real, intent(out)           :: aval       ! Value of a
    real, intent(out)           :: bval       ! Value of b
    logical, intent(out)        :: suces      ! Success?

    ! Make sure all values are positive
    if ((ds%has_np(1) .eqv. .true.) .or. (ds%has_np(2) .eqv. .true.)) then
      suces = .false.
      return
    end if
    
    ! Exponential least squares is done by reducing the problem to a linear
    ! one, which is more easily solved.
    
    ! We want    y  = b x^a
    ! ==>     ln(y) = ln(b x^a)
    ! ==>     ln(y) = ln(b) + ln(x^a)
    ! ==>     ln(y) = ln(b) + a ln(x)
    ! ==>     ln(y) = a ln(x) + ln(b)
    
    ! Use linear least squares on ln(y) = a ln(x) + ln(b)
    call LinearLeastSquares(ds, 3, 4, aval, bval)
    ! The variable bval now contains ln(b); make it contain b instead
    bval = exp(bval)
    
    ! Now     ln(y) = a ln(x) + ln(b)
    !
    !                  a ln(x) + ln(b)
    ! ==>        y  = e
    !
    !                  a ln(x)  ln(b)
    ! ==>        y  = e        e
    !
    !                   a ln(x)
    ! ==>        y  = (e )      b
    !
    !                   ln(x) a
    ! ==>        y  = (e     )  b
    !
    !                  a
    ! ==>        y  = x  b
    !
    !                    a
    ! ==>        y  = b x
    
    suces = .true.
    return
  
  end subroutine PowerLeastSquares
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! int2str - Convert an integer i to a string
  !
  ! The string will need to be trim()ed before displaying
  !
  ! Shouldn't something like this be built in?!
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  character(len=32) function int2str(i) result(s)
  
    ! Declare all variables
    implicit none
    
    ! Type                 Name       Description
    integer, intent(in) :: i        ! The number to convert
    integer             :: ival     ! The number, with leading digits removed
    integer             :: ndigits  ! The number of digits in the number
    integer             :: indx     ! Loop counter
    integer             :: is_neg   ! 0 if negative, 1 if positive

    ! Start with a 0 string, overwriting the 0 if necessary
    s = '0'
    
    ! Handle the cases where i is <= 0
    if (i .eq. 0) then
      return
    else if (i .lt. 0) then
      s(1:1) = '-'  ! Prefix output with a minus sign
      is_neg = 1    ! There is one extra character at the start of the number
      ival = -i     ! Now write out the absolute value of i
    else
      is_neg = 0    ! It's not negative
      ival = i      ! We want to write out i
    end if
    
    ! The number of digits needed to write out an integer i in base b is
    ! |           |
    ! |  log (i)  | + 1
    ! |_    b    _|
    !
    ndigits = floor(log10(ival*1.0)) + 1
    
    ! Write out the digits in order, starting at the left
    do indx = 1, ndigits
      ! What is the first digit's ASCII code?  Append it to the string
      s(indx+is_neg:indx+is_neg) = char(48 + ival/10**(ndigits-indx))
      ! Now remove the first digit from ival
      ival = ival - (ival/10**(ndigits-indx))*10**(ndigits-indx)
    end do
  
    return
    
  end function int2str
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! real2str - Convert a real number r to a string
  !
  ! The string will need to be trim()ed before displaying
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  character(len=128) function real2str(r) result(s)
  
    ! Declare all variables
    implicit none
    
    ! Type                 Name       Description
    real, intent(in)    :: r        ! The number to convert
    real                :: rval     ! Value of |r|
    integer             :: intpt    ! The integer part of the number
    integer             :: ifracpt  ! The fractional part of the number
                                    !   (rounded somewhat)
    character(len=1)    :: is_neg   ! '-' if negative, ' ' otherwise
    character(len=32)   :: sintpt   ! The integer part, as a string
    character(len=32)   :: slead0   ! Leading 0s for fractional part
    character(len=32)   :: sfracpt  ! The fractional part, as a string
    integer             :: nfdigs=4 ! Number of fractional digits
    
    ! If the number's negative, remember this, but make it positive
    if (r .lt. 0) then
      is_neg = '-'
      rval = -r
    else
      is_neg = ' '
      rval = r
    end if

    ! Integer part = |_ r _|
    intpt = floor(rval)
    
    ! To get the fractional part, subtract off the integer part, multiply
    ! by 10**(# decimal places), and take the floor to drop off extra digits
    ifracpt = floor((rval-intpt) * 10**nfdigs)
    
    ! Now convert the integer and fractional parts to strings
    sintpt = int2str(intpt)
    sfracpt = int2str(ifracpt)
    slead0 = repeat('0', nfdigs-len_trim(sfracpt))
    
    ! Now just concatenate the number together
    s = trim(is_neg) // &  ! Negative sign (if any)
        trim(sintpt) // &  ! Integer part
        '.' // &           ! Decimal point
        trim(slead0) // &  ! Fractional part leading zeros
        int2str(ifracpt)   ! Fractional part
    
    return
    
  end function real2str
  
end program BstFitProj