lxsd / fvn

Download as
Browse Code »

Commit e711bb807cee3a51f1eae75b9041469bcb96f2af

Authored by cwaterkeyn
1 parent 41811905dd
Exists in master and in 2 other branches geevx_2020, multi

ChW 02/2010 for typing errors between sp_kind and ip_kind 

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@68 b657c933-2333-4658-acf2-d3c7c2708721

Showing 1 changed file with 20 additions and 19 deletions Inline Diff

fvn_interpol/fvn_interpol.f90
Diff comments   View file @ e711bb8
module fvn_interpol 1 1 module fvn_interpol
use fvn_common 2 2 use fvn_common
implicit none 3 3 implicit none
4 4
5 5
! Utility procedure find interval 6 6 ! Utility procedure find interval
interface fvn_find_interval 7 7 interface fvn_find_interval
module procedure fvn_s_find_interval,fvn_d_find_interval 8 8 module procedure fvn_s_find_interval,fvn_d_find_interval
end interface fvn_find_interval 9 9 end interface fvn_find_interval
10 10
! Quadratic 1D interpolation 11 11 ! Quadratic 1D interpolation
interface fvn_quad_interpol 12 12 interface fvn_quad_interpol
module procedure fvn_s_quad_interpol,fvn_d_quad_interpol 13 13 module procedure fvn_s_quad_interpol,fvn_d_quad_interpol
end interface fvn_quad_interpol 14 14 end interface fvn_quad_interpol
15 15
! Quadratic 2D interpolation 16 16 ! Quadratic 2D interpolation
interface fvn_quad_2d_interpol 17 17 interface fvn_quad_2d_interpol
module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol 18 18 module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol
end interface fvn_quad_2d_interpol 19 19 end interface fvn_quad_2d_interpol
20 20
! Quadratic 3D interpolation 21 21 ! Quadratic 3D interpolation
interface fvn_quad_3d_interpol 22 22 interface fvn_quad_3d_interpol
module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol 23 23 module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol
end interface fvn_quad_3d_interpol 24 24 end interface fvn_quad_3d_interpol
25 25
! Akima interpolation 26 26 ! Akima interpolation
interface fvn_akima 27 27 interface fvn_akima
module procedure fvn_s_akima,fvn_d_akima 28 28 module procedure fvn_s_akima,fvn_d_akima
end interface fvn_akima 29 29 end interface fvn_akima
30 30
! Akima evaluation 31 31 ! Akima evaluation
interface fvn_spline_eval 32 32 interface fvn_spline_eval
module procedure fvn_s_spline_eval,fvn_d_spline_eval 33 33 module procedure fvn_s_spline_eval,fvn_d_spline_eval
end interface fvn_spline_eval 34 34 end interface fvn_spline_eval
35 35
contains 36 36 contains
37 37
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 38 38 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 39 39 !
! Quadratic interpolation of tabulated function of 1,2 or 3 variables 40 40 ! Quadratic interpolation of tabulated function of 1,2 or 3 variables
! 41 41 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 42 42 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43 43
subroutine fvn_s_find_interval(x,i,xdata,n) 44 44 subroutine fvn_s_find_interval(x,i,xdata,n)
implicit none 45 45 implicit none
! This routine find the indice i where xdata(i) <= x < xdata(i+1) 46 46 ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
! xdata(n) must contains a set of increasingly ordered values 47 47 ! xdata(n) must contains a set of increasingly ordered values
! if x < xdata(1) i=0 is returned 48 48 ! if x < xdata(1) i=0 is returned
! if x > xdata(n) i=n is returned 49 49 ! if x > xdata(n) i=n is returned
! special case is where x=xdata(n) then n-1 is returned so 50 50 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 51 51 ! we will not exclude the upper limit
! a simple dichotomy method is used 52 52 ! a simple dichotomy method is used
53 53
real(kind=sp_kind), intent(in) :: x 54 54 real(kind=sp_kind), intent(in) :: x
integer(kind=sp_kind), intent(in) :: n 55 55 integer(kind=ip_kind), intent(in) :: n
real(kind=sp_kind), intent(in), dimension(n) :: xdata 56 56 real(kind=sp_kind), intent(in), dimension(n) :: xdata
integer(kind=sp_kind), intent(out) :: i 57 57 integer(kind=ip_kind), intent(out) :: i
58 58
integer(kind=sp_kind) :: imin,imax,imoyen 59 59 integer(kind=ip_kind) :: imin,imax,imoyen
60 60
! special case is where x=xdata(n) then n-1 is returned so 61 61 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 62 62 ! we will not exclude the upper limit
if (x == xdata(n)) then 63 63 if (x == xdata(n)) then
i=n-1 64 64 i=n-1
return 65 65 return
end if 66 66 end if
67 67
! if x < xdata(1) i=0 is returned 68 68 ! if x < xdata(1) i=0 is returned
if (x < xdata(1)) then 69 69 if (x < xdata(1)) then
i=0 70 70 i=0
return 71 71 return
end if 72 72 end if
73 73
! if x > xdata(n) i=n is returned 74 74 ! if x > xdata(n) i=n is returned
if (x > xdata(n)) then 75 75 if (x > xdata(n)) then
i=n 76 76 i=n
return 77 77 return
end if 78 78 end if
79 79
! here xdata(1) <= x <= xdata(n) 80 80 ! here xdata(1) <= x <= xdata(n)
imin=0 81 81 imin=0
imax=n+1 82 82 imax=n+1
83 83
do while((imax-imin) > 1) 84 84 do while((imax-imin) > 1)
imoyen=(imax+imin)/2 85 85 imoyen=(imax+imin)/2
if (x >= xdata(imoyen)) then 86 86 if (x >= xdata(imoyen)) then
imin=imoyen 87 87 imin=imoyen
else 88 88 else
imax=imoyen 89 89 imax=imoyen
end if 90 90 end if
end do 91 91 end do
92 92
i=imin 93 93 i=imin
94 94
end subroutine 95 95 end subroutine
96 96
97 97
subroutine fvn_d_find_interval(x,i,xdata,n) 98 98 subroutine fvn_d_find_interval(x,i,xdata,n)
implicit none 99 99 implicit none
! This routine find the indice i where xdata(i) <= x < xdata(i+1) 100 100 ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
! xdata(n) must contains a set of increasingly ordered values 101 101 ! xdata(n) must contains a set of increasingly ordered values
! if x < xdata(1) i=0 is returned 102 102 ! if x < xdata(1) i=0 is returned
! if x > xdata(n) i=n is returned 103 103 ! if x > xdata(n) i=n is returned
! special case is where x=xdata(n) then n-1 is returned so 104 104 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 105 105 ! we will not exclude the upper limit
! a simple dichotomy method is used 106 106 ! a simple dichotomy method is used
107 107
real(kind=dp_kind), intent(in) :: x 108 108 real(kind=dp_kind), intent(in) :: x
integer(kind=sp_kind), intent(in) :: n 109 109 integer(kind=ip_kind), intent(in) :: n
real(kind=dp_kind), intent(in), dimension(n) :: xdata 110 110 real(kind=dp_kind), intent(in), dimension(n) :: xdata
integer(kind=sp_kind), intent(out) :: i 111 111 integer(kind=ip_kind), intent(out) :: i
112 112
integer(kind=sp_kind) :: imin,imax,imoyen 113 113 integer(kind=ip_kind) :: imin,imax,imoyen
114 114
! special case is where x=xdata(n) then n-1 is returned so 115 115 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 116 116 ! we will not exclude the upper limit
if (x == xdata(n)) then 117 117 if (x == xdata(n)) then
i=n-1 118 118 i=n-1
return 119 119 return
end if 120 120 end if
121 121
! if x < xdata(1) i=0 is returned 122 122 ! if x < xdata(1) i=0 is returned
if (x < xdata(1)) then 123 123 if (x < xdata(1)) then
i=0 124 124 i=0
return 125 125 return
end if 126 126 end if
127 127
! if x > xdata(n) i=n is returned 128 128 ! if x > xdata(n) i=n is returned
if (x > xdata(n)) then 129 129 if (x > xdata(n)) then
i=n 130 130 i=n
return 131 131 return
end if 132 132 end if
133 133
! here xdata(1) <= x <= xdata(n) 134 134 ! here xdata(1) <= x <= xdata(n)
imin=0 135 135 imin=0
imax=n+1 136 136 imax=n+1
137 137
do while((imax-imin) > 1) 138 138 do while((imax-imin) > 1)
imoyen=(imax+imin)/2 139 139 imoyen=(imax+imin)/2
if (x >= xdata(imoyen)) then 140 140 if (x >= xdata(imoyen)) then
imin=imoyen 141 141 imin=imoyen
else 142 142 else
imax=imoyen 143 143 imax=imoyen
end if 144 144 end if
end do 145 145 end do
146 146
i=imin 147 147 i=imin
148 148
end subroutine 149 149 end subroutine
150 150
151 151
function fvn_s_quad_interpol(x,n,xdata,ydata) 152 152 function fvn_s_quad_interpol(x,n,xdata,ydata)
implicit none 153 153 implicit none
! This function evaluate the value of a function defined by a set of points 154 154 ! This function evaluate the value of a function defined by a set of points
! and values, using a quadratic interpolation 155 155 ! and values, using a quadratic interpolation
! xdata must be increasingly ordered 156 156 ! xdata must be increasingly ordered
! x must be within xdata(1) and xdata(n) to actually do interpolation 157 157 ! x must be within xdata(1) and xdata(n) to actually do interpolation
! otherwise extrapolation is done 158 158 ! otherwise extrapolation is done
integer(kind=sp_kind), intent(in) :: n 159 159 integer(kind=ip_kind), intent(in) :: n
real(kind=sp_kind), intent(in), dimension(n) :: xdata,ydata 160 160 real(kind=sp_kind), intent(in), dimension(n) :: xdata,ydata
real(kind=sp_kind), intent(in) :: x 161 161 real(kind=sp_kind), intent(in) :: x
real(kind=sp_kind) :: fvn_s_quad_interpol 162 162 real(kind=sp_kind) :: fvn_s_quad_interpol
163 163
integer(kind=sp_kind) :: iinf,base,i,j 164 164 integer(kind=ip_kind) :: iinf,base,i,j
real(kind=sp_kind) :: p 165 165 real(kind=sp_kind) :: p
166 166
call fvn_s_find_interval(x,iinf,xdata,n) 167 167 call fvn_s_find_interval(x,iinf,xdata,n)
168 168
! Settings for extrapolation 169 169 ! Settings for extrapolation
if (iinf==0) then 170 170 if (iinf==0) then
! TODO -> Lower bound extrapolation warning 171 171 ! TODO -> Lower bound extrapolation warning
iinf=1 172 172 iinf=1
end if 173 173 end if
174 174
if (iinf==n) then 175 175 if (iinf==n) then
! TODO -> Higher bound extrapolation warning 176 176 ! TODO -> Higher bound extrapolation warning
iinf=n-1 177 177 iinf=n-1
end if 178 178 end if
179 179
! The three points we will use are iinf-1,iinf and iinf+1 with the 180 180 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 181 181 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (iinf==1) then 182 182 if (iinf==1) then
base=0 183 183 base=0
else 184 184 else
base=iinf-2 185 185 base=iinf-2
end if 186 186 end if
187 187
! The three points we will use are : 188 188 ! The three points we will use are :
! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3) 189 189 ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
190 190
! Straight forward Lagrange polynomial 191 191 ! Straight forward Lagrange polynomial
fvn_s_quad_interpol=0. 192 192 fvn_s_quad_interpol=0.
do i=1,3 193 193 do i=1,3
! polynome i 194 194 ! polynome i
p=ydata(base+i) 195 195 p=ydata(base+i)
do j=1,3 196 196 do j=1,3
if (j /= i) then 197 197 if (j /= i) then
p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j)) 198 198 p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
end if 199 199 end if
end do 200 200 end do
fvn_s_quad_interpol=fvn_s_quad_interpol+p 201 201 fvn_s_quad_interpol=fvn_s_quad_interpol+p
end do 202 202 end do
203 203
end function 204 204 end function
205 205
206 206
function fvn_d_quad_interpol(x,n,xdata,ydata) 207 207 function fvn_d_quad_interpol(x,n,xdata,ydata)
implicit none 208 208 implicit none
! This function evaluate the value of a function defined by a set of points 209 209 ! This function evaluate the value of a function defined by a set of points
! and values, using a quadratic interpolation 210 210 ! and values, using a quadratic interpolation
! xdata must be increasingly ordered 211 211 ! xdata must be increasingly ordered
! x must be within xdata(1) and xdata(n) to actually do interpolation 212 212 ! x must be within xdata(1) and xdata(n) to actually do interpolation
! otherwise extrapolation is done 213 213 ! otherwise extrapolation is done
integer(kind=sp_kind), intent(in) :: n 214 214 integer(kind=ip_kind), intent(in) :: n
real(kind=dp_kind), intent(in), dimension(n) :: xdata,ydata 215 215 real(kind=dp_kind), intent(in), dimension(n) :: xdata,ydata
real(kind=dp_kind), intent(in) :: x 216 216 real(kind=dp_kind), intent(in) :: x
real(kind=dp_kind) :: fvn_d_quad_interpol 217 217 real(kind=dp_kind) :: fvn_d_quad_interpol
218 218
integer(kind=sp_kind) :: iinf,base,i,j 219 219 integer(kind=ip_kind) :: iinf,base,i,j
real(kind=dp_kind) :: p 220 220 real(kind=dp_kind) :: p
221 221
call fvn_d_find_interval(x,iinf,xdata,n) 222 222 call fvn_d_find_interval(x,iinf,xdata,n)
223 223
! Settings for extrapolation 224 224 ! Settings for extrapolation
if (iinf==0) then 225 225 if (iinf==0) then
! TODO -> Lower bound extrapolation warning 226 226 ! TODO -> Lower bound extrapolation warning
iinf=1 227 227 iinf=1
end if 228 228 end if
229 229
if (iinf==n) then 230 230 if (iinf==n) then
! TODO Higher bound extrapolation warning 231 231 ! TODO Higher bound extrapolation warning
iinf=n-1 232 232 iinf=n-1
end if 233 233 end if
234 234
! The three points we will use are iinf-1,iinf and iinf+1 with the 235 235 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 236 236 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (iinf==1) then 237 237 if (iinf==1) then
base=0 238 238 base=0
else 239 239 else
base=iinf-2 240 240 base=iinf-2
end if 241 241 end if
242 242
! The three points we will use are : 243 243 ! The three points we will use are :
! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3) 244 244 ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
245 245
! Straight forward Lagrange polynomial 246 246 ! Straight forward Lagrange polynomial
fvn_d_quad_interpol=0. 247 247 fvn_d_quad_interpol=0.
do i=1,3 248 248 do i=1,3
! polynome i 249 249 ! polynome i
p=ydata(base+i) 250 250 p=ydata(base+i)
do j=1,3 251 251 do j=1,3
if (j /= i) then 252 252 if (j /= i) then
p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j)) 253 253 p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
end if 254 254 end if
end do 255 255 end do
fvn_d_quad_interpol=fvn_d_quad_interpol+p 256 256 fvn_d_quad_interpol=fvn_d_quad_interpol+p
end do 257 257 end do
258 258
end function 259 259 end function
260 260
261 261
function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata) 262 262 function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
implicit none 263 263 implicit none
! This function evaluate the value of a two variable function defined by a 264 264 ! This function evaluate the value of a two variable function defined by a
! set of points and values, using a quadratic interpolation 265 265 ! set of points and values, using a quadratic interpolation
! xdata and ydata must be increasingly ordered 266 266 ! xdata and ydata must be increasingly ordered
! the couple (x,y) must be as x within xdata(1) and xdata(nx) and 267 267 ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
! y within ydata(1) and ydata(ny) to actually do interpolation 268 268 ! y within ydata(1) and ydata(ny) to actually do interpolation
! otherwise extrapolation is done 269 269 ! otherwise extrapolation is done
integer(kind=sp_kind), intent(in) :: nx,ny 270 270 integer(kind=ip_kind), intent(in) :: nx,ny
real(kind=sp_kind), intent(in) :: x,y 271 271 real(kind=sp_kind), intent(in) :: x,y
real(kind=sp_kind), intent(in), dimension(nx) :: xdata 272 272 real(kind=sp_kind), intent(in), dimension(nx) :: xdata
real(kind=sp_kind), intent(in), dimension(ny) :: ydata 273 273 real(kind=sp_kind), intent(in), dimension(ny) :: ydata
real(kind=sp_kind), intent(in), dimension(nx,ny) :: zdata 274 274 real(kind=sp_kind), intent(in), dimension(nx,ny) :: zdata
real(kind=sp_kind) :: fvn_s_quad_2d_interpol 275 275 real(kind=sp_kind) :: fvn_s_quad_2d_interpol
276 276
integer(kind=sp_kind) :: ixinf,iyinf,basex,basey,i 277 277 integer(kind=ip_kind) :: ixinf,iyinf,basex,basey,i
real(kind=sp_kind),dimension(3) :: ztmp 278 278 real(kind=sp_kind),dimension(3) :: ztmp
!real(kind=4), external :: fvn_s_quad_interpol 279 279 !real(kind=4), external :: fvn_s_quad_interpol
280 280
call fvn_s_find_interval(x,ixinf,xdata,nx) 281 281 call fvn_s_find_interval(x,ixinf,xdata,nx)
call fvn_s_find_interval(y,iyinf,ydata,ny) 282 282 call fvn_s_find_interval(y,iyinf,ydata,ny)
283 283
! Settings for extrapolation 284 284 ! Settings for extrapolation
if (ixinf==0) then 285 285 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 286 286 ! TODO -> Lower x bound extrapolation warning
ixinf=1 287 287 ixinf=1
end if 288 288 end if
289 289
if (ixinf==nx) then 290 290 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 291 291 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 292 292 ixinf=nx-1
end if 293 293 end if
294 294
if (iyinf==0) then 295 295 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 296 296 ! TODO -> Lower y bound extrapolation warning
iyinf=1 297 297 iyinf=1
end if 298 298 end if
299 299
if (iyinf==ny) then 300 300 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 301 301 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 302 302 iyinf=ny-1
end if 303 303 end if
304 304
! The three points we will use are iinf-1,iinf and iinf+1 with the 305 305 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 306 306 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 307 307 if (ixinf==1) then
basex=0 308 308 basex=0
else 309 309 else
basex=ixinf-2 310 310 basex=ixinf-2
end if 311 311 end if
312 312
if (iyinf==1) then 313 313 if (iyinf==1) then
basey=0 314 314 basey=0
else 315 315 else
basey=iyinf-2 316 316 basey=iyinf-2
end if 317 317 end if
318 318
! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3) 319 319 ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
! stored in ztmp(1:3) 320 320 ! stored in ztmp(1:3)
do i=1,3 321 321 do i=1,3
ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i)) 322 322 ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i))
end do 323 323 end do
324 324
! Then we make an interpolation for y using previous interpolations 325 325 ! Then we make an interpolation for y using previous interpolations
fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp) 326 326 fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
end function 327 327 end function
328 328
329 329
function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata) 330 330 function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
implicit none 331 331 implicit none
! This function evaluate the value of a two variable function defined by a 332 332 ! This function evaluate the value of a two variable function defined by a
! set of points and values, using a quadratic interpolation 333 333 ! set of points and values, using a quadratic interpolation
! xdata and ydata must be increasingly ordered 334 334 ! xdata and ydata must be increasingly ordered
! the couple (x,y) must be as x within xdata(1) and xdata(nx) and 335 335 ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
! y within ydata(1) and ydata(ny) to actually do interpolation 336 336 ! y within ydata(1) and ydata(ny) to actually do interpolation
! otherwise extrapolation is done 337 337 ! otherwise extrapolation is done
integer(kind=sp_kind), intent(in) :: nx,ny 338 338 integer(kind=ip_kind), intent(in) :: nx,ny
real(kind=dp_kind), intent(in) :: x,y 339 339 real(kind=dp_kind), intent(in) :: x,y
real(kind=dp_kind), intent(in), dimension(nx) :: xdata 340 340 real(kind=dp_kind), intent(in), dimension(nx) :: xdata
real(kind=dp_kind), intent(in), dimension(ny) :: ydata 341 341 real(kind=dp_kind), intent(in), dimension(ny) :: ydata
real(kind=dp_kind), intent(in), dimension(nx,ny) :: zdata 342 342 real(kind=dp_kind), intent(in), dimension(nx,ny) :: zdata
real(kind=dp_kind) :: fvn_d_quad_2d_interpol 343 343 real(kind=dp_kind) :: fvn_d_quad_2d_interpol
344 344
integer(kind=sp_kind) :: ixinf,iyinf,basex,basey,i 345 345 integer(kind=ip_kind) :: ixinf,iyinf,basex,basey,i
real(kind=dp_kind),dimension(3) :: ztmp 346 346 real(kind=dp_kind),dimension(3) :: ztmp
!real(kind=8), external :: fvn_d_quad_interpol 347 347 !real(kind=8), external :: fvn_d_quad_interpol
348 348
call fvn_d_find_interval(x,ixinf,xdata,nx) 349 349 call fvn_d_find_interval(x,ixinf,xdata,nx)
call fvn_d_find_interval(y,iyinf,ydata,ny) 350 350 call fvn_d_find_interval(y,iyinf,ydata,ny)
351 351
! Settings for extrapolation 352 352 ! Settings for extrapolation
if (ixinf==0) then 353 353 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 354 354 ! TODO -> Lower x bound extrapolation warning
ixinf=1 355 355 ixinf=1
end if 356 356 end if
357 357
if (ixinf==nx) then 358 358 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 359 359 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 360 360 ixinf=nx-1
end if 361 361 end if
362 362
if (iyinf==0) then 363 363 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 364 364 ! TODO -> Lower y bound extrapolation warning
iyinf=1 365 365 iyinf=1
end if 366 366 end if
367 367
if (iyinf==ny) then 368 368 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 369 369 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 370 370 iyinf=ny-1
end if 371 371 end if
372 372
! The three points we will use are iinf-1,iinf and iinf+1 with the 373 373 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 374 374 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 375 375 if (ixinf==1) then
basex=0 376 376 basex=0
else 377 377 else
basex=ixinf-2 378 378 basex=ixinf-2
end if 379 379 end if
380 380
if (iyinf==1) then 381 381 if (iyinf==1) then
basey=0 382 382 basey=0
else 383 383 else
basey=iyinf-2 384 384 basey=iyinf-2
end if 385 385 end if
386 386
! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3) 387 387 ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
! stored in ztmp(1:3) 388 388 ! stored in ztmp(1:3)
do i=1,3 389 389 do i=1,3
ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i)) 390 390 ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i))
end do 391 391 end do
392 392
! Then we make an interpolation for y using previous interpolations 393 393 ! Then we make an interpolation for y using previous interpolations
fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp) 394 394 fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
end function 395 395 end function
396 396
397 397
function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata) 398 398 function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
implicit none 399 399 implicit none
! This function evaluate the value of a 3 variables function defined by a 400 400 ! This function evaluate the value of a 3 variables function defined by a
! set of points and values, using a quadratic interpolation 401 401 ! set of points and values, using a quadratic interpolation
! xdata, ydata and zdata must be increasingly ordered 402 402 ! xdata, ydata and zdata must be increasingly ordered
! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 403 403 ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually
! perform an interpolation, otherwise extrapolation is done 404 404 ! perform an interpolation, otherwise extrapolation is done
integer(kind=sp_kind), intent(in) :: nx,ny,nz 405 405 integer(kind=ip_kind), intent(in) :: nx,ny,nz
real(kind=sp_kind), intent(in) :: x,y,z 406 406 real(kind=sp_kind), intent(in) :: x,y,z
real(kind=sp_kind), intent(in), dimension(nx) :: xdata 407 407 real(kind=sp_kind), intent(in), dimension(nx) :: xdata
real(kind=sp_kind), intent(in), dimension(ny) :: ydata 408 408 real(kind=sp_kind), intent(in), dimension(ny) :: ydata
real(kind=sp_kind), intent(in), dimension(nz) :: zdata 409 409 real(kind=sp_kind), intent(in), dimension(nz) :: zdata
real(kind=sp_kind), intent(in), dimension(nx,ny,nz) :: tdata 410 410 real(kind=sp_kind), intent(in), dimension(nx,ny,nz) :: tdata
real(kind=sp_kind) :: fvn_s_quad_3d_interpol 411 411 real(kind=sp_kind) :: fvn_s_quad_3d_interpol
412 412
integer(kind=sp_kind) :: ixinf,iyinf,izinf,basex,basey,basez,i,j 413 413 integer(kind=ip_kind) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
!real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol 414 414 !real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol
real(kind=sp_kind),dimension(3,3) :: ttmp 415 415 real(kind=sp_kind),dimension(3,3) :: ttmp
416 416
call fvn_s_find_interval(x,ixinf,xdata,nx) 417 417 call fvn_s_find_interval(x,ixinf,xdata,nx)
call fvn_s_find_interval(y,iyinf,ydata,ny) 418 418 call fvn_s_find_interval(y,iyinf,ydata,ny)
call fvn_s_find_interval(z,izinf,zdata,nz) 419 419 call fvn_s_find_interval(z,izinf,zdata,nz)
420 420
! Settings for extrapolation 421 421 ! Settings for extrapolation
if (ixinf==0) then 422 422 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 423 423 ! TODO -> Lower x bound extrapolation warning
ixinf=1 424 424 ixinf=1
end if 425 425 end if
426 426
if (ixinf==nx) then 427 427 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 428 428 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 429 429 ixinf=nx-1
end if 430 430 end if
431 431
if (iyinf==0) then 432 432 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 433 433 ! TODO -> Lower y bound extrapolation warning
iyinf=1 434 434 iyinf=1
end if 435 435 end if
436 436
if (iyinf==ny) then 437 437 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 438 438 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 439 439 iyinf=ny-1
end if 440 440 end if
441 441
if (izinf==0) then 442 442 if (izinf==0) then
! TODO -> Lower z bound extrapolation warning 443 443 ! TODO -> Lower z bound extrapolation warning
izinf=1 444 444 izinf=1
end if 445 445 end if
446 446
if (izinf==nz) then 447 447 if (izinf==nz) then
! TODO -> Higher z bound extrapolation warning 448 448 ! TODO -> Higher z bound extrapolation warning
izinf=nz-1 449 449 izinf=nz-1
end if 450 450 end if
451 451
! The three points we will use are iinf-1,iinf and iinf+1 with the 452 452 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 453 453 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 454 454 if (ixinf==1) then
basex=0 455 455 basex=0
else 456 456 else
basex=ixinf-2 457 457 basex=ixinf-2
end if 458 458 end if
459 459
if (iyinf==1) then 460 460 if (iyinf==1) then
basey=0 461 461 basey=0
else 462 462 else
basey=iyinf-2 463 463 basey=iyinf-2
end if 464 464 end if
465 465
if (izinf==1) then 466 466 if (izinf==1) then
basez=0 467 467 basez=0
else 468 468 else
basez=izinf-2 469 469 basez=izinf-2
end if 470 470 end if
471 471
! We first make 9 one dimensional interpolation on variable x. 472 472 ! We first make 9 one dimensional interpolation on variable x.
! results are stored in ttmp 473 473 ! results are stored in ttmp
do i=1,3 474 474 do i=1,3
do j=1,3 475 475 do j=1,3
ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j)) 476 476 ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
end do 477 477 end do
end do 478 478 end do
479 479
! We then make a 2 dimensionnal interpolation on variables y and z 480 480 ! We then make a 2 dimensionnal interpolation on variables y and z
fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, & 481 481 fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &
3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp) 482 482 3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
end function 483 483 end function
484 484
485 485
function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata) 486 486 function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
implicit none 487 487 implicit none
! This function evaluate the value of a 3 variables function defined by a 488 488 ! This function evaluate the value of a 3 variables function defined by a
! set of points and values, using a quadratic interpolation 489 489 ! set of points and values, using a quadratic interpolation
! xdata, ydata and zdata must be increasingly ordered 490 490 ! xdata, ydata and zdata must be increasingly ordered
! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 491 491 ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually
! perform an interpolation, otherwise extrapolation is done 492 492 ! perform an interpolation, otherwise extrapolation is done
integer(kind=sp_kind), intent(in) :: nx,ny,nz 493 493 integer(kind=ip_kind), intent(in) :: nx,ny,nz
real(kind=dp_kind), intent(in) :: x,y,z 494 494 real(kind=dp_kind), intent(in) :: x,y,z
real(kind=dp_kind), intent(in), dimension(nx) :: xdata 495 495 real(kind=dp_kind), intent(in), dimension(nx) :: xdata
real(kind=dp_kind), intent(in), dimension(ny) :: ydata 496 496 real(kind=dp_kind), intent(in), dimension(ny) :: ydata
real(kind=dp_kind), intent(in), dimension(nz) :: zdata 497 497 real(kind=dp_kind), intent(in), dimension(nz) :: zdata
real(kind=dp_kind), intent(in), dimension(nx,ny,nz) :: tdata 498 498 real(kind=dp_kind), intent(in), dimension(nx,ny,nz) :: tdata
real(kind=dp_kind) :: fvn_d_quad_3d_interpol 499 499 real(kind=dp_kind) :: fvn_d_quad_3d_interpol
500 500
integer(kind=sp_kind) :: ixinf,iyinf,izinf,basex,basey,basez,i,j 501 501 integer(kind=ip_kind) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
!real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol 502 502 !real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol
real(kind=dp_kind),dimension(3,3) :: ttmp 503 503 real(kind=dp_kind),dimension(3,3) :: ttmp
504 504
call fvn_d_find_interval(x,ixinf,xdata,nx) 505 505 call fvn_d_find_interval(x,ixinf,xdata,nx)
call fvn_d_find_interval(y,iyinf,ydata,ny) 506 506 call fvn_d_find_interval(y,iyinf,ydata,ny)
call fvn_d_find_interval(z,izinf,zdata,nz) 507 507 call fvn_d_find_interval(z,izinf,zdata,nz)
508 508
! Settings for extrapolation 509 509 ! Settings for extrapolation
if (ixinf==0) then 510 510 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 511 511 ! TODO -> Lower x bound extrapolation warning
ixinf=1 512 512 ixinf=1
end if 513 513 end if
514 514
if (ixinf==nx) then 515 515 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 516 516 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 517 517 ixinf=nx-1
end if 518 518 end if
519 519
if (iyinf==0) then 520 520 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 521 521 ! TODO -> Lower y bound extrapolation warning
iyinf=1 522 522 iyinf=1
end if 523 523 end if
524 524
if (iyinf==ny) then 525 525 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 526 526 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 527 527 iyinf=ny-1
end if 528 528 end if
529 529
if (izinf==0) then 530 530 if (izinf==0) then
! TODO -> Lower z bound extrapolation warning 531 531 ! TODO -> Lower z bound extrapolation warning
izinf=1 532 532 izinf=1
end if 533 533 end if
534 534
if (izinf==nz) then 535 535 if (izinf==nz) then
! TODO -> Higher z bound extrapolation warning 536 536 ! TODO -> Higher z bound extrapolation warning
izinf=nz-1 537 537 izinf=nz-1
end if 538 538 end if
539 539
! The three points we will use are iinf-1,iinf and iinf+1 with the 540 540 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 541 541 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 542 542 if (ixinf==1) then
basex=0 543 543 basex=0
else 544 544 else
basex=ixinf-2 545 545 basex=ixinf-2
end if 546 546 end if
547 547
if (iyinf==1) then 548 548 if (iyinf==1) then
basey=0 549 549 basey=0
else 550 550 else
basey=iyinf-2 551 551 basey=iyinf-2
end if 552 552 end if
553 553
if (izinf==1) then 554 554 if (izinf==1) then
basez=0 555 555 basez=0
else 556 556 else
basez=izinf-2 557 557 basez=izinf-2
end if 558 558 end if
559 559
! We first make 9 one dimensional interpolation on variable x. 560 560 ! We first make 9 one dimensional interpolation on variable x.
! results are stored in ttmp 561 561 ! results are stored in ttmp
do i=1,3 562 562 do i=1,3
do j=1,3 563 563 do j=1,3
ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j)) 564 564 ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
end do 565 565 end do
end do 566 566 end do
567 567
! We then make a 2 dimensionnal interpolation on variables y and z 568 568 ! We then make a 2 dimensionnal interpolation on variables y and z
fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, & 569 569 fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &
3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp) 570 570 3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
end function 571 571 end function
572 572
573 573
574 574
575 575
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 576 576 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 577 577 !
! Akima spline interpolation and spline evaluation 578 578 ! Akima spline interpolation and spline evaluation
! 579 579 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 580 580 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581 581
! Single precision 582 582 ! Single precision
subroutine fvn_s_akima(n,x,y,br,co) 583 583 subroutine fvn_s_akima(n,x,y,br,co)
implicit none 584 584 implicit none
integer, intent(in) :: n 585 585 integer, intent(in) :: n
real, intent(in) :: x(n) 586 586 real, intent(in) :: x(n)
real, intent(in) :: y(n) 587 587 real, intent(in) :: y(n)
real, intent(out) :: br(n) 588 588 real, intent(out) :: br(n)
real, intent(out) :: co(4,n) 589 589 real, intent(out) :: co(4,n)
590 590
real, allocatable :: var(:),z(:) 591 591 real, allocatable :: var(:),z(:)
real :: wi_1,wi 592 592 real :: wi_1,wi
integer :: i 593 593 integer :: i
real :: dx,a,b 594 594 real :: dx,a,b
595 595
! br is just a copy of x 596 596 ! br is just a copy of x
br(:)=x(:) 597 597 br(:)=x(:)
598 598
allocate(var(n+3)) 599 599 allocate(var(n+3))
allocate(z(n)) 600 600 allocate(z(n))
! evaluate the variations 601 601 ! evaluate the variations
do i=1, n-1 602 602 do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) 603 603 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do 604 604 end do
var(n+2)=2.e0*var(n+1)-var(n) 605 605 var(n+2)=2.e0*var(n+1)-var(n)
var(n+3)=2.e0*var(n+2)-var(n+1) 606 606 var(n+3)=2.e0*var(n+2)-var(n+1)
var(2)=2.e0*var(3)-var(4) 607 607 var(2)=2.e0*var(3)-var(4)
var(1)=2.e0*var(2)-var(3) 608 608 var(1)=2.e0*var(2)-var(3)
609 609
do i = 1, n 610 610 do i = 1, n
wi_1=abs(var(i+3)-var(i+2)) 611 611 wi_1=abs(var(i+3)-var(i+2))
wi=abs(var(i+1)-var(i)) 612 612 wi=abs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.e0) then 613 613 if ((wi_1+wi).eq.0.e0) then
z(i)=(var(i+2)+var(i+1))/2.e0 614 614 z(i)=(var(i+2)+var(i+1))/2.e0
else 615 615 else
z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) 616 616 z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
end if 617 617 end if
end do 618 618 end do
619 619
do i=1, n-1 620 620 do i=1, n-1
dx=x(i+1)-x(i) 621 621 dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd 622 622 a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd
b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd 623 623 b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i) 624 624 co(1,i)=y(i)
co(2,i)=z(i) 625 625 co(2,i)=z(i)
!co(3,i)=-(a-3.*b)/dx**2 ! méthode wd 626 626 !co(3,i)=-(a-3.*b)/dx**2 ! méthode wd
!co(4,i)=(a-2.*b)/dx**3 ! méthode wd 627 627 !co(4,i)=(a-2.*b)/dx**3 ! méthode wd
co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx ! méthode JP Moreau 628 628 co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx ! méthode JP Moreau
co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2 ! 629 629 co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2 !
! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6 630 630 ! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6
! etrangement la fonction csval corrige et donne la bonne valeur ... 631 631 ! etrangement la fonction csval corrige et donne la bonne valeur ...
end do 632 632 end do
co(1,n)=y(n) 633 633 co(1,n)=y(n)
co(2,n)=z(n) 634 634 co(2,n)=z(n)
co(3,n)=0.e0 635 635 co(3,n)=0.e0
co(4,n)=0.e0 636 636 co(4,n)=0.e0
637 637
deallocate(z) 638 638 deallocate(z)
deallocate(var) 639 639 deallocate(var)
640 640
end subroutine 641 641 end subroutine
642 642
! Double precision 643 643 ! Double precision
subroutine fvn_d_akima(n,x,y,br,co) 644 644 subroutine fvn_d_akima(n,x,y,br,co)
645 645
implicit none 646 646 implicit none
integer, intent(in) :: n 647 647 integer, intent(in) :: n