!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
! D.J. Jeffery, 2020jan06
! Use the quadratic Lagrange interpolation formula:
!  https://en.wikipedia.org/wiki/Lagrange_polynomial 
!  http://mathonline.wikidot.com/quadratic-lagrange-interpolating-polynomials
!  https://math.stackexchange.com/questions/889569/finding-a-parabola-from-three-points-algebraically
! Also the numerically good quadratic solution of Numerical Recipies (1992, p. 178).
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
      subroutine xne_quadratic(pair,x_q) 
      use precision_set
      use contants
      include 'module_implicit.f'
!      implicit none
      real (kind=np) :: a,b,c
      real (kind=np) :: aa(3)
      real (kind=np), intent(in) :: pair(2,3)  
      real (kind=np) :: q
      real (kind=np) :: ya,yc
      real (kind=np), intent(out) :: x_q  
!
      aa(1)=pair(2,1)/( (pair(1,1)-pair(1,2))*(pair(1,1)-pair(1,3)) )
      aa(2)=pair(2,2)/( (pair(1,2)-pair(1,1))*(pair(1,2)-pair(1,3)) )
      aa(3)=pair(2,3)/( (pair(1,3)-pair(1,1))*(pair(1,3)-pair(1,2)) )
      a=sum(aa(:))
      b=-(pair(1,2)+pair(1,3))*aa(1)                                    &
     &  -(pair(1,1)+pair(1,3))*aa(2)                                    &
     &  -(pair(1,1)+pair(1,2))*aa(3)                                    &
     &  -1.0_np
      c= (pair(1,2)*pair(1,3))*aa(1)                                     &
     &  +(pair(1,1)*pair(1,3))*aa(2)                                     &
     &  +(pair(1,1)*pair(1,2))*aa(3)
!      
      q=-sign(1.0_np,b)*0.5_np*(abs(b)+sqrt(b**2-4.0_np*a*c))
      ya=q/a
      yc=c/q
      if(a .gt. 0.0_np) then
          x_q=min(ya,yc)    ! This upward parabola case.
        else
          x_q=max(ya,yc)    ! This downward parabola case.
      end if
!      b=b+1.0_np
!      print*,'In xne_quadratic'
!      print*,a,b,c
!      print*,a,ya,yc
!      print*,a*x_q**2+b*x_q+c
!      print*,a*ya**2+b*ya+c
!      print*,a*yc**2+b*yc+c
!
      end subroutine xne_quadratic 
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12