| 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 |
|
|
| 72 |
|
t = max( t, tlower ) |
| 73 |
|
t = min( t, tupper ) |
| 74 |
|
|
| 75 |
|
return |
| 76 |
|
end |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch