Matlab: Punkt auf ellipsoidalen Oberflächen

(German) This post is focussing on the solution of calculating a point of an ellipsoidal surface in dependency of to angles like in spherical coordinates.

Mode of operation:

This function calculates the closest point on an ellipsoidal surface from an point outside or inside the body in respect to its center. It needs the mesh data of the ellipsoidal body X, Y, Z and the point, where the point with direct connection to the center should be found on the surface. The output parameters are [point, phi, the] = PointOnSurfaceAnalytic(X, Y, Z, x, y, z); the surface point and the angles phi and theta of the ellipsoidal equations. It can be downloaded here.  

To derive the formula we start with the ellipsoidal equation:
(x/a)^2 + (y/b)^2 + (z/c)^2 =1 
If we consider 𝝋 to be between the x and y axis and 
𝝑 be between the y and z axis it leads to:
x=y/tan 𝝋  
z=y/tan 𝝑  
Substitute this in the ellipsoidal equation one gets:
(y/(a tan 𝝋))^2 + (y/b)^2 + (y/(tan 𝝑c))^2 =1  
y^2* [(1/(a tan 𝝋))^2 + (1/b)^2 + (1/(c tan 𝝑))^2] =1  
solve this for y gives:
y= [1/(1/(a tan 𝝋)^2 + (1/b)^2 + 1/(c tan 𝝑)^2)]^1/2 
so one can calculate x, y, z and therefor one gets the coordinates of the point on the surface for a given 𝝋 and 𝝑.

Source code:

function [p, varargout] = PointOnSurfaceAnalytic(X, Y, Z, x, y, z)
% This function calculates the nearest point on the given surface [X, Y, Z]
% from the point [x, y, z]. The surface are mesh data with dimension NxN.

[x_mid, a] = mat_max(X, 'mid', 'dist');
[y_mid, b] = mat_max(Y, 'mid', 'dist');
[z_mid, c] = mat_max(Z, 'mid', 'dist');

phi = atan2(y-y_mid, x-x_mid);
the = atan2(z-z_mid, y-y_mid);

y_surf1 = ( 1/(1/(tan(the)^2/c^2 + 1/b^2 + 1/(tan(phi)^2*a^2) ) )^(1/2);
x_surf1 = y_surf1/tan(phi);
z_surf1 = y_surf1*tan(the);

y_surf2 = -( 1
/(1/(tan(the)^2/c^2 + 1/b^2 + 1/(tan(phi)^2*a^2) ) )^(1/2);
x_surf2 = y_surf2/tan(phi);
z_surf2 = y_surf2*tan(the);

p1 = [x_surf1, y_surf1, z_surf1];
p2 = [x_surf2, y_surf2, z_surf2];
pd = [x, y, z];

d1 = sum((p1-pd).^2);
d2 = sum((p2-pd).^2);

if d1 < d2
    p = p1;
else
    p = p2;
end

varargout{1} = phi;
varargout{2} = the;

end


%% ______________Max Matrix___________________%%
function [varargout] = mat_max(mat, varargin)
% calculates certain values of a nxn mesh matrix (mat) like min, max value,
% the mid value, the sice of the (ellipsoidal) matrix

varargout = cell(length(varargin),1);

logic = strcmp(varargin,'mid');
if any(logic)~= 0
    v_min = min(min(mat));
    v_max = max(max(mat));
    pos = find(logic,1);
    varargout{pos} = (v_min+v_max)/2;
end

logic = strcmp(varargin,'dist');
if any(logic)~= 0
    v_min = min(min(mat));
    v_max = max(max(mat));
    pos = find(logic,1);
    varargout{pos} = (v_max-v_min)/2;
end

end

Kommentare

Beliebte Posts aus diesem Blog

Matlab: Fehlergeraden

Matlab: 3D Kartesisches Grid, Vektoren mit Transformation/Rotation

Matlab: Farbspektrum für Plots