function myrizLSMand(gamma)
nx = 60;
ny = 160;
ColorMset = zeros(ny,nx,3);
if gamma
    xmin=-10;
    xmax=20;
    ymin=-40;
    ymax=40;
else
    xmin=-25;
    xmax=5;
    ymin=-40;
    ymax=40;
end    
maxiter = 450;
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)
zz = complex(cx,cy);
a=-0.2959050055752;
b=0.5+14.1347i; 
iter = 0;
%z=-2.6;
%z=b;
%z=zz;
%z=999;
%z=-15.34;
if gamma
    z=1/4;
else
    z=999;
end
while (iter < maxiter)&&(abs(real(z))<1000)%&&(~isnan(z))
   if gamma
       z=mygamma(z)+zz;
   else
       z=myzeta(z)+zz;
   end
   iter = iter+1;
end
if iter < maxiter
    if (real(z) >= 100)%||isnan(z)
        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)
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