function [p,t] = bendp(filename, ax, ay)

% BENDP bends initial mesh element saved in FILENAME over 
% a parabola with focal length AX (XZ plane) and/or AY (YZ plane).
% AX = AY means bending over a paraboloid surface.

% Usage:  bendp('plate.mat', AX, AY) or
%           bendp('plate', AX, AY) 


%Load the data.
load(filename);

% Bend the mesh over a parabola where theta1 and theta2 correspond 
%to the angles in XZ and YZ plane respectively.

if (ax == 0 | ay == 0)
    error('Focal length has to be a non-zero number.');
else
    theta1 = atan(p(1,:)/ax);
    p(1,:) = 2*ax*tan(theta1/2);
    rho = ax*(sec(theta1/2)).^2;
    p(3,:) = p(3,:) + ax- sqrt(rho.^2 - (p(1,:)).^2);
    
    theta2 = atan(p(2,:)/ay);
    p(2,:) = 2*ay*tan(theta2/2);
    rho = ay*(sec(theta2/2)).^2;
    p(3,:) = p(3,:) + ay - sqrt(rho.^2 - (p(2,:)).^2);
end


%Save and view the result.
save bendp p t;
viewer bendp