function myrizLSM2c(maxiter,huge,v,ww,tiny)
%example myrizLSM2c(400,100,0,1,0.001)
%maxiter=max iterations in loop, huge=size of large target for divergence
%v=parameter value (i.e. c), ww=flag to write wait bar or not
%tiny=size of small target distance for checking periodicity
%maxiter=400; %huge=100; %tiny=.001; %suitable values 
%Adjust for speed or accuracy
nx = 320;
ny = 320;
fprintf('%4.5f %4.5f\n',real(v),imag(v));
ColorMset = zeros(nx,ny,3);
b=0.5+14.1347i; 
% xmin = -1;
% xmax = 1;
% ymin = 15.1347;
% ymax = 13.1347;
xmin=-40;
xmax=40;
ymin=-40;
ymax=40;
%maxiter=400;
kk=2/10;
orbit=zeros(2,maxiter);
if ww
    wb = waitbar(0,'Please wait...');
end
for iy = 1:ny
  cy = ymin + iy*(ymax - ymin)/(ny - 1);
  for ix= 1:nx
    cx = xmin + ix*(xmax - xmin)/(nx - 1);
   [k l] = Mlevel(cx,cy,maxiter,v,huge,tiny);
    if k == 0 
        ColorMset(iy,ix,:) = 0;
    else
        if l==1
            ColorMset(iy,ix,1) = abs(sin(kk*k));
        end
        if l==2
            ColorMset(iy,ix,2) = abs(sin(kk*k));
        end
        if l==3
            ColorMset(iy,ix,3) = abs(sin(kk*k));
        end
        if l==0
            ColorMset(iy,ix,:) = abs(sin(kk*k))/3;
        end
     end
  end
  if ww
    waitbar(iy/ny,wb)
  end
end
if ww
    close(wb);
end
image(ColorMset);
text(20,460,num2str(v),'color','white');
imwrite(ColorMset,strcat(num2str(v),'.jpg'),'jpg','Quality',100);

function [potential color] = Mlevel(cx,cy,maxiter,v,huge,tiny)
z = complex(cx,cy);
a=-0.2959050055752;
b=0.5+14.1347i; 
%huge=100;
%tiny=.001;
color=1;
myflag=0;
for iter=1:maxiter
   z=myzeta(z)+v;
   orbit(1,iter)=real(z);
   orbit(2,iter)=imag(z);
   if abs(real(z))>huge||isnan(z)
       break
   end
   myloop=min(10,iter-1);
   re=real(z);
   im=imag(z);
   if real(z)>10
       color=2;
   end
   pflag=1;
   for i=iter-1:-1:iter-myloop
       re2=re-orbit(1,i);
       if orbit(1,i)<0
           pflag=0;
       end
       im2=im-orbit(2,i);
       if re2*re2+im2*im2<tiny
           if pflag==1
               color=2;
           end
           myflag=1;
           break
       end
   end
   if myflag==1
       break
   end
end
potential = iter;
if myflag==0
    color=0;
end
if abs(real(z))>huge||isnan(z)
    color=3;
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