@@ -10,27 +10,31 @@ program dbint4_test
10
10
11
11
implicit none
12
12
13
- integer (ip),parameter :: nx = 7 ! ! number of points in x
13
+ ! integer(ip),parameter :: nx = 7 !! number of points in x
14
14
integer (ip),parameter :: kx = 4 ! ! order in x
15
15
16
- logical ,parameter :: extrap = .true .
16
+ logical ,parameter :: extrap = .false .
17
17
real (wp),parameter :: rad2deg = 180.0_wp / acos (- 1.0_wp ) ! ! deg. to radians conversion factor
18
18
19
- real (wp) :: x(nx)
20
- real (wp),dimension (:),allocatable :: xval,fval ! (nx+2)
21
- real (wp) :: fcn_1d(nx)
22
- real (wp),dimension (:),allocatable :: tx ! (nx+2+kx)
23
- real (wp),dimension (:),allocatable :: bcoef ! (nx+2)
24
- real (wp) :: tol,val,tru,err,errmax,fbcr,fbcl
19
+ integer :: nx ! ! number of points in x
20
+ real (wp),dimension (:),allocatable :: x ! ! (nx)
21
+ real (wp),dimension (:),allocatable :: fcn_1d ! ! (nx)
22
+ real (wp),dimension (:,:),allocatable :: w ! ! (5,nx+2)
23
+ real (wp),dimension (:),allocatable :: xval,fval ! ! (nx+2)
24
+ real (wp),dimension (:),allocatable :: tx ! ! (nx+2+kx)
25
+ real (wp),dimension (:),allocatable :: bcoef ! ! (nx+2)
26
+ real (wp) :: tol,val,tru,err,errmax,tmp
25
27
logical :: fail
26
- integer (ip) :: n,i,idx,iflag,inbvx,ibcl,ibcr,kntopt
27
- real (wp),dimension (5 ,nx+2 ) :: w
28
+ integer (ip) :: i,idx,iflag,inbvx,ibcl,ibcr,kntopt
28
29
type (pyplot) :: plt
29
30
integer :: istat
30
31
integer :: icase
31
32
real (wp),dimension (3 ) :: tleft, tright
32
33
real (wp),dimension (3 * kx) :: w1_1d
33
34
35
+ real (wp),parameter :: fbcl = 1.0_wp ! ! left derivative
36
+ real (wp),parameter :: fbcr = 1.0_wp ! ! right derivative
37
+
34
38
character (len=* ),dimension (7 ),parameter :: labels = [' not-a-knot [db1ink] ' , &
35
39
' 2nd der=0, kntopt=1 [dbint4] ' , &
36
40
' 1st der=0, kntopt=1 [dbint4] ' , &
@@ -49,15 +53,19 @@ program dbint4_test
49
53
fail = .false.
50
54
tol = 100 * epsilon (1.0_wp )
51
55
idx = 0
52
- x = real ([1 ,2 ,3 ,4 ,5 ,6 ,7 ], wp) ! nx points
53
- fcn_1d = f1(x)
54
56
55
- ! initialize the plot:
56
- call plt% initialize(grid= .true. ,xlabel= ' x (deg)' ,ylabel= ' f(x)' ,&
57
- title= ' B-Spline End Conditions' ,legend= .true. )
58
- call plt% add_plot(x* rad2deg,fcn_1d,&
59
- label= ' Function $f(x) = \\sin(x)$' ,&
60
- linestyle= ' ko' ,markersize= 5 ,linewidth= 2 ,istat= istat)
57
+ ! create x and fcn_1d array: [0.1, ..., 7.1]
58
+ allocate (x(0 ), fcn_1d(0 ))
59
+ tmp = 0.0_wp
60
+ nx = 0
61
+ do
62
+ tmp = tmp + 0.1_wp
63
+ if (tmp>= 7.2_wp ) exit
64
+ nx = nx + 1
65
+ x = [x, tmp]
66
+ fcn_1d = [fcn_1d, f1(tmp)]
67
+ end do
68
+ allocate (w(5 ,nx+2 ) )
61
69
62
70
if (extrap) then
63
71
! points to evaluate [with extrapolation]:
@@ -68,6 +76,13 @@ program dbint4_test
68
76
end if
69
77
allocate (fval(size (xval)))
70
78
79
+ ! initialize the plot:
80
+ call plt% initialize(grid= .true. ,xlabel= ' x (deg)' ,ylabel= ' f(x)' ,&
81
+ title= ' B-Spline End Conditions' ,legend= .true. )
82
+ call plt% add_plot(x* rad2deg,fcn_1d,&
83
+ label= ' Function $f(x) = \\sin(x)$' ,&
84
+ linestyle= ' ko' ,markersize= 5 ,linewidth= 2 ,istat= istat)
85
+
71
86
do icase = 1 , 7
72
87
73
88
write (* ,* ) ' '
@@ -107,8 +122,6 @@ program dbint4_test
107
122
ibcl = 1
108
123
ibcr = 1
109
124
end if
110
- fbcl = 0.0_wp
111
- fbcr = 0.0_wp
112
125
113
126
write (* ,* ) ' kntopt: ' , kntopt
114
127
@@ -125,16 +138,16 @@ program dbint4_test
125
138
kntopt = 3
126
139
w = 0.0_wp
127
140
! WARNING: the knot values seem to make no difference in the result
128
- tleft = [- 999 .0_wp ,- 999 .0_wp ,- 999 .0_wp ]
129
- tright = [999 .0_wp , 999 .0_wp , 999 .0_wp ]
141
+ tleft = [0 .0_wp ,0 .0_wp ,0 .0_wp ]
142
+ tright = [8 .0_wp , 8 .0_wp , 8 .0_wp ]
130
143
call db1ink(x,nx,fcn_1d,kx,ibcl,ibcr,fbcl,fbcr,tleft,tright,tx,bcoef,iflag)
131
144
end select
132
145
133
- write (* ,* ) ' '
134
- write (* ,* ) ' x: ' , x
135
- write (* ,* ) ' '
136
- write (* ,* ) ' tx: ' , tx
137
- write (* ,* ) ' '
146
+ ! write(*,*) ''
147
+ ! write(*,*) 'x: ', x
148
+ ! write(*,*) ''
149
+ ! write(*,*) 'tx: ', tx
150
+ ! write(*,*) ''
138
151
139
152
end if
140
153
if (iflag/= 0 ) then
0 commit comments