

%phlft5 component

function U=biharmNormal(Dev)
%               b2
%  ^          +----+
% m|       b1 |    | b3 
%  + - >      +----+
%    n          b4


[m,n] = size(Dev);

dx = 1/(n-1);
dy = 1/(m-1);
cory =dy:dy:1-dy;
corx =dx:dx:1-dx;

[x,y] = meshgrid(corx,cory);



b1 = Dev(2:end-1,1);
b2 = Dev(end,2:end-1);
b3 = Dev(2:end-1,end);
b4 = Dev(1,2:end-1);

yScale = (cory-1).*cory;
xScale = (corx-1).*corx;
% Now for the actual "solution" of the poisson equation
db1 = dst(b1./yScale');
db2 = dst(b2./xScale);
db3 = dst(b3./yScale');
db4 = dst(b4./xScale);

lambdan = pi*[1:n-2];
lambdam = pi*[1:m-2]';

u = repmat(lambdam,1,n-2);
v = repmat(lambdan,m-2,1);



%h3 = (x-1).*x.*sinh(u.*x)./sinh(u);
h3 = (x-1).*x.*exp(u.*(x-1)) .* ( 1 - exp(-2*u.*x))./(1-exp(-2*u));
%h1 = x.*(x-1).*sinh(u.*(1-x))./sinh(u);
h1 = -(x-1).*x.*exp(u.*(-x)) .* ( 1 - exp(-2*u.*(1-x)))./(1-exp(-2*u));
%h2 = (y-1).*y.*sinh(v.*y)./sinh(v);
h2 = (y-1).*y.*exp(v.*(y-1)) .* ( 1 - exp(-2*v.*y))./(1-exp(-2*v));
%h4 = y.*(y-1).*sinh(u.*(1-y))./sinh(u);
h4 = -(y-1).*y.*exp(v.*(-y)) .* ( 1 - exp(-2*v.*(1-y)))./(1-exp(-2*v));

h1 = repmat(db1,1,n-2).*h1;
h3 = repmat(db3,1,n-2).*h3;
h2 = repmat(db2,m-2,1).*h2;
h4 = repmat(db4,m-2,1).*h4;
U=zeros(m,n);
U(2:end-1,2:end-1) = idst(h1+h3).*repmat(yScale',1,n-2) + idst(h2'+h4')'.*repmat(xScale,m-2,1);

