function [p,t] = paraboloid(filename, focus, a)

% PARABOLOID projects the initial mesh element saved in 
% FILENAME over the surface of a paraboloid 
% Coordinates of the focus and the focal length are  
% specified by the user.

% Usage:  paraboloid('plate.mat', [0 0 -2], a) or
%         paraboloid('plate', [0 0 -2], a)


% load the data.
load (filename);

if(a==0 | a==inf) 
    error('The focal length cannot be zero or infinity.');
elseif(focus(1,1)==inf | focus(1,2)==inf | focus(1,3)==inf | focus(:,:)==0)
    error('Focal Point of the paraboloid cannot be at origin or at infinity.')
else
    % Get the center of the mesh at origin.
    p(1,:)=p(1,:)-mean(p(1,:));
    p(2,:)=p(2,:)-mean(p(2,:));
    p(3,:)=p(3,:)-mean(p(3,:));
    
    % Project the mesh over the surface of a paraboloid. 
    Q=p-repmat(focus',[1,length(p)]);
    n=sqrt(sum(Q.*Q,1));
    X = (focus + [0 0 a])';
    cos_theta = (dot(Q, repmat(X, [1,length(Q)])))./(n*a);
    rho = (2*a)./(1+cos_theta);
    p=(Q./repmat(n,[3,1])).*repmat(rho,[3,1]);
end

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