1 | |
---|
2 | ! |
---|
3 | ! This program may be freely redistributed under the |
---|
4 | ! condition that the copyright notices (including this |
---|
5 | ! entire header) are not removed, and no compensation |
---|
6 | ! is received through use of the software. Private, |
---|
7 | ! research, and institutional use is free. You may |
---|
8 | ! distribute modified versions of this code UNDER THE |
---|
9 | ! CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE |
---|
10 | ! TO IT IN THE SAME FILE REMAIN UNDER COPYRIGHT OF THE |
---|
11 | ! ORIGINAL AUTHOR, BOTH SOURCE AND OBJECT CODE ARE |
---|
12 | ! MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR |
---|
13 | ! NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution |
---|
14 | ! of this code as part of a commercial system is |
---|
15 | ! permissible ONLY BY DIRECT ARRANGEMENT WITH THE |
---|
16 | ! AUTHOR. (If you are not directly supplying this |
---|
17 | ! code to a customer, and you are instead telling them |
---|
18 | ! how they can obtain it for free, then you are not |
---|
19 | ! required to make any arrangement with me.) |
---|
20 | ! |
---|
21 | ! Disclaimer: Neither I nor: Columbia University, the |
---|
22 | ! National Aeronautics and Space Administration, nor |
---|
23 | ! the Massachusetts Institute of Technology warrant |
---|
24 | ! or certify this code in any way whatsoever. This |
---|
25 | ! code is provided "as-is" to be used at your own risk. |
---|
26 | ! |
---|
27 | ! |
---|
28 | |
---|
29 | ! |
---|
30 | ! ROOT1D.f90: find the "roots" of degree-k polynomials. |
---|
31 | ! |
---|
32 | ! Darren Engwirda |
---|
33 | ! 25-Mar-2019 |
---|
34 | ! de2363 [at] columbia [dot] edu |
---|
35 | ! |
---|
36 | ! |
---|
37 | |
---|
38 | pure subroutine roots_2(aa,bb,cc,xx,haveroot) |
---|
39 | |
---|
40 | ! |
---|
41 | ! solve:: aa * xx**2 + bb * xx**1 + cc = +0.0 . |
---|
42 | ! |
---|
43 | |
---|
44 | implicit none |
---|
45 | |
---|
46 | !------------------------------------------- arguments ! |
---|
47 | real*8 , intent( in) :: aa,bb,cc |
---|
48 | real*8 , intent(out) :: xx(1:2) |
---|
49 | logical, intent(out) :: haveroot |
---|
50 | |
---|
51 | !------------------------------------------- variables ! |
---|
52 | real*8 :: sq,ia,a0,b0,c0,x0 |
---|
53 | |
---|
54 | real*8, parameter :: rt = +1.d-14 |
---|
55 | |
---|
56 | a0 = abs(aa) |
---|
57 | b0 = abs(bb) |
---|
58 | c0 = abs(cc) |
---|
59 | |
---|
60 | sq = bb * bb - 4.0d+0 * aa * cc |
---|
61 | |
---|
62 | if (sq .ge. 0.0d+0) then |
---|
63 | |
---|
64 | sq = sqrt (sq) |
---|
65 | |
---|
66 | xx(1) = - bb + sq |
---|
67 | xx(2) = - bb - sq |
---|
68 | |
---|
69 | x0 = max(abs(xx(1)), & |
---|
70 | & abs(xx(2))) |
---|
71 | |
---|
72 | if (a0 .gt. (rt*x0)) then |
---|
73 | |
---|
74 | !-------------------------------------- degree-2 roots ! |
---|
75 | |
---|
76 | haveroot = .true. |
---|
77 | |
---|
78 | ia = 0.5d+0 / aa |
---|
79 | |
---|
80 | xx(1) = xx(1) * ia |
---|
81 | xx(2) = xx(2) * ia |
---|
82 | |
---|
83 | else & |
---|
84 | & if (b0 .gt. (rt*c0)) then |
---|
85 | |
---|
86 | !-------------------------------------- degree-1 roots ! |
---|
87 | |
---|
88 | haveroot = .true. |
---|
89 | |
---|
90 | xx(1) = - cc / bb |
---|
91 | xx(2) = - cc / bb |
---|
92 | |
---|
93 | else |
---|
94 | |
---|
95 | haveroot = .false. |
---|
96 | |
---|
97 | end if |
---|
98 | |
---|
99 | else |
---|
100 | |
---|
101 | haveroot = .false. |
---|
102 | |
---|
103 | end if |
---|
104 | |
---|
105 | return |
---|
106 | |
---|
107 | end subroutine |
---|
108 | |
---|
109 | |
---|
110 | |
---|