fractalsurface.m

来自「分形计算的几个程序原码」· M 代码 · 共 47 行

clc
clear
%x=1:.1:10;
%y=1:.1:10;
t=0:0.1:2*pi;
r=10;
x=r*cos(t);
y=r*sin(t);
if x==0
    x=eps;
end
if y==0
    y=eps;
end
[x,y]=meshgrid(x,y);

D=2.4;
G=1.36;
r=1.5;
M=10;
Ls=1;
L=5;
z=zeros(size(x));
nmax=fix(log10(L/Ls)/log10(r));
z1=L*(G/L).^(D-2).*(log(r)/M)^(1/2);
ff=rand(M,nmax)*(2*pi);
for m=1:M
    for n=1:nmax
        z2=r^((D-3)*n);
        z3=cos(ff(m,n));
        z4=cos((2*pi*r^n*(x.^2+y.^2).^(1/2))/L*cos(abs(atan(abs(y./x))-pi*m/M))+ff(m,n));
        z=z1*z2*(z3-z4);
    end
end
for i=1:length(x)
    for j=1:length(y)
        if (x(i,j).^2+y(i,j).^2)>100
            z(i,j)=0;
        end
        if z(i,j)<0
            z(i,j)=0;
        end
    end
end
view(45,45);
surface(x,y,z);