[13926] | 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 | ! |
---|
[14212] | 30 | ! ROOT1D.h90: find the "roots" of degree-k polynomials. |
---|
[13926] | 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 | |
---|
[14387] | 54 | real*8, parameter :: rt = +1.e-14 |
---|
[13926] | 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 | |
---|