%       program hd2dPnse
%       Current version: August 11, 2017
%       Previous version: October 4, 2016
%       Computes Hausdorff dimension of a special example where
%       the exact solution is known.
%         Let z = x + iy
%         Let D:= {(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 theta_b(z) = 1/(z+b).
%         Let I = \{1 \pm i, 2 \pm i, 3 \pm i\}
%       To compute the Hausdorff dimension of K, 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, Im(b) >=0}
%		   [g_b(z)]^s v_s(1/(\bar z + \bar b))
%                + \sum_{b in I_j, Im(b) < 0}
%		 [g_b(z)]^s   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 polynomial of degree deg in each
%       variable separately.
%   
%      Uses the routines eval2dPnse.m and thetabz.m
%
%    Key variables are described below.
%    Iset = integer (1,2,or 3) corresponding to I_1, I_2, or I_3     
%    zlen = number of vertices in the rectangular mesh
%    zlent = 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.
%    aa is the approx matrix
%    muaa is the largest eigenvalue of aa
%  
%      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
Iset = 3;
prompt = 'Input N '
N = input(prompt)
%
%      N = 200
%      N = 100
%      N = 50
%
h = 1.0/double(N)
%     deg = the degree (in each variabel separately)
%     of the piecewise polynomial approximation
      prompt = 'Input degree of piecewise polynomial ';  
      deg = input(prompt)
N2 = N/2;
N21 = N2-1;
%
%
%    Construct a domain D_h made up of a union of squares of side h
%    such that D_h contains D.  D_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 D_h.  Although we could simply
%    choose D_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 D_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 D_h that contains D.
%
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 D_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 D_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 = sum(2*J + 1);
% display(tf)
% display(J)
display(zlen)
Jr = double(J);
%
%    Next, compute the mesh points in D_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 
%  In order to define the piecewise polynomial approximation, 
%  in addition to the values at the corners of each mesh rectangle,
%  which we term "mesh points,"
%  we also need  values of the approximating function at 
%  (deg -1) interpolation points on each mesh rectangle edge and 
%  (deg -1) x (deg -1) additional interior points in each mesh rectangle.
%  First compute the total number of interpolation points 
%  in our approximate vector
%  
%  Recall, for each y = jh, 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.
%  So the number of interpolation points at level y=jh is 2*deg*J(j) + 1.
%  For y = jh, j = 2, ... N/2-1, and deg >1, we add
%  2*deg*J(j) + 1 interpolation points at levels
%  (j -i/deg)*h, i = 1, ... deg -1.
%  Then the total number of interpolation points is:
		    zlent =  deg*sum(2*deg*J +1) - (deg-1)*(2*deg*J(1) +1)
%    Compute all these points in D_h, indexing them
%    by a single index k of length zlent starting at 1.
%    The index runs over the points in each row from left to right
%    starting with the points at level y=0.
%    Compute the points so they are symmetric around x = 1/2.
%
xt=zeros(zlent,1);     yt=zeros(zlent,1);
ll=1;
for jl = 1:N2+1
  jlr = double(jl);
  jlJ = 2*J(jl);
  up = min(deg*(N2+1-jl), deg-1);
  for jj= 0:up
    for i = 0:J(jl)-1
    ir = double(i);
      for ii = 0:deg-1
      xt(ll + deg*i + ii) = .5 - J(jl)*h + ir*h + ii*h/deg;
      xt(ll + deg*2*J(jl) - deg*i - ii) = .5 + J(jl)*h - ir*h - ii*h/deg;
      yt(ll + deg*i+ ii) = (jlr-1)*h + jj*h/deg;
      yt(ll + deg*2*J(jl) - deg*i - ii) = (jlr-1)*h +jj*h/deg;
      end
    end
    xt(ll + deg*J(jl)) = .5;
    yt(ll + deg*J(jl)) = (jlr-1)*h + jj*h/deg;
		    if (jj ==0)   ll = ll + deg*2*J(jl) +1;
                    else ll = ll + deg*2*J(jl+1) +1;
                    end
  end
end 

%    display([xt, yt])
% 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 D and mesh domain D_h')
 axis('equal')
 axis([-.1, 1.1, -.1, .6])
 hold on
 scatter(x,y,5,'s','fill')
 hold off % reset hold state
   blen = 3;
   bm = zeros(blen,1);  bn= zeros(blen,1);
   %    These are
   bm(1) = 3;    bn(1) = 1;
   bm(2) = 2;    bn(2) = 1;
   bm(3) = 1;    bn(3) = 1;
   bmr = zeros(blen,1);
   bnr = zeros(blen,1);
   for k  = 1:blen
      bmr(k) = double(bm(k));
      bnr(k) = double(bn(k));
   end
%
%    For each of the interpolation points, compute the functions
%    aphibp = (1/6)|1/(z +1)|^2 |(z+b+1)/(z+b)|^2
%   = (1/6) (1/((xt + 1)^2 + yt^2)((xt+m+1)^2 + (yt+1)^2)/((xt+m)^2 + (yt+1)^2)
%    aphibm = (1/6)|1/(z +1)|^2 |(z+ \bar b+1)/(z+ \bar b)|^2
%   = (1/6) (1/((xt + 1)^2 + yt^2)((xt+m+1)^2 + (yt-1)^2)/((xt+m)^2 + (yt-1)^2)
aphibp= zeros(zlent,blen);
aphibm= zeros(zlent,blen);

for k = 1:zlent
  for l = 1:blen
	     aphibp(k,l) = (1/6)*(1/((xt(k) + 1)^2 + yt(k)^2)) ...
	  *((xt(k)+bmr(l)+1)^2 + (yt(k)+1)^2)/((xt(k)+bmr(l))^2 + (yt(k)+1)^2);
	     aphibm(k,l) =  (1/6)*(1/((xt(k) + 1)^2 + yt(k)^2)) ...
     *((xt(k) + bmr(l)+1)^2 + (yt(k)-1)^2)/((xt(k)+bmr(l))^2 + (yt(k)-1)^2);
   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.
ixp=zeros(zlent,blen);   iyp=zeros(zlent,blen);
   ixm=zeros(zlent,blen);   iym=zeros(zlent,blen);
%
%  When (thetabx,thetaby) = (thetabxm,thetabym), we only consider bn(l) > 0,
%  not bn(l) >=0.  Note we have replaced bn(l) by -bn(l) in the formula
%  for theta_b, so this corresponds to the case Im[b(z)] < 0
%
for k = 1:zlent
  for l = 1:blen
    [thetabxpkl,thetabypkl]  = thetabz(1, xt(k),yt(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);
    for ideg = 1:deg+1
     Lxp(k,l,ideg) = 1.;
     Lyp(k,l,ideg) = 1.;
        for isum = 0:deg
           if (isum ~= ideg-1)
     Lxp(k,l,ideg) = Lxp(k,l,ideg)*(thetabxpkl/h - double(ixklp) - isum/deg);
     Lyp(k,l,ideg) = Lyp(k,l,ideg)*(thetabypkl/h - double(iyklp) - isum/deg);
           end
        end
     denom = factorial(ideg-1)*factorial(deg-ideg+1)*(-1)^(deg-ideg+1)/deg^deg;
        Lxp(k,l,ideg) = Lxp(k,l,ideg)/denom;
        Lyp(k,l,ideg) = Lyp(k,l,ideg)/denom;
    end
       if (bn(l) > 0)  & (Iset == 1 | Iset ==3)
	 [thetabxmkl,thetabymkl]  = thetabz(-1, xt(k),yt(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);
    for ideg = 1:deg+1
     Lxm(k,l,ideg) = 1.;
     Lym(k,l,ideg) = 1.;
        for isum = 0:deg
           if (isum ~= ideg-1)
     Lxm(k,l,ideg) = Lxm(k,l,ideg)*(thetabxmkl/h - double(ixklm) - isum/deg);
     Lym(k,l,ideg) = Lym(k,l,ideg)*(thetabymkl/h - double(iyklm) - isum/deg);
           end
        end
     denom = factorial(ideg-1)*factorial(deg-ideg+1)*(-1)^(deg-ideg+1)/deg^deg;
        Lxm(k,l,ideg) = Lxm(k,l,ideg)/denom;
        Lym(k,l,ideg) = Lym(k,l,ideg)/denom;
      end
    end
end
end
%
%    display([ixp,iyp,ixm,iym])
%
% Compute the root of log(max evalue(s)) = 0 using the Matlab routine fzero.
%  
fun = @(s) log(eval2dPnse(Iset, s, N, J, ixp,ixm, iyp, iym, aphibp, ...
aphibm, zlent, blen,deg, Lxp, Lxm, Lyp, Lym, bn));
 s = fzero(fun, 1.5);
mua = eval2dPnse(Iset, s, N, J, ixp,ixm, iyp, iym, aphibp, ...
aphibm, zlent, blen,deg, Lxp, Lxm, Lyp, Lym, bn);
mua1 = mua - 1;
s
%
%
%       end program hd2dPnse
