function mandelset(f,c_1,r,s,color)
%Function to iterate original: f(z,c)=z^2+c
%c0 = complex initial value
%r = radius of window
%s = Steps (should be an even number)
%This function will plot the Mandelbrot set
%For z^2+c or other functions.
%z_0=c_0
%color parameter 0-13; (0=Default, 1=Jet, 2=HSV, 
%3=Hot, 4=Cool, 5=Spring, 6=Summer, 7=Autumn, 8=Winter,
%9=Gray, 10=Bone, 11=Copper, 12=Pink, 13=Lines)

%Make sure that step number is even.
if (s/2) ~= floor(s/2)
    s=s-1;
end
%Make sure the radius is positive and real
r=abs(r);

%calculate step value for s steps based on r
m=abs(r/(s/2));
%calculate position of lower, left point in window
c_min=c_1-(r+r*1i);

%Initialize plot and colormap
       
switch color
    case 0; colormap('default')
    case 1; colormap('Jet')
    case 2; colormap('HSV')
    case 3; colormap('Hot')
    case 4; colormap('Cool')
    case 5; colormap('Spring')
    case 6; colormap('Summer')
    case 7; colormap('Autumn')
    case 8; colormap('Winter')
    case 9; colormap('Gray')
    case 10;colormap('Bone')
    case 11;colormap('Copper')
    case 12;colormap('Pink')
    case 13;colormap('Lines')   
end
newplot
reset(gca);
cm=colormap;
%Determin Maximum length - default colormap=64, others are 256
Max=length(cm);
hold on;
%Initialize N(epsilon) for box-counting dimension
N_eps=0;

%Loop the y values of x+yi
for y=0:s
    z_0=c_min+y*m*1i;

    %Loop the x values of x+yi
    for x=0:s
        k=y*(s+1)+x+1;

        %initialize variables for iteration
        z=z_0+x*m;
        c=z;
        %initialize P to record orbit points
        P(1)=z;
        n=2;
        d=abs(P(1));
        
        %Store initial point in C
        C(k,1)=z; 
        
        %While loop to count the number of iterations to leave the 
        %Max radius(d) of 2, or have a orbit greater than 64 
        %(greatest color value)
        while n<=Max && d<=2
            z=f(z,c);
            P(n)=z;
            d=abs(P(n));
            n=n+1;         
        end
        %uncomment to view orbit of specific points
        %text = ['point: ',num2str(z_1),' length: ',num2str(length(P))];
        %disp (text)
        %store length of orbit
        C(k,2)=length(P);
        clear P;

        %increment N(epsilon) if region attracted to 0.
        if C(k,2)==Max
            N_eps=N_eps+1;
        end
        
        %determine color of patch
        orbs=C(k,2);
        color=[cm(orbs,1),cm(orbs,2),cm(orbs,3)];    
        %Determine location and color of point
        x2=real(C(k,1));
        y2=imag(C(k,1));        
        %Generate a square of apothem (ap) of m/2
        ap=m/2;
        X=[x2-ap;x2-ap;x2+ap;x2+ap];
        Y=[y2-ap;y2+ap;y2+ap;y2-ap]; 
        patch(X,Y,color,'edgeColor',color) 
    end 
end    
hold off;

%Calculate and display approximate box-counting dimension
boxdim=(log(N_eps)/log(1/m));
disp('Approximate Box-Counting Dimension')
boxdim
fname=func2str(f);
%%%[ret,fname]=f(0,0);
%format plot window
axis image;
title (['Function: ',fname, ';  Center= ',num2str(c_0),';   Radius= ',num2str(r),';  Steps= ', num2str(s), ';BoxDim= ',num2str(boxdim)]);

%save tiff
%    name = ['z^2+c_',num2str(c_0),'_',num2str(r),'_',num2str(s),'.tiff'];
%   print ('-dtiff', name)
%save jpeg
%    name = ['z^2+c_',num2str(c_0),'_',num2str(r),'_',num2str(s),'.jpg'];
%    print ('-djpeg', name)
end

% % function [ret,fname]=f(z,c)
% % ret=z^2+c;
% % fname='z^2+c';
% % end
