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 |