| 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.h90: 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 | |
|---|
