function as = matas(Iset, mtype, s, N, J, R, ixp, ixm, iyp, iym, x, y, bm, bn, aphibp, aphibm, zlen, blen,  Lxp, Lxm, Lyp, Lym)

%
  h = 1.0/double(N);
  gam = min(bm);
%    We compute either the lower matrix (mtype = -1), the approximate
%    matrix (mtype = 0), or the upper matrix (mtype = 1) of size (zlen, zlen)
%
%    For Iset = 1 or 3, 
%    the computation is done by finding matrices corresponding to
%    values of b = m + in first for which n >=0 (the matrix ap),
%    and then the matrix for which n <0, (the matrix am), and 
%    finally adding the two.  Only the matrix am is needed for Iset = 2.
%    The main issue here is to sum over the values of b which give 
%    contributions from smallest to largest, to avoid loss of accuracy. 
%
%  We use the facts that for all z_1, z_2 in D,
%  v_s(z_1) <= v_s(z_2) exp(sqrt{5}*s |z_1 -z_2|),
%  -s/[4(s+1]*v_s(z) <= D_{xx} v_s(z) <= 2*s*(2*s +1)*v_s(z),
%  -2s v_s(z) <= D_{yy} v_s(z) <= 2*s*(2*s +1)/4*v_s(z).
%
% If Theta_b(z) \in R_{k,l} with lower left corner x_k,y_l, define the 
% bilinear interpolant v_s^I of v_s by:
% v_s^I(\theta_b(z)) 
%  = (x_{k+1} - Re[theta_b(z)])*(y_{l+1} - Im[theta_b(z)])/h^2*v_s(z_{k,l})
%  + (Re[theta_b(z)]-x_k)*(y_{l+1} - Im[theta_b(z)])/h^2*v_s(z_{k+1,l})
%  + (x_{k+1} - Re[theta_b(z)])*(Im[theta_b(z)]- y_l)/h^2*v_s(z_{k,l+1})
%  + (Re[theta_b(z)]-x_k)*(Im[theta_b(z)]- y_l)/h^2*v_s(z_{k+1,l+1})
%
%  From the error for bilinear interpolation, we have for
%      x_k \le Re[theta_b(z)] \le x_{k+1}},
%      y_l \le Im[theta_b(z)] \le y_{l+1},
% 
%  [1 - err_b^1(z_{i,j})] v_s^I(theta_b(z_{i,j})) <= v_s(theta_b(z_{i,j})) 
%  <= [1 + err_b^@(z_{i,j})] v_s^I(theta_b(z_{i,j}))
%  where
%  err_b^1(theta_b(z)) = 
%  [(x_{k+1} - Re[theta_b(z)])*(Re[theta_b(z)]-x_k)
%  + (y_{l+1} - Im[theta_b(z)])*(Im[theta_b(z)]-y_l)]
%  *s(2s+1) exp(sqrt{10}sh/gam)/gam^2
%  err_b^2(theta_b(z)) = 
%  [(x_{k+1} - Re[theta_b(z)])*(Re[theta_b(z)]-x_k)
%  + (y_{l+1} - Im[theta_b(z)])*(Im[theta_b(z)]-y_l)]*s(9+8s)/(8+8s)
%  *exp(sqrt{10}sh/gam)/gam^2
%    Initialize matrices
as = zeros(zlen, zlen);
ap = zeros(zlen, zlen);
am = zeros(zlen, zlen);
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));
    ixklp = ixp(k,l);     iyklp = iyp(k,l);
    w1p = Lyp(k,l,2);    w0p = Lyp(k,l,1);
    v1p =  Lxp(k,l,2);   v0p = Lxp(k,l,1);
%   display([v0p, v1p, w0p, w1p])
%
%   The following calculations use the function kindex = kindx(il,jl,N,J)
%   residing in a matlab file of the same name. 
%   Suppose the mesh points are ordered in each row from left to right
%   starting with the mesh points at level y=0.
%   Given a mesh point of the form (il*h,jl*h),
%   with il = N/2-J(jl+1), N/2 + J(jl+1) and 
%   jl =  0, N/2, we wish to determine the index k of that mesh point 
%   in an ordered list beginning with 1.
%
    kk00p = kindx(ixklp,iyklp,N,J);
    kk10p = kindx(ixklp +1,iyklp,N,J);
    kk01p = kindx(ixklp,iyklp+1,N,J);
    kk11p = kindx(ixklp+1,iyklp+1,N,J);      
%  
%   The following statements check that we have not exceeded the array size
%
     if (kk00p > zlen) display(['kk00p > zlen', kk00p, zlen, ixklp, iyklp])
     end
     if (kk10p > zlen) display(['kk10p > zlen', kk10p, zlen, ixklp, iyklp])
     end
     if (kk01p > zlen) display(['kk01p > zlen', kk01p, zlen, ixklp, iyklp])
     end
     if (kk11p > zlen) display(['kk11p > zlen', kk11p, zlen, ixklp, iyklp])
     end
%     display([kk00p,kk10p,kk01p,kk11p])
%     display(kk11p)
%
% 
if (mtype == -1)
   err = -h^2*s*(2*s+1)*exp(s*sqrt(10.0)*h/gam)*(v0p*v1p + w0p*w1p)/gam^2;
elseif (mtype == 0)
   err = 0.;
elseif (mtype ==1)
err = h^2*s*((9+8*s)/(8+8*s))*exp(s*sqrt(10.0)*h/gam)*(v0p*v1p + w0p*w1p)/gam^2;
else
  fprintf('mtype is not one of the allowed values');  
end
%
    aphibpskl = aphibp(k,l)^s;
    ap(k,kk00p) = ap(k,kk00p) + aphibpskl*v0p*w0p*(1.0+err);
    ap(k,kk10p) = ap(k,kk10p) + aphibpskl*v1p*w0p*(1.0+err);
    ap(k,kk01p) = ap(k,kk01p) + aphibpskl*v0p*w1p*(1.0+err);
    ap(k,kk11p) = ap(k,kk11p) + aphibpskl*v1p*w1p*(1.0+err);
end
    if (bn(l) > 0)
%     display(bn(l))
      [thetabxmkl,thetabymkl]  = thetabz(-1, x(k),y(k),bm(l),bn(l));
      ixklm = ixm(k,l);     iyklm = iym(k,l);
      w1m = thetabymkl/h - double(iyklm);   w0m = 1.0 - w1m;
      v1m = thetabxmkl/h - double(ixklm);   v0m = 1.0 - v1m;
%       display([v0m, v1m, w0m, w1m])
      kk00m = kindx(ixklm,iyklm,N,J);
      kk10m = kindx(ixklm +1,iyklm,N,J);
      kk01m = kindx(ixklm,iyklm+1,N,J);
      kk11m = kindx(ixklm+1,iyklm+1,N,J);       
       if (kk00m > zlen) display(['kk00m > zlen', kk00m, zlen, ixklm, iyklm])
       end
       if (kk10m > zlen) display(['kk10m > zlen', kk10m, zlen, ixklm, iyklm])
       end
       if (kk01m > zlen) display(['kk01m > zlen', kk01m, zlen, ixklm, iyklm])
       end
       if (kk11m > zlen) display(['kk11m > zlen', kk11m, zlen, ixklm, iyklm])
       end
%     display(kk11m)
if (mtype == -1)
   err = -h^2*s*(2*s+1)*exp(s*sqrt(10.0)*h/gam)*(v0m*v1m + w0m*w1m)/gam^2;
elseif (mtype == 0)
   err = 0.;
elseif (mtype ==1)
err = h^2*s*((9+8*s)/(8+8*s))*exp(s*sqrt(10.0)*h/gam)*(v0m*v1m + w0m*w1m)/gam^2;
else
  fprintf('mtype is not one of the allowed values');  
 end
      aphibmskl = aphibm(k,l)^s;
      am(k,kk00m) = am(k,kk00m) + aphibmskl*v0m*w0m*(1.0+err);
      am(k,kk10m) = am(k,kk10m) + aphibmskl*v1m*w0m*(1.0+err);
      am(k,kk01m) = am(k,kk01m) + aphibmskl*v0m*w1m*(1.0+err);
      am(k,kk11m) = am(k,kk11m) + aphibmskl*v1m*w1m*(1.0+err);
    end
%       display([v0m, v1m, w0m, w1m])
  end
end
    if (Iset == 1 | Iset == 3)
      as = ap + am;
    else as = am;
    end

%
% Corrections to the matrices for |b| > R
% 
if (mtype == 1)
   if (Iset == 1)
      delRs = exp(s/sqrt(R^2-R))*(R/(R-1))^s*((1/(2*s-1))*(1/(R-1))^(2*s-1) ...
      + (pi/2)*(1/(s-1))*(1/(R - sqrt(2)))^(2*s-2));
    elseif  (Iset == 2)
      delRs = exp(s/sqrt(R^2-R))*(R/(R-1))^s* ...
       (pi/4)*(1/(s-1))*(1/(R - sqrt(2)))^(2*s-2);
    elseif (Iset == 3)
      delRs = 0;
    else
       fprintf('Iset is not one of the allowed values');
     end
    for ll = 1:zlen
       as(ll,1) = as(ll,1) + delRs;
    end
end
if (mtype == -1)
   if (Iset == 1)
     etaRs = exp(-sqrt(5)*s/sqrt(R^2-R))*(R/(R+sqrt(5)+5/(4*R)))^s ...
       *(pi - 2*asin(1/(R +sqrt(2))))*(1/(2*s-2))*(1/(R+sqrt(2)))^(2*s-2);
    elseif  (Iset == 2)
      etaRs = exp(-sqrt(5)*s/sqrt(R^2-R))*(R/(R+sqrt(5)+5/(4*R)))^s ...
      *(pi/2 - 2*asin(1/(R +sqrt(2))))*(1/(2*s-2))*(1/(R+sqrt(2)))^(2*s-2);
    elseif (Iset == 3)
      etaRs = 0;
    else
       fprintf('Iset is not one of the allowed values');
    end
    for ll = 1: zlen
       as(ll,1) = as(ll,1) + etaRs;
    end
end

% Computation of the largest eigenvalue of the matrix a and corresponding
% eigenvector by the Matlab routine eigs.  
%    [va,mua] = eigs(a,1);
% end of function matas


