! $Id: sqfit.f90 2022-11-02 09:51:28Z ychen $

!****s* fsi/sqfit *
!
! NAME
!   sqfit  - This subroutine performs quadratic least squares fit
!
! SYNOPSIS
!    call sqfit(x, y, a, b, c)
!
! DESCRIPTION
!    This subroutine performs quadratic least squares fit
!              y=a*x^2 + b*x +c
!
! INPUTS
!    REAL(wp), DIMENSION(:)   :: x    array of input argument data
!    REAL(wp), DIMENSION(:)   :: y    array of input function data
!
! OUTPUT
!    REAL(wp)                 :: a
!    REAL(wp)                 :: b
!    REAL(wp)                 :: c
!
! AUTHOR
!   S.V.Sokolovskiy, UCAR
!   Update: Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov
!
! COPYRIGHT
!   Copyright (c) 2022-2023 Yong Chen <yong.chen@noaa.gov>
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE sqfit(x, y, a, b, c)

  USE typesizes, ONLY: wp => EightByteReal

  IMPLICIT NONE

  REAL(wp), DIMENSION(:),    INTENT(IN)    :: x
  REAL(wp), DIMENSION(:),    INTENT(IN)    :: y
  REAL(wp),                  INTENT(OUT)   :: a
  REAL(wp),                  INTENT(OUT)   :: b
  REAL(wp),                  INTENT(OUT)   :: c

  INTEGER  :: n, i
  REAL(wp) :: a11, a12, a13, a21, a22, a23, a31, a32, a33
  REAL(wp) :: b1,  b2,  b3
  REAL(wp) :: det, det1, det2, det3

  n = SIZE(x)

  a11 = 0.0_wp
  a21 = 0.0_wp
  a31 = 0.0_wp
  a32 = 0.0_wp
  a33 = 0.0_wp
  b1 = 0.0_wp
  b2 = 0.0_wp
  b3 = 0.0_wp

  DO i = 1, n
     a11 = a11 + x(i)**4
     a21 = a21 + x(i)**3
     a31 = a31 + x(i)**2
     a32 = a32 + x(i)
     a33 = a33 + 1.0_wp
     b1 = b1 + y(i) * x(i)**2
     b2 = b2 + y(i) * x(i)
     b3 = b3 + y(i)
  ENDDO

  a12 = a21
  a13 = a31
  a22 = a31
  a23 = a32

  det  = a11 * a22 * a33 + a13 * a21 * a32 + a12 * a23 * a31 - a13 * &
        a22 * a31 - a11 * a23 * a32 - a12 * a21 * a33

  det1 = b1 * a22 * a33 + a13 * b2 * a32 + a12 * a23 * b3 - a13 *   &
        a22 * b3 - b1 * a23 * a32 - a12 * b2 * a33

  det2 = a11 * b2 * a33 + a13 * a21 * b3 + b1 * a23 * a31 - a13 *   &
         b2 * a31 - a11 * a23 * b3 - b1 * a21 * a33

  det3 = a11 * a22 * b3 + b1 * a21 * a32 + a12 * b2 * a31 - b1 *    &
         a22 * a31 - a11 * b2 * a32 - a12 * a21 * b3

  a = det1 / det
  b = det2 / det
  c = det3 / det

END SUBROUTINE sqfit