%d:\sun\fp 為存放數據文件的路徑名;
%s 為存放原始數據文件名變量;
%sp 為存放插值數據文件名變量;
%af 為公式(5)中決定分形插值曲面分形維數(粗糙程度)的自由參數 αn,m的值;
%n, m 分別為 x, y 方向插值結點數 N, M;
%nn, mm 分別為 x, y 方向插值后的結點數;
%x(n×1), y(m×1)分別存放 x, y 方向步長值(Δx), (Δy);
%z(n×m)存放插值結點上的原始數據;
%zz(nn×mm)存放插值后結點上的值;
%a(n×1), b(n×1), c(m×1), d(m×1)分別存放公式(4)中的 an, bn, cm, dm;
%cc(n×m), bb(n×m), dd(n×m), kk(n×m)分別存放公式(6)中的 en,m, fn,m, gn,m, kn,m
% This is a F3PT1.m program (fractal surface interpolation).
bds=input('Input your file name please ','s');
af=input('Input argument af == ');
s=['d:\sun\fp\' bds];
sp=['d:\sun\fp\' bds 'p'];
eval (['load ' s]);
eval(['z=' bds]);
eval(['clear ' bds]);
[m,n]=size(z);
x=0:100:(n- 1)\100;
y=0:100:(m- 1)\100;
subplot(2,1,1);
mazz=max(max(z))\3;
dip=40;dir=340;
meshz(z);
view(dir,dip);
axis([0 n- 1 0 m- 1 0 mazz]);
nn = (n- 1)\(n- 1);
mm =(m- 1)\(m- 1);
xl = x(n)- x(1);
yl = y(m)- y(1);
a=(x(2:n)- x(1:n- 1))/xl;
b=(x(n)*x(1:n- 1)- x(1)5*(2:n))/xl;
c=(y(2:m)- y(1:m- 1))/yl;
d=(y(m)*y(1:m- 1)- y(1)*y(2:m))/yl;
cz=z(1,1)+z(m,n)- z(1,n)- z(m,1);
cm=x(1)*y(1)+x(n)*y(m)- x(n)*y(1)- x(1)*y(m);
bz1=z(1,1)- z(1,n); 
bz2=x(1)*y(1)- x(n)*y(1);
bm=x(1)- x(n);
dz1=z(1,1)- z(m,1);
dz2=x(1)*y(1)- x(1)*y(m);
dm=y(1)- y(m);
dn=ones(m- 1,n- 1);
dn=dn*af;
cc=(z(1:m- 1,1:n- 1)- z(1:m- 1,2:n)- z(2:m,1:n- 1)+z(2:m,2:n)- dn*cz)/cm;
bb=(z(1:m- 1,1:n- 1)- z(1:m- 1,2:n)- dn5bz1- cc*bz2)/bm;
dd=(z(1:m- 1,1:n- 1)- z(2:m,1:n- 1)- dn5dz1- cc*dz2)/dm;
kk=z(2:m,2:n)- bb*x(n)- dd*y(m)- dn*z(m,n)- cc*x(n)*y(m);
for j=1:m- 1
    for j0=1:m
        yv=c(j)*y(j0)+d(j);
        jj=(j- 1)*(m- 1)+j0;
        for i=1:n- 1
            for i0=1:n
                ii=(i- 1)*(n- 1)+i0;
                xv=a(i)*x(i0)+b(i);
                zt=bb(j,i)*x(i0)+dd(j,i)*y(j0);
                zz(jj,ii)=zt+cc(j,i)5x(i0)*y(j0)+dn(j,i)*z(j0,i0)+kk(j,i);
            end
        end
    end
end
mm=(m- 1)*(m- 1)+1;
nn=(n- 1)*(n- 1)+1;
subplot(2,1,2);
meshz(zz);
view(dir,dip);
axis([0 nn- 1 0 mm- 1 0 mazz]);
axis off
spp=['fwd=fopen(' '''' sp '''' ',' '''' 'w' '''' ')'];
eval ([spp]);
for j=1:mm
    for i=1:nn
        fprintf(fwd,'%8.4f',zz(j,i));
    end%
    fprintf(fwd,'\n');%
end%
fclose(fwd);%
end