clf;
clear;
%f=0.000001 N=100 A=1/N^2 tmax=5000 with tall vertical axis.
%F=0.00001 N=50 A=0.33/N^2 tmax=5000 with short axis.
N=100; %100 is good, but slow.
k=0.005; %scaling factor for random position changes. .002 is ok.
A=0.33/N^2;%scaling factor for changes in the Ls.
M=zeros(3,N,N); % creates an NxN array of 3D vectors representing grid points.
f=0.00001; %surface bending energy parameter. 
u=2.5; %power on surface bending term
nt=0;
for i=1:N
    for j=1:N
        M(1,i,j)=i; %sets the initial position vectors to fall on a grid.
        M(2,i,j)=j;
    end
end
LH=zeros(N,N-1);
LV=zeros(N-1,N);
for i=1:N
    for j=1:N-1
        LH(i,j)=norm(M(:,i,j)-M(:,i,j+1))+A*i^2; %sets the initial distances between horizontally adjacent points.
    end
end
for i=1:N-1
    for j=1:N
        LV(i,j)=norm(M(:,i,j)-M(:,i+1,j));%+A*i^2; %sets the initial distances between vertically adjacent points.
    end
end
t=0;
while t<50000
    t=t+1;
    for i=1:N
        for j=1:N
            %computes the initial local energy (the energy for each adjacency for a given point).
            if i==1
                if j==1
                    E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2;
                elseif j==N
                    E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2;
                else
                    E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i,j-1)),(M(:,i,j+1)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i,j-1))*norm(M(:,i,j+1)-M(:,i,j)))))^u; 
                end
            elseif i==N
                if j==1
                    E=(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2;
                elseif j==N
                    E=(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2;
                else
                    E=(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i,j-1)),(M(:,i,j+1)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i,j-1))*norm(M(:,i,j+1)-M(:,i,j)))))^u; 
                end
            elseif j==1
                E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i-1,j)),(M(:,i+1,j)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i-1,j))*norm(M(:,i+1,j)-M(:,i,j)))))^u; 
            elseif j==N
                E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i-1,j)),(M(:,i+1,j)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i-1,j))*norm(M(:,i+1,j)-M(:,i,j)))))^u; 
            else
                E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i-1,j)),(M(:,i+1,j)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i-1,j))*norm(M(:,i+1,j)-M(:,i,j)))))^u+...
                        f*(acosd( dot((M(:,i,j)-M(:,i,j-1)),(M(:,i,j+1)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i,j-1))*norm(M(:,i,j+1)-M(:,i,j)))))^u;
            end
            Eold=E;
            Xtemp=M(1,i,j);
            Ytemp=M(2,i,j);
            Ztemp=M(3,i,j);
            %Modifies all the coordinates of each point by a small random value between k and negative k.
            M(1,i,j)=M(1,i,j)+(2*rand()-1)*k;
            M(2,i,j)=M(2,i,j)+(2*rand()-1)*k;
            M(3,i,j)=M(3,i,j)+(2*rand()-1)*k;
            %Computes the new energy.
           if i==1
                if j==1
                    E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2;
                elseif j==N
                    E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2;
                else
                    E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i,j-1)),(M(:,i,j+1)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i,j-1))*norm(M(:,i,j+1)-M(:,i,j)))))^u; 
                end
            elseif i==N
                if j==1
                    E=(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2;
                elseif j==N
                    E=(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2;
                else
                    E=(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i,j-1)),(M(:,i,j+1)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i,j-1))*norm(M(:,i,j+1)-M(:,i,j)))))^u; 
                end
            elseif j==1
                E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i-1,j)),(M(:,i+1,j)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i-1,j))*norm(M(:,i+1,j)-M(:,i,j)))))^u; 
            elseif j==N
                E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i-1,j)),(M(:,i+1,j)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i-1,j))*norm(M(:,i+1,j)-M(:,i,j)))))^u; 
            else
                E=(norm(M(:,i,j)-M(:,i+1,j))-LV(i,j))^2+(norm(M(:,i,j)-M(:,i-1,j))-LV(i-1,j))^2+(norm(M(:,i,j)-M(:,i,j-1))-LH(i,j-1))^2+(norm(M(:,i,j)-M(:,i,j+1))-LH(i,j))^2+...
                        f*(acosd( dot((M(:,i,j)-M(:,i-1,j)),(M(:,i+1,j)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i-1,j))*norm(M(:,i+1,j)-M(:,i,j)))))^u+...
                        f*(acosd( dot((M(:,i,j)-M(:,i,j-1)),(M(:,i,j+1)-M(:,i,j)))/(norm(M(:,i,j)-M(:,i,j-1))*norm(M(:,i,j+1)-M(:,i,j)))))^u;
            end
            %Comapres the new energy with the inital energy.
            %If the new energy is great it undoes the move.
            if E>Eold
                M(1,i,j)=Xtemp;
                M(2,i,j)=Ytemp;
                M(3,i,j)=Ztemp;
            end
        end
    end
    if mod(t+499,500)==0
        nt=nt+1;
        X=squeeze(M(1,:,:));  %not sure if squeeze is the right thing.
        Y=squeeze(M(2,:,:));
        Z=squeeze(M(3,:,:));
        colormap cool
        surf(X,Y,Z)
        axis([-N/2 N+N/2 -N/2 N+N/2 -2 2])
        view(40, 30);
        saveas(gcf, ['Diff_growth_',sprintf('%03d',nt),'.png'])
        drawnow
    end
    sprintf('%d',t)
end






