Optimal trajectory over the surface of a sphere

Question 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;

Question 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

Question 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.

Question 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

Question 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;
}