(Satellite on orbit) - Create the tangent cone from point to sphere

Question

I am trying to draw a cone, connected to the sphere in Matlab. I have the point [x1,y1,z1] outside of the sphere [x2,y2,z2] with R radius and I want it to be the top of the cone, created out of tangents.

On this pictures you can see what I have in mind:

Below you can see what I have already done. I am using it in order to mark the part of the Earth's surface, visible from the satellite position in orbit. Unfortunately, the cone in this picture is approximate, I need to create accurate one, connected with surface. For now, it is not only inaccurate, but also goes under it.

I am creating the sphere with this simple code (I am skipping the part of putting the map on it, it's just an image):

r = 6371.0087714;
[X,Y,Z] = sphere(50);
X2 = X * r;
Y2 = Y * r;
Z2 = Z * r;
surf(X2,Y2,Z2)
props.FaceColor= 'texture';
props.EdgeColor = 'none';
props.FaceLighting = 'phong';
figure();
globe = surface(X2,Y2,Z2,props);

Let's assume that I have the single point in 3D:

plot3(0,0,7000,'o');

How can I create such a cone?

saastn · Answer 1 · 2020-11-23T21:36:41.940

There are two different questions here:

How to calculate cone dimensions?
How to plot lateral faces of a 3D cone?

Calculating Cone Dimensions

Assuming that center of sphere is located on [0 0 0]:

d = sqrt(Ax^2+Ay^2+Az^2);
l = sqrt(d^2-rs^2);
t = asin(rs/d);
h = l * cos(t);
rc = l * sin(t);

Plotting the Cone

The following function returns coordinates of lateral faces of cone with give apex point, axis direction, base radius and height, and the number of lateral faces.

function [X, Y, Z] = cone3(A, V, r, h, n)
% A: apex, [x y z]
% V: axis direction, [x y z]
% r: radius, scalar
% h: height, scalar
% n: number of lateral surfaces, integer
% X, Y, Z: coordinates of lateral points of the cone, all (n+1) by 2. You draw the sphere with surf(X,Y,Z) or mesh(X,Y,Z)
v1 = V./norm(V);
B = h*v1+A;
v23 = null(v1);
th = linspace(0, 2*pi, n+1);
P = r*(v23(:,1)*cos(th)+v23(:,2)*sin(th));
P = bsxfun(@plus, P', B);
zr = zeros(n+1, 1);
X = [A(1)+zr P(:, 1)];
Y = [A(2)+zr P(:, 2)];
Z = [A(3)+zr P(:, 3)];
end

The Results

rs = 6371.0087714; % globe radius
A = rs * 2 * [1 1 1]; % sattelite location
V = -A; % vector from sat to the globe center
% calculating cone dimensions
d = norm(A); % distance from cone apex to sphere center
l = (d^2-rs^2)^.5; % length of generating line of cone
sint = rs/d; % sine of half of apperture
cost = l/d; % cosine of half of apperture
h = l * cost; % cone height
rc = l * sint; % cone radius

% globe surface points
[XS,YS,ZS] = sphere(32);
% cone surface points
[XC, YC, ZC] = cone3(A, V, rc, h, 32);

% plotting results
hold on
surf(XS*rs,YS*rs,ZS*rs, 'facecolor', 'b', 'facealpha', 0.5, 'edgealpha', 0.5)
surf(XC, YC, ZC, 'facecolor', 'r', 'facealpha', 0.5, 'edgealpha', 0.5);

axis equal
grid on

Animating the satellite

The simplest way to animate objects is to clear the whole figure by clf and plot objects again in new positions. But a way more efficient method is to plot all objects once and in each frame, only update positioning data of moving objects:

clc; close all; clc
rs = 6371.0087714; % globe radius
r = rs * 1.2;
n = 121;
t = linspace(0, 2*pi, n)';
% point on orbit
Ai = [r.*cos(t) r.*sin(t) zeros(n, 1)];

[XS,YS,ZS] = sphere(32);
surf(XS*rs,YS*rs,ZS*rs, 'facecolor', 'b', 'facealpha', 0.5, 'edgealpha', 0.5)
hold on
[XC, YC, ZC] = cone3(Ai(1, :), Ai(1, :), 1, 1, 32);
% plot a cone and store handel of surf object
hS = surf(XC, YC, ZC, 'facecolor', 'r', 'facealpha', 0.5, 'edgealpha', 0.5);

for i=1:n
    % calculating new point coordinates of cone
    A = Ai(i, :);
    V = -A;
    d = norm(A);
    l = (d^2-rs^2)^.5;
    sint = rs/d;
    cost = l/d;
    h = l * cost;
    rc = l * sint;
    
    [XC, YC, ZC] = cone3(A, V, rc, h, 32);
    % updating surf object
    set(hS, 'xdata', XC, 'ydata', YC, 'zdata', ZC);
    pause(0.01); % wait 0.01 seconds
    drawnow(); % repaint figure
end

Another sample with 3 orbiting satellites:

Thank you! Why is there `l = (d^2-rs^2)^.5;` Why to you need this `.^5` ? — KS0232, Nov 23 '20 at 17:13
@KamilSerafin That's the Pythagorean theorem: `d^2 = rs^2+l^2`. — saastn, Nov 23 '20 at 17:25
@KamilSerafin maybe you are confused with the syntax, that's powered by 0.5, not 5. The zero can be eliminated in most languages. — saastn, Nov 23 '20 at 17:28
Right, I can see now! Do you have any idea how can I delete the created cone? I am creating animation, and after each iteration I have to create the new cone and delete the previous one from figure. I've tried using `delete [XC,YC,ZC]` but old cone is still in there (after a few seconds I have hundreds of them). How can I destroy previous cone? — KS0232, Nov 23 '20 at 18:37
@KamilSerafin actually I made an animation the other day after posting this. I cleaned up my code and added an example. — saastn, Nov 23 '20 at 21:38

(Satellite on orbit) - Create the tangent cone from point to sphere

1 Answers1