clear
clf

% RESOLUTION/INITIAL FACTORS
x= 0:.02:2*pi;   % smaller increment = greater resolution
y= 0:.02:2*pi;   % .1 for decent resolution, .02 for high resolution
[Y,X] = meshgrid(y,x);  % high resolution = LONG processing time
Z= zeros(length(x),length(y));
r=1;    % can set a radius to planet
d=1;    % can set an elevation increase/decrease factor
f=0;    % counts frames for gif image.

for t=1:10000           % 13+ HOURS TO PROCESS @ 10,000 ITERATIONS!!!
    
    % GENERATE RANDOME SINE CURVE
    sign= randomsign(); % randomsign generates a number -1 or 1
    a= rand*sign*pi;    % 'a' is the amplitude from -pi to pi
    sign= randomsign(); 
    b= rand*sign*2*pi;      % 'b' is the phase from 0 to 2*pi
    line= a*sin(x-b)+pi;    % generates a sine curve with random amplitude and phase
    sign= randomsign();
    elev1= sign*d*rand();   % generates elevation for values above sine curve
    sign= randomsign();
    elev2= sign*d*rand();   % generates elevation for values below sine curve
    
    % DETERMINE ELEVATION CHANGES
    for i=1:length(x)       % 'for' loop goes over all values of Z
        for j=1:length(y)   % then uses the distance formula with 'i,j' and 'i,line'
            if sqrt((2*pi*i/length(x))^2+(2*pi*j/length(y))^2) > sqrt((2*pi*i/length(x))^2+line(i)^2)
                Z(i,j)=Z(i,j)+elev1;    % these values are above the sine curve
            else
                Z(i,j)=Z(i,j)+elev2;    % these values are below the sine curve
            end
        end
    end

    if mod(t,5)==0              % tells me how far along in the code is
        fprintf('t= %d\n',t);   % by printing 't' every 5th iteration
    end
    
    Z2= Z/t;    % 'Z2' is now the normalized version of Z.  Could use 'Z', but didn't
    for i=1:length(x)   % goes over every element of Z2, and multiplies by 'r'
        for j=1:length(y)
            Z2(i,j)= (Z2(i,j))*r;   % this way you can have a planet with a determined radius
        end
    end
    
    % MAKE MOVIE
    if t<5 || (mod(t+9,10)==0 && t<100) || (mod(t+25,25)==0 && t<500) || ... % when to draw/save an image
            (mod(t+50,50)==0 && t<1000) || (mod(t+120,120)==0 && t<10000) || t==10000
        f=f+1;
        subplot(1,2,1);     % 1 by 2 subplot for topological view and spherical view
        surf(X,Y,Z2, 'EdgeColor', 'none')
        axis tight
        view(2)     % topological view goes on first plot with a top-down view
        subplot(1,2,2);     % second subplot shows the spherical view
        [xx,yy,z] = sphere(100);
        cla reset
        axis square off
        props.AmbientStrength = 0.1;       % this is all graphical modifications
        props.DiffuseStrength = 1;
        props.SpecularColorReflectance = .5;
        props.SpecularExponent = 20;
        props.SpecularStrength = 1;
        props.FaceColor= 'texture';
        props.EdgeColor = 'none';
        props.FaceLighting = 'phong';
        props.Cdata = Z2';
        surface(xx,yy,z,props);
        %light('position',[-1 0 1]);   % adds lighting for a more 3D effect
        %light('position',[-1.5 0.5 -0.5]);%, 'color', [.6 .2 .2]);
        view(-37.5+4*360/t, 30)     % view will rotate 4 times over course of animation
        set(figure(1), 'position', [10, 540, 1200, 400], 'PaperPositionMode', 'auto')   % sets the size and location of figure
        demcmap(Z2)
        drawnow
        saveas(gcf, ['fractal_planet',sprintf('%04d',f),'.png'])
    end
end

save('Zvalue2.mat', 'Z2')
 