| 1 |
|
| 2 |
subroutine cubic( t, f, fp, ta, fa, fpa, tlower, tupper ) |
| 3 |
|
| 4 |
c----------------------------------------- |
| 5 |
c arguments |
| 6 |
c----------------------------------------- |
| 7 |
double precision t, f, fp, ta, fa, fpa, tlower, tupper |
| 8 |
|
| 9 |
c----------------------------------------- |
| 10 |
c local variables |
| 11 |
c----------------------------------------- |
| 12 |
double precision sign, den, anum |
| 13 |
double precision z1, b, discri |
| 14 |
|
| 15 |
c----------------------------------------- |
| 16 |
c Using f and fp at t and ta, |
| 17 |
c computes new t by cubic formula |
| 18 |
c safeguarded inside [tlower,tupper]. |
| 19 |
c----------------------------------------- |
| 20 |
z1 = dble(fp) + dble(fpa) - 3.d0*dble(fa-f)/dble(ta-t) |
| 21 |
b = z1 + dble(fp) |
| 22 |
|
| 23 |
c----------------------------------------- |
| 24 |
c first compute the discriminant |
| 25 |
c (without overflow) |
| 26 |
c----------------------------------------- |
| 27 |
if (abs(z1).le.1.) then |
| 28 |
discri = z1*z1-dble(fp)*dble(fpa) |
| 29 |
if (discri .lt. 0.d0) then |
| 30 |
if (fp.lt.0.) t = tupper |
| 31 |
if (fp.ge.0.) t = tlower |
| 32 |
return |
| 33 |
else |
| 34 |
discri = dsqrt(discri) |
| 35 |
end if |
| 36 |
else |
| 37 |
discri = dble(fp)/z1 |
| 38 |
discri = discri*dble(fpa) |
| 39 |
discri = z1-discri |
| 40 |
if (z1.ge.0.d0 .and. discri.ge.0.d0) then |
| 41 |
discri = dsqrt(z1)*dsqrt(discri) |
| 42 |
else if (z1.le.0.d0 .and. discri.le.0.d0) then |
| 43 |
discri = dsqrt(-z1)*dsqrt(-discri) |
| 44 |
else |
| 45 |
if (fp.lt.0.) t = tupper |
| 46 |
if (fp.ge.0.) t = tlower |
| 47 |
return |
| 48 |
end if |
| 49 |
end if |
| 50 |
|
| 51 |
c----------------------------------------- |
| 52 |
c discriminant nonnegative, |
| 53 |
c compute solution (without overflow) |
| 54 |
c----------------------------------------- |
| 55 |
if (t-ta .lt. 0.0) then |
| 56 |
discri = -discri |
| 57 |
end if |
| 58 |
|
| 59 |
sign = (t-ta)/abs(t-ta) |
| 60 |
if (sngl(b)*sign .gt. 0.0) then |
| 61 |
t = t + fp*(ta-t)/sngl(b+discri) |
| 62 |
else |
| 63 |
den = sngl(z1+b+dble(fpa)) |
| 64 |
anum = sngl(b-discri) |
| 65 |
if (abs((t-ta)*anum).lt.(tupper-tlower)*abs(den)) then |
| 66 |
t = t + anum*(ta-t)/den |
| 67 |
else |
| 68 |
t = tupper |
| 69 |
end if |
| 70 |
end if |
| 71 |
t = max( t, tlower ) |
| 72 |
t = min( t, tupper ) |
| 73 |
|
| 74 |
return |
| 75 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph