function mygrizLSM(gamma)
%gamma=1 draws the julia set of the gamma function 
%gamma=0 draws the julia set of the Riemann zeta function
%level set method for the Julia set of zeta separately plotting basin of
%attraction of infinity and the fixed attractor a
nx = 640;
ny = 640;
ColorMset = zeros(nx,ny,3);
if gamma
    xmin=-5;
    xmax=15;
    ymin=-10;
    ymax=10;
else
    xmin = -40;
    xmax = 40;
    ymin = -40;
    ymax = 40;
end
maxiter = 100;
wb = waitbar(0,'Please wait...');
for iy = 1:ny
  cy = ymin + iy*(ymax - ymin)/(ny - 1);
  for ix= 1:nx
    cx = xmin + ix*(xmax - xmin)/(nx - 1);
    k = Mlevel(cx,cy,maxiter,gamma);
    if k == 0 
        ColorMset(iy,ix,:) = 0;
    else    
        ColorMset(iy,ix,1) = abs(sin(2*k/10));
        ColorMset(iy,ix,2) = abs(sin(2*k/10+pi/4));
        ColorMset(iy,ix,3) = abs(cos(2*k/10));
    end
  end
  waitbar(iy/ny,wb)
end
close(wb);
image(ColorMset);
imwrite(ColorMset,'mand.jpg','jpg','Quality',100);

function [potential] = Mlevel(cx,cy,maxiter,gamma)
z = complex(cx,cy);
if gamma
    a=1;
else
    a=-0.2959050055752;
end
iter = 0;
huge=1000;
%while (iter < maxiter)&&(real(z)<huge)&&(abs(z-0.4971)>.01)
while (iter < maxiter)&&(real(z)<huge)&&(abs(z-a)>.01)
%while (iter < maxiter)&&(abs(z)<huge)&&(abs(z)>.01)
   %z=myzetas(z,1002)*myzetas(z-1,1002);
   if gamma
       z=mygamma(z);
   else
       z=myzeta(z,2002);
   end
   %z=mygamma(z/2+1)*(z-1)*pi^(-z/2)*myzeta(z);
   iter = iter+1;
end
if iter < maxiter
    if real(z) >= huge
        potential = iter;
    else
        potential = maxiter-iter;
    end
else
    potential = 0;
end

function f=mygamma2(z,its)
%alternative gamma function
f=1/z;
for n=1:its
    f=f*((1+1/n)^z)/(1+z/n);
end

function [f]=mygamma(z)
% complex gamma function using Lanczos approximation to extend to the
% entire plane. 

pi=3.14159;
twopi=pi+pi;
c = [   1.000000000000000174663;
     5716.400188274341379136;
   -14815.30426768413909044;
    14291.49277657478554025;
    -6348.160217641458813289;
     1301.608286058321874105;
     -108.1767053514369634679;
        2.605696505611755827729;
       -0.7423452510201416151527e-2;
        0.5384136432509564062961e-7;
       -0.4023533141268236372067e-8];
g=9;
s=0;
zz=z;
%beep
if ((round(zz)==zz)&&(imag(zz)==0)&&(real(zz)<=0))
	%f=10^308;
	f=Inf;
%end
else
% if ((round(zz)==zz)&&(imag(zz)==0)&&(real(zz)>=0))
% 	fac=1;
% 	for k=1:zz-1
% 		fac=fac*k;
%     end
% 	f=fac;
% end
%if ((round(zz)~=zz)||(imag(zz)~=0))
    if real(z)<0
        z=-z;
    end
    t=z+g;
    for k=g+2:-1:2
        s=s+c(k)/t;
        t=t-1;
    end
    s=s+c(1);
    ss=(z+g-0.5);
    s=log(s*sqrt(twopi))+(z-0.5)*log(ss)-ss;
    LogofGamma = s;
    f = exp(LogofGamma);
    if real(zz)<0
        f=-pi/(z*f*sin(pi*z));
    end
end

function f=myzeta(z,steps)
zz=z;
if real(z)<0
    z=1-z;
end
%steps=2002;
f=0;
sw=-1;
for i=1:steps
    sw=-sw;
    hold=sw/i^z;
    f=f+hold;
    if abs(hold)<.001
        break
    end
end
sw=-sw;
f=f+sw/(i+1)^z/2;
f=f/(1-2^(1-z));
if real(zz)<0
    f=2*(2*pi)^(-z)*cos(z*pi/2)*mygamma(z)*f;
end