   % calculate the value of Eigenfunction at the surface
function Zsurf=cal_Zsurf(cn,W,P,H,ssp,N,h)
    kn=W/cn;
    Z0=cal_z0(W,P,kn,ssp(1),H);
    Z1=cal_z1(Z0,W,P,kn,ssp(1),ssp(2),H,h);
    Z=cal_z(Z0,Z1,W,kn,ssp,H,N);
    Zsurf=real(Z(N+1));       % the real part            
    