program log_sinh_params
      implicit none

      integer :: i, n, iter
      real :: x(100), y(100), a, b, f, g, h, sumx, sumy, sumxy, sumxx
      real :: epsilon, delta, max_iter

      ! Define the size of the data series
      n = 100

      ! Read in the data series
      do i = 1, n
         read (unit=*, fmt=*) x(i), y(i)
      end do

      ! Initialize the parameter estimates
      a = 0.0
      b = 1.0

      ! Set the stopping criteria
      epsilon = 1.0e-6
      max_iter = 100

      ! Perform the least squares optimization
      do iter = 1, max_iter

         ! Initialize the sums
         sumx = 0.0
         sumy = 0.0
         sumxy = 0.0
         sumxx = 0.0

         ! Compute the sums
         do i = 1, n
            f = sinh(a + b*x(i))
            g = 1/b * log(f)
            h = f / (f + 1)
            sumx = sumx + (y(i) - g) * h
            sumy = sumy - (y(i) - g) / b
            sumxy = sumxy + (y(i) - g) * x(i) * h
            sumxx = sumxx + x(i)**2 * h
         end do

         ! Update the parameter estimates
         delta = sumx / sumxx
         a = a - delta * sumy
         b = b - delta * sumxy

         ! Check for convergence
         if (abs(delta) < epsilon) exit

      end do

      ! Output the estimated parameters
      write (unit=*, fmt=*) "Estimated a: ", a
      write (unit=*, fmt=*) "Estimated b: ", b

      end program log_sinh_params