29
30
31
32
33
34
35#include "implicit_f.inc"
36
37
38
39
40
41
42
43
44
45
46
52
53
54
55
56 my_real :: dk21,dk32,dk43,dk31,dk42,kt1,kt2,kt3,kt4
58 my_real :: a1x,a1y,a2x,a2y,a3x,a3y,b1x,b1y,b2x,b2y
59 my_real :: cx,cy,c_dx,c_dy,c_ddx,c_ddy
60 my_real :: a1px,a1py,a2px,a2py,a3px,a3py,b1px,b1py,b2px,b2py
61 my_real :: b1ppx,b1ppy,b2ppx,b2ppy
62
63
64
65 dk21 = knots(2) - knots(1)
66 dk32 = knots(3) - knots(2)
67 dk43 = knots(4) - knots(3)
68 dk31 = knots(3) - knots(1)
69 dk42 = knots(4) - knots(2)
70 tt = knots(2) + t * dk32
71 kt1 = knots(1) - tt
72 kt2 = knots(2) - tt
73 kt3 = knots(3) - tt
74 kt4 = knots(4) - tt
75
76 a1x = (kt2*ptx(1) - kt1*ptx(2)) / dk21
77 a1y = (kt2*pty(1) - kt1*pty(2)) / dk21
78 a2x = (kt3*ptx(2) - kt2*ptx(3)) / dk32
79 a2y = (kt3*pty(2) - kt2*pty(3)) / dk32
80 a3x = (kt4*ptx(3) - kt3*ptx(4)) / dk43
81 a3y = (kt4*pty(3) - kt3*pty(4)) / dk43
82
83 b1x = (kt3*a1x - kt1*a2x) / dk31
84 b1y = (kt3*a1y - kt1*a2y) / dk31
85 b2x = (kt4*a2x - kt2*a3x) / dk42
86 b2y = (kt4*a2y - kt2*a3y) / dk42
87
88 nx = (kt3*b1x - kt2*b2x) / dk32
89 ny = (kt3*b1y - kt2*b2y) / dk32
90
91 a1px = (ptx(2) - ptx(1)) / dk21
92 a1py = (pty(2) - pty(1)) / dk21
93 a2px = (ptx(3) - ptx(2)) / dk32
94 a2py = (pty(3) - pty(2)) / dk32
95 a3px = (ptx(4) - ptx(3)) / dk43
96 a3py = (pty(4) - pty(3)) / dk43
97 b1px = (a2x-a1x)/ dk31 + kt3 / dk31*a1px - kt1 / dk31*a2px
98 b1py = (a2y-a1y)/ dk31 + kt3 / dk31*a1py - kt1 / dk31*a2py
99 b2px = (a3x-a2x)/ dk42 + kt4 / dk42*a2px - kt2 / dk42*a3px
100 b2py = (a3y-a2y)/ dk42 + kt4 / dk42*a2py - kt2 / dk42*a3py
101
102
103 b1ppx= two*(a2px - a1px) / dk31
104 b1ppy= two*(a2py - a1py) / dk31
105 b2ppx= two*(a3px - a2px) / dk42
106 b2ppy= two*(a3py - two*a2py) / dk42
107
108
109
110 RETURN