function fz = linear_beltrami_solver_rect(face,z,mu,corner,landmark,target,s)
nv = size(z,1);
bd = compute_bd(face);
A = generalized_laplacian(face,[real(z),imag(z)],mu);
b = complex(zeros(size(z)));

b(landmark) = target;
b(corner) = z(corner);
landmark = [landmark;corner];
Al = A(landmark,:);
ind = length(landmark)*(landmark-1)+(1:length(landmark))';
Al(ind) = Al(ind) - 1;
A(landmark,:) = A(landmark,:) - Al;

i = find(bd==corner(1),1,'first');
bd = bd([i:end,1:i]);
corner = corner([1:end,1]);
zc = z(corner);

if exist('s','var')
    b = real(b)*s(1)+1i*imag(b)*s(2);
    zc = real(zc)*s(1)+1i*imag(zc)*s(2);
end

ck = zeros(size(corner));
k = 1;
for i = 1:length(bd)
    if(bd(i) == corner(k))
        ck(k) = i;
        k = k+1;
    end
end
nl = length(bd)-length(corner);
I = zeros(nl,1);
Vx = zeros(nl,1);
Vy = zeros(nl,1);
d = zeros(size(z));
l = 1;
for k = 1:length(ck)-1
    z1 = zc(k);
    z2 = zc(k+1);
    for i = ck(k)+1:ck(k+1)-1
        I(l) = bd(i);
        Vx(l) = imag(z2)-imag(z1);
        Vy(l) = real(z1)-real(z2);
        d(bd(i)) = real(z1)*imag(z2) - imag(z1)*real(z2);
        l = l+1;
    end
end

B = sparse(I,I,Vx,nv,nv);
C = sparse(I,I,Vy,nv,nv);

[I,~,Vx] = find(B);
Ax = A;
bx = real(b);
bx(I) = d(I);
Al = Ax(I,:);
ind = length(I)*(I-1)+(1:length(I))';
Al(ind) = Al(ind) - Vx;
Ax(I,:) = Ax(I,:) - Al;
x = Ax\bx;

[I,~,Vy] = find(C);
Ay = A;
by = imag(b);
by(I) = d(I);
Al = Ay(I,:);
ind = length(I)*(I-1)+(1:length(I))';
Al(ind) = Al(ind) - Vy;
Ay(I,:) = Ay(I,:) - Al;
y = Ay\by;

if exist('s','var')
    fz = x+1i*y;
else
    dxy = sum(x(face).*(z(face(:,[2 3 1]))-z(face(:,[3 1 2]))),2);
    ux = imag(dxy);
    uy = -real(dxy);
    af = (1-2*real(mu)+abs(mu).^2)./(1.0-abs(mu).^2);
    bf = -2*imag(mu)./(1.0-abs(mu).^2);
    cf = (1+2*real(mu)+abs(mu).^2)./(1.0-abs(mu).^2);
    fa = face_area(face,z);
    h = 0.25*sum((af.*ux.^2 + 2*bf.*ux.*uy + cf.*uy.^2)./fa)/sum(fa);
    fz = x+h*1i*y; 
end