二维三次卷积插值算法及Fortran代码

cyhwjie

==============================J*e*r*k*w*i*n*@*g*m*a*i*l*.*c*o*m=============================
  1    Program CubConv
  2        integer, parameter:: NgrdX=5, NgrdY=5
  3        integer i, j, Nx, Ny
  4        real*8  Xmin, Ymin, Xmax, Ymax, Xgrd, Ygrd, x, y, dx, dy, Rtmp, &
  5        &       Zxy(0:NgrdX+1, 0:NgrdY+1),dZx(NgrdX, NgrdY), dZy(NgrdX, NgrdY), &
  6        &       F, CubConvInterp
  7   
  8        F(x, y) = 0.5D0*(x*x+y*y/20)
  9   
10        Xgrd = 1.D0
11        Ygrd = 1.D0
12        Xmin = -2.D0
13        Ymin = -2.D0
14        Xmax = Xmin+(NgrdX-1)*Xgrd
15        Ymax = Ymin+(NgrdY-1)*Ygrd
16   
17        do i=1, NgrdX
18            x = Xmin+dble(i-1)*Xgrd
19            do j=1, NgrdY
20                y = Ymin+dble(j-1)*Ygrd
21                Zxy(i, j) = F(x, y)
22            end do
23        end do
24   
25        Zxy(0, :) = 2.D0*Zxy(1, :)-Zxy(2, :)
26        Zxy(NgrdX+1, :) = 2.D0*Zxy(NgrdX, :)-Zxy(NgrdX-1, :)
27   
28        Zxy(:, 0) = 2.D0*Zxy(:, 1)-Zxy(:, 2)
29        Zxy(:, NgrdY+1) = 2.D0*Zxy(:, NgrdY)-Zxy(:, NgrdY-1)
30   
31    ! 节点导数
32        ! do i=1, NgrdX
33            ! do j=1, NgrdY
34                ! dZx(i, j) = 0.5D0*( Zxy(i+1, j)-Zxy(i-1, j) )/Xgrd
35                ! dZy(i, j) = 0.5D0*( Zxy(i, j+1)-Zxy(i, j-1) )/Ygrd
36                ! write(*, '(2I4, 2F8.3)') i, j, dZx(i, j), dZy(i,j)
37            ! end do
38        ! end do
39   
40    ! 改变网格，测试插值
41        Nx=50
42        Ny=50
43        dx = (Xmax-Xmin)/dble(Nx-1)
44        dy = (Ymax-Ymin)/dble(Ny-1)
45        do i=1, Nx
46            x = Xmin+dble(i-1)*dx
47            do j=1, Ny
48                y = Ymin+dble(j-1)*dy
49                Rtmp = CubConvInterp(x, y, Xmin, Ymin, Xgrd, Ygrd, NgrdX, NgrdY, Zxy)
50                write(1, '(4F10.3)') x, y, Rtmp, F(x, y)
51            end do
52        end do
53   
54    End Program CubConv
55   
56    Real*8  Function CubConvInterp(Xpos, Ypos, Xmin, Ymin, Xgrd, Ygrd, Nx, Ny, Zxy)
57        integer i, j, Nx, Ny
58        real*8  Xpos, Ypos, Xmin, Ymin, Xgrd, Ygrd, u, v, Zxy(0:Nx+1, 0:Ny+1), &
59        &       X(4), Y(4), C(4, 4), D(4, 4), DX(4), DY(4), CDX(4)
60   
61        data D(1, 1:4) / -1,  2, -1, 0 /
62        data D(2, 1:4) /  3, -5,  0, 2 /
63        data D(3, 1:4) / -3,  4,  1, 0 /
64        data D(4, 1:4) /  1, -1,  0, 0 /
65   
66        i = int((Xpos-Xmin)/Xgrd)+1
67        j = int((Ypos-Ymin)/Ygrd)+1
68        u = Xpos-(Xmin+dble(i-1)*Xgrd)
69        v = Ypos-(Ymin+dble(j-1)*Ygrd)
70        ! print*, Xpos, Ypos, i, j, u, v
71        ! pause
72   
73        C(1, 1:4) = [ Zxy(i-1, j-1), Zxy(i, j-1), Zxy(i+1, j-1), Zxy(i+2, j-1) ]
74        C(2, 1:4) = [ Zxy(i-1, j  ), Zxy(i, j  ), Zxy(i+1, j  ), Zxy(i+2, j  ) ]
75        C(3, 1:4) = [ Zxy(i-1, j+1), Zxy(i, j+1), Zxy(i+1, j+1), Zxy(i+2, j+1) ]
76        C(4, 1:4) = [ Zxy(i-1, j+2), Zxy(i, j+2), Zxy(i+1, j+2), Zxy(i+2, j+2) ]
77   
78        X(1:4) = [u*u*u, u*u, u, 1.D0]
79        Y(1:4) = [v*v*v, v*v, v, 1.D0]
80   
81        DX = matmul(D, X)
82        DY = matmul(D, Y)
83        CDX = matmul(C, DX)
84        CubConvInterp = 0.25D0*( DY(1)*CDX(1)+DY(2)*CDX(2)+DY(3)*CDX(3)+DY(4)*CDX(4) )
85    End Function CubConvInterp
============================================================================================

本文引用地址：http://blog.sciencenet.cn/blog-548663-595297.html

topmad

没看懂那个  1.DO 是啥....
插值的算法好多  一个也没弄明白

郭小小

{:5_235:}{:5_235:}{:5_235:}{:5_235:}{:5_235:}{:5_235:}

最爱小小罗

顶个~~~~~~~

最爱小小罗

顶个~~~~~~~

最爱小小罗

顶个~~~~~~~

Anonymous888

谢谢分享！！！