Geometry/parabola_ex.f90 at master · lguez/Geometry · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
subroutine parabola_ex ( x1, y1, x2, y2, x3, y3, x, y, ierror )

!*****************************************************************************80
!
!! PARABOLA_EX: extremal point of a parabola determined by three points.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    02 November 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision X1, Y1, X2, Y2, X3, Y3, the coordinates of
!    three points on the parabola.  X1, X2 and X3 must be distinct.
!
!    Output, double precision X, Y, the X coordinate of the extremal point
!    of the parabola, and the value of the parabola at that point.
!
!    Output, integer IERROR, error flag.
!    0, no error.
!    1, two of the X values are equal.
!    2, the data lies on a straight line; there is no finite extremal point.
!
  implicit none

  double precision bot
  integer ierror
  double precision x
  double precision x1
  double precision x2
  double precision x3
  double precision y
  double precision y1
  double precision y2
  double precision y3

  ierror = 0

  if ( x1 == x2 .or. x2 == x3 .or. x3 == x1 ) then
    ierror = 1
    return
  end if

  if ( y1 == y2 .and. y2 == y3 .and. y3 == y1 ) then
    x = x1
    y = y1
    return
  end if

  bot = ( x2 - x3 ) * y1 - ( x1 - x3 ) * y2 + ( x1 - x2 ) * y3

  if ( bot == 0.0D+00 ) then
    ierror = 2
    return
  end if

  x = 0.5D+00 * ( x1 * x1 * ( y3 - y2 ) &
                + x2 * x2 * ( y1 - y3 ) &
                + x3 * x3 * ( y2 - y1 ) ) / bot

  y = (    ( x  - x2 ) * ( x  - x3 ) * ( x2 - x3 ) * y1 &
         - ( x  - x1 ) * ( x  - x3 ) * ( x1 - x3 ) * y2 &
         + ( x  - x1 ) * ( x  - x2 ) * ( x1 - x2 ) * y3 ) / &
         ( ( x1 - x2 ) * ( x2 - x3 ) * ( x1 - x3 ) )

  return
end