function griezeta(gamma,col)
%gamma=1 does the gamma function gamma=0 does Riamann zeta function
%col=0 3 colour real imaginary and argument 
%col=1 2 colour modulus and argument
nx = 640;
ny = 640;
ColorMset = zeros(ny,nx,3);
b=0.5+14.1347i;
xmin = -40;
xmax = 40;
ymin = -40;
ymax = 40;
kk=1;
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);
   c=complex(cx,cy);
   
   if gamma
        k=mygamma(c);
   else
        k=myzeta(c);
   end
    if k == 0 
        ColorMset(iy,ix,:) = 0;
    else
        if col
           ColorMset(iy,ix,1) = abs(sin(kk*abs(k)));
           ColorMset(iy,ix,2) = (1+sin(angle(k)))/2;
        else
           ColorMset(iy,ix,1) = (1+cos(kk*real(k)))/2;
           ColorMset(iy,ix,2) = (1+sin(angle(k)))/2;
           ColorMset(iy,ix,3) = (1+cos(kk*imag(k)))/2;
        end
    end
  end
  waitbar(iy/ny,wb)
end
close(wb);
image(ColorMset);
imwrite(ColorMset,'mand.jpg','jpg','Quality',100);

% 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