%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);
 x1=x(n)-x(1);
 y1=y(m)-y(1);
 a=(x(2:n)-x(1:n-1))/x1;
 b=(x(n)*x(1:n-1)-x(1)*x(2:n))/x1;
 c=(y(2:m)-y(1:m-1))/y1;
 d=(y(m)*y(1:m-1)-y(1)*y(2:m))/y1;
 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)-dn*bz1-cc*bz2)/bm;
 dd=(z(1:m-1,1:n-1)-z(2:m,1:n-1)-dn*dz1-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)*x(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