%       program hd2di
%       Current version: August 17, 2017
%       Previous versions: November 30, 2015, May 18, 2017
%       Computes Hausdorff dimension of one of several invariant sets
%       described below:
%         Let z = x + iy
%         Let H:= {(x,y) in R^2: (x-1/2)^2 + y^2 <= 1/4, y >=0},
%            which we consider a subset of the complex plane.
%        Let I_1 = {b=m+ni: m is a positive integer and n is an integer.
%        Let I_2 = {b=m-ni: m and n are positive integers.
%        Let I_3 = {b=m+ni: m =1 or 2 and n = 0,1,2,-1,-2}.
%         Let theta_b(z) = 1/(z+b).
%       For each j=1,2,3, there exists a unique invariant set
%       K_j contained in D such that
%       K_j = closure [Union_{b in I_j} \theta_b(K)].
%       To compute the Hausdorff dimension of K_j, we use the nontrivial 
%       result that it is the unique value of s > 1 for which mu_s =1, where
%       mu_s is the spectral radius of the operator L_s given by
%       L_s v_s = \sum_{b in I_j, Im(b) >=0}
%		   (1/|z+b|^{2s} v_s(1/(\bar z + \bar b))
%                + \sum_{b in I_j, Im(b) < 0}
%		   (1/|z+b|^{2s} v_s(1/(z + b))
%       where \bar z denotes the complex conjugate of z.
%       We reduce this to a matrix eigenvalue problem by approximating the
%       function v_s by a piecewise bilinear function and
%       restricting the set I_1 to values of b for which
%       |b| <= R, for a sufficiently large value of r. 
%   
%      Uses the routines matas.m and thetabz.m
%
%      Change from hd2d.m is to compute bilinear functions only once
%
%    Key variables are described below.
%    Iset = integer (1,2,or 3) corresponding to I_1, I_2, or I_3     
%    zlen = length of eigenvector
%    N = number of subintervals along line y=0;  Choose N to be divisable by 2
%    So N+1 = number of mesh points along the line y=0.
%    h = 1/N
%    R = maximum modulus of b used in the computation
%    Many variables below have the last letter p or m.
%    p is used in the case when Im(b) >=0 and m is used when Im(b) <0.
%    If au is the upper matrix, al is the lower matrix, aa is the approx matrix
%    muau is the largest eigenvalue of au, mual is the largest eigenvalue of al,
%    muaa is the largest eigenvalue of aa
%    If we find values s_1 for which mual >= 1 and s_2 for which muau <= 1,
%    then the true largest eigenvalue (giving the Hausdorff dimension of K)
%    satisfies  s_1 <= s <= s_2.
%  
%      integer scalar variables
%      integer N, N2, N21, R, P, zlen, blen, i, jl, k, l, m
%      integer nn, mm, ll, jlJ, mtype, Iset
%
%      integer vector variables
%      integer J, nrange, nmax, nmin, bm, bn
%
%      integer matrix variables
%      integer ixp,iyp,ixm,iym
%
%      real scalar variables
%      real(dp) h, thetabxpkl, thetabypkl, thetabxmkl, thetabymkl
%      real(dp) aphibpskl, aphibmskl, ir, jlr
%
%      real vector variables
%      real(dp) x,y, Jr, bmr, bnr, tf
%      real(dp) muaa, mual, muau
%      real(dp) muaam1, mualm1, muaum1
%
%      real matrix variables
%      real(dp) aphibp, aphibm
%      real(dp) aphibps, aphibms
%
format long
prompt = 'Input Iset '
Iset = input(prompt)
prompt = 'Input N '
N = input(prompt)
%
%      N = 200
%      N = 100
%      N = 50
%
h = 1.0/double(N)
N2 = N/2;
N21 = N2-1;
if (Iset == 1 | Iset == 2)
     prompt = 'Input R '
     R = input(prompt)
elseif (Iset == 3)
     R=0
end
%
%      R = 100
%      R = 200
%      R = 300
%
%
%    Construct a domain H_h made up of a union of squares of side h
%    such that H_h contains D.  H_h should have the property that
%    if (x,y) is in D, then the four corners of the mesh square of side h
%    containing (x,y) are mesh points in H_h.  Although we could simply
%    choose H_h to be the rectangle [0,1] x [0,.5], that choice would add 
%    unknowns we do not use.
%
%    Compute the number of mesh points inside the closure of H_h
%    for each y = jh, in the form 2*J(j) + 1. So J(j) is the number of 
%    subintervals starting from the left and ending at x=1/2, where j runs from 
%    1, ... N/2 + 1.
%    These mesh points will be of the form ih. 
%    Also compute the total number of mesh points, storing the
%    result in zlen.
%    The method first computes an initial approximation to J(j) based on 
%    finding the closest mesh point outside or on the circle.  This value is
%    then adjusted to give a mesh domain H_h that contains H.
%
J = zeros(N2+1,1);
tf = zeros(N2+1,1);
J(1) = N2;
J(N2+1) = 1;
for jl = 1:N21
  jlr = double(jl);
  tf(jl+1) = sqrt(.25-(jlr*h)^2)/h;
  J(jl+1) = ceil(tf(jl+1));	   
end
%    display(J)
%    This next loop adjusts to give a mesh domain H_h with the property that
%    if (x,y) is in D, then the four corners of the mesh square of side h
%    containing (x,y) are mesh points in H_h.
%    
for jjl = 2:N2+1
  jl = N2+3-jjl;
  J(jl) = J(jl-1);
end
% display(J);
%
%  Recall zlen = sum_j (2*J(j) +1)
%
zlen = 2*sum(J) + length(J);
% display(tf)
% display(J)
display(zlen)
Jr = double(J);
%
%    Next, compute the mesh points in H_h, indexing them
%    by a single index k of length zlen starting at 1.
%    The index runs over the mesh points in each row from left to right
%    starting with the mesh points at level y=0.
%    Compute the mesh points so they are symmetric around x = 1/2.
%
x=zeros(zlen,1);     y=zeros(zlen,1);
ll=1;
for jl = 0:N2
  jlr = double(jl);
  jlJ = 2*J(jl+1);
  for i = 0:J(jl+1)
    ir = double(i);
    x(ll + i) = .5 - J(jl+1)*h + ir*h;
    x(ll + 2*J(jl+1) -i) = .5 + J(jl+1)*h - ir*h;
    y(ll + i) = jlr*h;
    y(ll + 2*J(jl+1) -i) = jlr*h;
  end
%         display([jl, x(ll), y(ll), x(ll+jlJ), y(ll+jlJ)])
  ll = ll + jlJ + 1;
end 
% The following is used to plot the mesh domain
 circx = linspace(0.0, 1.0, 200);
 circy = sqrt(.25 - (circx-.5).^2);
% figure % new figure window
 plot(circx,circy);
 title('Domain H and mesh domain H_h')
 axis('equal')
 axis([-.1, 1.1, -.1, .6])
 hold on
 scatter(x,y,5,'s','fill')
 hold off % reset hold state
%    If Iset = 1 or 2:
%    Compute the values of b = m + in for which m= 1,...,R and m^2 + n^2 <=R^2.
%    Index these points, which we write as b = bm(l) + i bn(l), by a 
%    single index l starting at 1.
%    The index uses an ordering in which values of (m,n) satisfy
%    (P-1)^2 < m^2 + n^2 <= P^2, where P = R, R-1, ... 1.
%    So for $P=R$, we have (m,n) = (R,0), (R-1, nrange(R-1)), (R-1, 1),
%    (R-2, ...)
%    No attempt is made to order the points within the various ranges.
%
%    For each m = 1, ... R, we first compute nrange(m), 
%    the largest value of n between 0 and R
%    for which b= m + ni satisfies  |b| <= R. Also compute the total number of 
%    such values, storing the result in blen.
if Iset == 3
   blen = 6;
   bm = zeros(blen,1);  bn= zeros(blen,1);
   %    These are
   bm(1) = 2;    bn(1) = 2;
   bm(2) = 2;    bn(2) = 1;
   bm(3) = 1;    bn(3) = 2;
   bm(4) = 2;    bn(4) = 0;
   bm(5) = 1;    bn(5) = 1;
   bm(6) = 1;    bn(6) = 0;
   bmr = zeros(blen,1);
   bnr = zeros(blen,1);
   for k  = 1:blen
      bmr(k) = double(bm(k));
      bnr(k) = double(bn(k));
   end
end
if Iset == 1 | Iset == 2
%%%   nrange = zeros(R,1);
   blen = 0;
   for m = 1:R-1
      nrange(m) = floor(sqrt(double(R^2 - m^2)));
      blen = blen + nrange(m) + 1;
   end
   if Iset == 1 
     nrange(R) = 0;
     blen = blen +1;
   end
   if Iset ==2
    blen = blen - R + 1;  % we do not allow the value n=0
   end
display(blen)
%
   bm = zeros(blen,1);  bn= zeros(blen,1);
   nmin=zeros(R,1); nmax=zeros(R,1);
   l=1;
   for PP = 2:R
      P = R - PP + 2;
       if Iset == 1
         bm(l) = P;
         bn(l) = 0;
         l = l+1;
       end
      for m = 1:P-1
         nmin(m) = ceil(sqrt(double((P-1)^2 - m^2)));
         nmax(m) = floor(sqrt(double(P^2-m^2)));
         for n = nmin(m):nmax(m)
            if (m^2 + n^2 > (P-1)^2)
	       bm(l) = m;   bn(l) = n;
%              display([l,bm(l),bn(l)])
	       l = l+1;
             end
         end
       end
   end  
    if Iset == 1
      bm(l) = 1;
      bn(l) = 0;
    end
%    display([l,bm(l),bn(l)])
   bmr = zeros(blen,1);
   bnr = zeros(blen,1);
   for k  = 1:blen
      bmr(k) = double(bm(k));
      bnr(k) = double(bn(k));
   end
end
%
%    For each of the mesh points, compute the functions
%    aphibp = |1/(z + b)|^2 = 1/((x + m)^2 + (y +n)^2)
%    aphibm = |1/(z + bar b)|^2 = 1/((x + m)^2 + (y - n)^2)
%
aphibm= zeros(zlen,blen);
if (Iset == 1 | Iset == 3)
   aphibp= zeros(zlen,blen);
elseif Iset == 2 
   aphibp = zeros(1,1);
end
for k = 1:zlen
  for l = 1: blen
       aphibm(k,l) = 1/((x(k) + bmr(l))^2 + (y(k) - bnr(l))^2);
   if (Iset == 1 | Iset == 3)
    aphibp(k,l) = 1/((x(k) + bmr(l))^2 + (y(k) + bnr(l))^2);
    end
  end
end
%
%    display(aphibp)
%    display(aphibm)
%
%    For each of the points thetab(z) = thetabx + i thetaby, 
%    find the vertices (ix*h,iy*h) of the lower left corner of the 
%    square to which it belongs.
%    This is done for the two choices (thetabx,thetaby) = (thetabxp,thetabyp)
%    and (thetabx,thetaby) = (thetabxm,thetabym) with the results called
%    (ixp*h, iyp*h) and (ixm*h, iym*h), respectively.
%    The functions thetabxp, thetabyp, thetabxm, thetabym are computed
%    in a separate m-file named thetabz.m.
%    The main complication is deciding how to choose the square when the points
%    lie on the boundary of two squares.
ixm=zeros(zlen,blen);   iym=zeros(zlen,blen);
if Iset == 1 | Iset == 3
   ixp=zeros(zlen,blen);   iyp=zeros(zlen,blen);
elseif Iset == 2
   ixp=zeros(1,1);   iyp=zeros(1,1);
end
%
%  When (thetabx,thetaby) = (thetabxm,thetabym), we only consider bn(l) > 0,
%  not bn(l) >=0.  Note for Iset = 1 or 3, we have replaced bn(l) 
%  by -bn(l) in the formula for theta_b, 
%  so this corresponds to the case Im[b(z)] < 0.
%
		    Lxp = zeros(zlen,blen,2);
                    Lyp = zeros(zlen,blen,2);
for k = 1:zlen
  for l = 1: blen
   if (Iset == 1 | Iset ==3)
    [thetabxpkl,thetabypkl]  = thetabz(1, x(k),y(k),bm(l),bn(l));
    if (thetabxpkl <= .5)
      ixp(k,l) = floor(thetabxpkl/h);
    else
      ixp(k,l) = ceil(thetabxpkl/h)-1;
    end
    if (thetabypkl <= .25)
      iyp(k,l) = floor(thetabypkl/h);
    else
      iyp(k,l) = ceil(thetabypkl/h)-1;
    end
%
%    display([k,l,thetabxpkl,thetabypkl,ixp(k,l),iyp(k,l)])
%   We next compute the Lagrange basis functions Lxp, Lyp, Lxm, and Lym  
%   corresponding to each of these degrees of freedom at the point
%   (thetabx,thetaby).
%
    ixklp = ixp(k,l);     iyklp = iyp(k,l);
     Lxp(k,l,2) = thetabxpkl/h - double(ixklp);
     Lyp(k,l,2) = thetabypkl/h - double(iyklp);
     Lxp(k,l,1) = 1.0 - Lxp(k,l,2);
     Lyp(k,l,1) = 1.0 - Lyp(k,l,2);
   end
       if (bn(l) > 0)
	 [thetabxmkl,thetabymkl]  = thetabz(-1, x(k),y(k),bm(l),bn(l));
       if (thetabxmkl <= .5)  
         ixm(k,l) = floor(thetabxmkl/h);
       else
         ixm(k,l) = ceil(thetabxmkl/h)-1;
       end
      if (thetabymkl <= .25)  
         iym(k,l) = floor(thetabymkl/h);
      else
         iym(k,l) = ceil(thetabymkl/h) -1;
      end
%    display([k,l,  thetabxmkl, thetabymkl, ixm(k,l), iym(k,l)])
    ixklm = ixm(k,l);     iyklm = iym(k,l);
    Lxm(k,l,2) = thetabxmkl/h - double(ixklm);
    Lym(k,l,2) = thetabymkl/h - double(iyklm);
    Lxm(k,l,1) = 1.0 - Lxm(k,l,2);
    Lym(k,l,1) = 1.0 - Lym(k,l,2);
    end
end
end
%
%    display([ixp,iyp,ixm,iym])
%
% Compute the root of log(max evalue(s)) = 0 using the Matlab routine fzero.
%  
% First do the computation for the approximate matrix.
mtype = 0;
funa = @(s) log(eigs(matas(Iset, mtype, s, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym),1));
%%% sa = fzero(funa, [1.5,2.0]);
%%% sa = fzero(funa, [1.2,2.0]);
sa = fzero(funa, 1.5);
[vaa, muaa] = eigs(matas(Iset, mtype, sa, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym),1);
muaa1 = muaa - 1;
disp ('Approximate value of the Hausdorff dimension.')
sa

%
% Next do the computation for the lower matrix.
		    mtype = -1;
funl = @(s) log(eigs(matas(Iset, mtype, s, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym),1));
sl = fzero(funl, sa);
As = matas(Iset, mtype, sl, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym);
[val, mual] = eigs(As,1);
mual1 = mual - 1;
%
% Add loop to ensure mual1 > 0 and min(As*val - val) >=0
% First make sure that val is a strictly positive eigenvector
if (max(val) <= 0) val = -val;
end
while (mual1 <= 0 || min(As*val - val) < 0)
  sl = sl - eps;
  As = matas(Iset, mtype, sl, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym);
  [val, mual] = eigs(As,1);
   mual1 = mual - 1;
if (max(val) <= 0) val = -val;
end
% min(As*val - val)
end
disp ('Lower bound on the Hausdorff dimension.')
sl
% disp ('The eigenvalue of the  lower matrix minus one.')
% disp ('The value should be >=0.')
mual1
% The following value should be >=0
minl = min(As*val - val)
%
% Next do the computation for the upper matrix
mtype = 1;
funu = @(s) log(eigs(matas(Iset, mtype, s, N, J, R, ixp, ixm, iyp, iym, x, y, bm,  bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym),1));
su = fzero(funu, sl);
Bs = matas(Iset, mtype, su, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym);
[vau, muau] = eigs(Bs,1);
muau1 = muau - 1;
%
% Add loop to ensure muau1 <0 and min(vau - Bs*vau) >=0
if (max(vau) <= 0) vau = -vau;
end
while (muau1 >= 0 || min(vau - Bs*vau) < 0)
su = su +  eps;
Bs = matas(Iset, mtype, su, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen, Lxp, Lxm, Lyp, Lym);
[vau, muau] = eigs(Bs,1);
muau1 = muau -1;
if (max(vau) <= 0) vau = -vau;
end
% min(vau - Bs*vau)
end
disp ('Upper bound on the Hausdorff dimension.')
su
% disp ('The eigenvalue of the  upper matrix minus one.')
% disp ('The value should be <=0.')
muau1
% The following value should be >=0
minu = min(vau - Bs*vau)
%       end program hd2di
