V3FIT
akima.f
1 
2 
3 !ref: http://www.femto-st.fr/~daniau/fvn.html
4  SUBROUTINE fvn_akima(n,x,y,br,co)
5  USE precision
6  IMPLICIT NONE
7  INTEGER, INTENT(in) :: n
8  REAL(rprec), INTENT(in) :: x(n)
9  REAL(rprec), INTENT(in) :: y(n)
10  REAL(rprec), INTENT(out) :: br(n)
11  REAL(rprec), INTENT(out) :: co(4,n)
12 
13  REAL(rprec), ALLOCATABLE :: var(:),z(:)
14  REAL(rprec) :: wi_1,wi
15  INTEGER :: i
16  REAL(rprec) :: dx,a,b
17 
18  ! br is just a copy of x
19  br(:)=x(:)
20 
21  ALLOCATE(var(n+3))
22  ALLOCATE(z(n))
23  ! evaluate the variations
24  DO i=1, n-1
25  var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
26  ENDDO
27  var(n+2)=2.d0*var(n+1)-var(n)
28  var(n+3)=2.d0*var(n+2)-var(n+1)
29  var(2)=2.d0*var(3)-var(4)
30  var(1)=2.d0*var(2)-var(3)
31 
32  DO i = 1, n
33  wi_1=abs(var(i+3)-var(i+2))
34  wi=abs(var(i+1)-var(i))
35  IF ((wi_1+wi) .EQ. 0.d0) THEN
36  z(i)=(var(i+2)+var(i+1))/2.d0
37  ELSE
38  z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
39  ENDIF
40  ENDDO
41 
42  DO i=1, n-1
43  dx=x(i+1)-x(i)
44  a=(z(i+1)-z(i))*dx ! coeff INTermediaires pour calcul wd
45  b=y(i+1)-y(i)-z(i)*dx ! coeff INTermediaires pour calcul wd
46  co(1,i)=y(i)
47  co(2,i)=z(i)
48  co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx
49  co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2 !
50  ENDDO
51  co(1,n)=y(n)
52  co(2,n)=z(n)
53  co(3,n)=0.d0
54  co(4,n)=0.d0
55 
56  DEALLOCATE(z)
57  DEALLOCATE(var)
58 
59  END SUBROUTINE
60 
61  FUNCTION fvn_spline_eval(x,n,br,co)
62 ! VALUES
63  USE precision
64  IMPLICIT NONE
65  REAL(rprec) fvn_spline_eval
66  REAL(rprec), INTENT(in) :: x
67 ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
68  INTEGER, INTENT(in) :: n ! number of INTervals
69  REAL(rprec), INTENT(in) :: br(n+1) ! breakpoints (abcissas)
70  REAL(rprec), INTENT(in) :: co(4,n+1)! spline coeeficients
71  REAL(rprec) :: fvn_d_spline_eval
72 
73  INTEGER :: i
74  REAL(rprec) :: dx
75 
76 
77  IF (x<=br(1)) THEN
78  i=1
79  ELSE IF (x>=br(n+1)) THEN
80  i=n
81  ELSE
82  i=1
83  DO WHILE(x>=br(i))
84  i=i+1
85  ENDDO
86  i=i-1
87  ENDIF
88 
89  dx=x-br(i)
90  fvn_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
91 
92  END FUNCTION fvn_spline_eval
93 
94 
95  FUNCTION fvn_spline_devaldx(x,n,br,co)
96 ! 1ST DERIVATIVE
97  USE precision
98  IMPLICIT NONE
99  REAL(rprec) fvn_spline_devaldx
100  REAL(rprec), INTENT(in) :: x
101 ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
102  INTEGER, INTENT(in) :: n ! number of INTervals
103  REAL(rprec), INTENT(in) :: br(n+1) ! breakpoints
104  REAL(rprec), INTENT(in) :: co(4,n+1)! spline coeeficients
105  REAL(rprec) :: fvn_d_spline_eval
106 
107  INTEGER :: i
108  REAL(rprec) :: dx
109 
110 
111  IF (x<=br(1)) THEN
112  i=1
113  ELSE IF (x>=br(n+1)) THEN
114  i=n
115  ELSE
116  i=1
117  DO WHILE(x>=br(i))
118  i=i+1
119  ENDDO
120  i=i-1
121  ENDIF
122 
123  dx=x-br(i)
124  fvn_spline_devaldx=co(2,i)+2*co(3,i)*dx+3*co(4,i)*dx**2
125 
126  END FUNCTION fvn_spline_devaldx
127 
128  FUNCTION fvn_spline_d2evaldx2(x,n,br,co)
129 ! 2ND DERIVATIVE
130  USE precision
131  IMPLICIT NONE
132  REAL(rprec) fvn_spline_d2evaldx2
133  REAL(rprec), INTENT(in) :: x
134 ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
135  INTEGER, INTENT(in) :: n ! number of INTervals
136  REAL(rprec), INTENT(in) :: br(n+1) ! breakpoints
137  REAL(rprec), INTENT(in) :: co(4,n+1)! spline coeeficients
138 
139  INTEGER :: i
140  REAL(rprec) :: dx
141 
142 
143  IF (x<=br(1)) THEN
144  i=1
145  ELSE IF (x>=br(n+1)) THEN
146  i=n
147  ELSE
148  i=1
149  DO WHILE(x>=br(i))
150  i=i+1
151  ENDDO
152  i=i-1
153  ENDIF
154 
155  dx=x-br(i)
156  fvn_spline_d2evaldx2=2*co(3,i)+6*co(4,i)*dx
157 
158  END FUNCTION fvn_spline_d2evaldx2
159 
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17