Optimal trajectory over the surface of a sphere

Question

I am trying to find the parametric equation of the trajectory of a point jumping over different spots on the surface of a unit sphere such that:

1. each jump is small (pi/4 < d < pi/2) and in a narrow interval, e.g. [1.33, 1.34]
2. the point visits most regions of the sphere as quickly and uniformly as possible
3. the point travels along "direction vectors" as different as possible

This is what I have tried

N = 3600;    % number of points
t = (1:N) * pi / 180;    % parameter
theta_sph = sqrt(2) * t * pi;    % first angle
phi_sph = sqrt(3) * t * pi;    % second angle
rho_sph = 1;    % radius
% Coordinates of a point on the surface of a sphere
x_sph = rho_sph * sin(phi_sph) .* cos(theta_sph);
y_sph = rho_sph * sin(phi_sph) .* sin(theta_sph);
z_sph = rho_sph * cos(phi_sph);

% Check length of jumps (it is intended that this is valid only for small jumps!!!)
aa = [x_sph(1:(N-1)); y_sph(1:(N-1)); z_sph(1:(N-1))];
bb = [x_sph(2:N); y_sph(2:N); z_sph(2:N)];
cc = cross(aa, bb);
d = rho_sph * atan2(arrayfun(@(n) norm(cc(:, n)), 1:size(cc,2)), dot(aa, bb));
figure
plot(d, '.')
figure
plot(diff(d), '.')

% Check trajectory on the surface of the sphere
figure
hh = 1;
h_plot3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), '-');
hold on
axis square
% axis off
set(gca, 'XLim', [-1 1])
set(gca, 'YLim', [-1 1])
set(gca, 'ZLim', [-1 1])
for hh = 1:N
  h_point3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), ...
      'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
  drawnow
  delete(h_point3)
  set(h_plot3, 'XData', x_sph(1:hh))
  set(h_plot3, 'YData', y_sph(1:hh))
  set(h_plot3, 'ZData', z_sph(1:hh))
end

EDIT --> Can anybody find a more regular trajectory, perhaps covering the sphere more quickly (i.e. with the smallest number of jumps) and more uniformly? Regular trajectory in the sense that it should change direction smoothly, not sharply. Aesthetic beauty is a bonus. The points should be spread on the surface of the sphere as uniformly as possible.

Mucky

Can you further clarify the criteria? Specifically points 2 and 3, and also what you mean by 'regular' and 'uniform' at the end? – Durman

Uniform means that the points should be spread all over the surface of the sphere, the farthest from each other (this can be measured alright). Improperly, of course, just as a generic intuition: by regular trajectory I meant a trajectory that doesn't change too much, i.e. it should change direction very smoothly, not sharply. And if it is aesthetically pleasing, that's bonus :-) – Mucky

Can you clarify point 3? If I ignore it, I expect a spiral down the sphere (which is the path you are using) to be optimal. Just adjust the number of full rotations the polar angle goes through to control jump size and coverage. If you control the polar rate you can 'stagger' your points so the don't all fall on the same latitude lines since it seems you are more concerned about the 'point' coverage than the continuous trajectory. – Heidyheifer

The reason why I don't like the case that I have illustrated is that the point keeps visiting the north and the south poles, i.e. at the two poles the concentration of points is way greater than anywhere else. Can the spiral be changed accordingly? How would you adjust the number of full rotations; how would you control the jump size; and how would you control the polar rate? Can you show an example? – Mucky

Answer 1

Modifying the beginning of your code to introduce a rotation of the underlying sphere. This gives a trajectory that does not return to the poles as frequently. It may take some tuning of the rotation speeds to look "nice" (and it probably looks better when it is only rotating around one axis, not all 3). rot_angle1 is rotation around the x-axis, and rot_angle2 and rot_angle3 are rotation around the y and z axes. Maybe this gives you an idea at least!

N = 3600;    % number of points
t = (1:N) * pi / 180;    % parameter
theta_sph = sqrt(2) * t * pi;    % first angle
phi_sph = sqrt(3) * t * pi;    % second angle
rho_sph = 1;    % radius
rot_angle1 = sqrt(2) * t * pi;
rot_angle2 = sqrt(2.5) * t * pi;
rot_angle3 = sqrt(3) * t * pi;
% Coordinates of a point on the surface of a sphere
x_sph0 = rho_sph * sin(phi_sph) .* cos(theta_sph);
y_sph0 = rho_sph * sin(phi_sph) .* sin(theta_sph);
z_sph0 = rho_sph * cos(phi_sph);

x_sph1 = x_sph0;
y_sph1 = y_sph0.*cos(rot_angle1)-z_sph0.*sin(rot_angle1);
z_sph1 = y_sph0.*sin(rot_angle1)+z_sph0.*cos(rot_angle1);

x_sph2 = x_sph1.*cos(rot_angle2)+z_sph1.*sin(rot_angle2);
y_sph2 = y_sph1;
z_sph2 = -x_sph1.*sin(rot_angle2)+z_sph1.*cos(rot_angle2);

x_sph = x_sph2.*cos(rot_angle3)-y_sph2.*sin(rot_angle3);
y_sph = x_sph2.*sin(rot_angle3)+y_sph2.*cos(rot_angle3);
z_sph = z_sph2;

Tiga

Answer 2

I don't have a copy of matlab handy, but I'll post the modifications I would make to your curve.

Just to be clear since there are n-finity+1 different definitions for spherical angles. I will use the following, its backwards from your definition but I'm bound to make mistakes if I try and switch.

• \phi - the angle from the z-axis
• \theta - the projected angle in the x-y plane.

The Parameterization

Let t be a discrete set of N evenly spaced points from 0 to pi (inclusive).

\phi(t) = t
\theta = 2 * c * t

Rather straight forward and simple, a spiral around a sphere is linear in \phi and theta. c is a constant that represents the number of full rotations in \theta, it need not be integer.

Neighboring Points

In your example, you calculate the angle between vectors with atan2(norm(cross....) which is fine, but doesn't give any insight into the problem. Your problem is on the surface of a sphere, use this fact. So I consider the distance between points using this formula

Now you find neighboring points, these occur at t +- dt and theta +- 2pi for whatever t this happens.

In the first case t +- dt, it is easy to calculate cos(gamma) = 1 - 2 c^2 sin^2(t) dt^2. The sin^2(t) dependence is why the poles are more dense. Ideally you want to choose theta(t) and phi(t) such that dtheta^2 * sin^2(phi) is constant and minimal to satisfy this case.

The second case is a little more difficult and brings up my comments about "staggering" your points. If we choose an N such that dtheta does not evenly divide 2pi, then after a full rotation around the sphere in theta I cannot end up directly beneath a previous point. To compare the distance between points in this case, use delta t so that c delta t = 1. Then you have delta phi = delta t and delta theta = 2 c delta t - 2pi. Depending on your choice of c, delta phi may or may not be small enough to use the small angle approximations.

Final notes

It should be apparent that c=0 is a straight line down the sphere. By increasing c you increase the "density of the spiral" gaining more coverage. However, you also increase the distance between neighboring points. You will want to choose a c for the chosen N that makes the two distance formulas above approximately equal.

EDIT moved some things to mathbin for cleanliness

Heidyheifer

I have posted a modified version of my original code - do you think it incorporates your hints, or is there anything that I have overlooked? To compile in Matlab/Octave you can use compileonline.com/execute_matlab_online.php – Mucky

I have also added the formula that you suggested for the distance between two points. – Mucky

Answer 3

I edit here because it's a long code. After the hints from David and kalhartt, I have tried this:

N = 3600;    % number of points
t = (1:N) * pi / 180;    % parameter
% theta_sph much faster than phi_sph to avoid overly visiting the poles
theta_sph = sqrt(20.01) * t * pi;    % first angle
phi_sph = sqrt(.02) * t * pi;    % second angle
rho_sph = 1;    % radius
% Coordinates of a point on the surface of a sphere
x_sph0 = rho_sph * sin(phi_sph) .* cos(theta_sph);
y_sph0 = rho_sph * sin(phi_sph) .* sin(theta_sph);
z_sph0 = rho_sph * cos(phi_sph);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use David's hint to rotate axes (but only the first now):
rot_angle1 = t * pi / 10;
rot_angle2 = 0;
rot_angle3 = 0;

x_sph1 = x_sph0;
y_sph1 = y_sph0.*cos(rot_angle1)-z_sph0.*sin(rot_angle1);
z_sph1 = y_sph0.*sin(rot_angle1)+z_sph0.*cos(rot_angle1);

x_sph2 = x_sph1.*cos(rot_angle2)+z_sph1.*sin(rot_angle2);
y_sph2 = y_sph1;
z_sph2 = -x_sph1.*sin(rot_angle2)+z_sph1.*cos(rot_angle2);

x_sph = x_sph2.*cos(rot_angle3)-y_sph2.*sin(rot_angle3);
y_sph = x_sph2.*sin(rot_angle3)+y_sph2.*cos(rot_angle3);
z_sph = z_sph2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check length of jumps
aa = [x_sph(1:(N-1)); y_sph(1:(N-1)); z_sph(1:(N-1))];
bb = [x_sph(2:N); y_sph(2:N); z_sph(2:N)];
cc = cross(aa, bb);
d = rho_sph * atan2(arrayfun(@(n) norm(cc(:, n)), 1:size(cc,2)), dot(aa, bb));
figure
plot(d, '.')

% Check trajectory on the surface of the sphere
figure
hh = 1;
h_plot3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), '-');
hold on
axis square
% axis off
set(gca, 'XLim', [-1 1])
set(gca, 'YLim', [-1 1])
set(gca, 'ZLim', [-1 1])
for hh = 1:N
  h_point3 = plot3(x_sph(hh), y_sph(hh), z_sph(hh), ...
      'o', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r');
  drawnow
  delete(h_point3)
  set(h_plot3, 'XData', x_sph(1:hh))
  set(h_plot3, 'YData', y_sph(1:hh))
  set(h_plot3, 'ZData', z_sph(1:hh))
end

I think it's much better than before! Two things that I have found to be important: 1. theta_sph must be much faster than phi_sph to avoid visiting two opposite poles too often; 2. if theta_sph goes faster than phi_sph, then you have to slowly rotate over either rot_angle1 or rot_angle2 in order to obtain a trajectory that is not too messy. I'm still open to any other hints to improve the result.

Mucky

Answer 4

clear all
close all

u = pi/2:(-pi/36):-pi/2;
v = 0:pi/36:2*pi;

nv = length(v);
nu = length(u);

f = myfigure(1);
ax = myaxes(f,1,1);

hold on

for aa = 1:nv
    tv = v(aa);
    for bb = 1:nu
        tu = u(bb);
        x = sin(tu)*cos(tv);
        y = cos(tu)*cos(tv);
        z = sin(tv);
        plot3(x,y,z,'*')
    end
end

edit: myfigure and myaxes are functions I have for creating a figure and an axes

Revelatory

Answer 5

I've written up a quick version in C that behaves very well given a fixed number of points. You can play with it at ideone. If you have WebGL enabled browser (Chrome, Firefox) you can paste those results here to see them plotted. The poles are a little bit off due to some integral approximations used in deriving the formulas, but other than that it's hard to see a flaw. There's no constants that should need tweaking other than the number of points you desire output for.

#include <stdio.h>
#include <math.h>

int main(void) {
    int i, numPoints = 200;
    double slope = sqrt(1.2) / sqrt(numPoints);
    for (i = 0; i < numPoints; i++) {
        double s = asin((double)i / (double)(numPoints - 1) * 2.0 - 1.0);
        double z = sin(s);
        double r = sqrt(1.0 - z * z);
        double ss = (2.0 * s + M_PI) / slope;
        double x = cos(ss) * r;
        double y = sin(ss) * r;
        printf("%lf,%lf,%lf,\"1\"\n", x, y, z);
    }
    return 0;
}

Mell