
%% ================================================================
% John Kerl
% kerl.john.r@gmail.com
% 2007-02-08
%
% A Matlab script to display the two-sphere along with a tangent plane.
%% ================================================================

%% ----------------------------------------------------------------
ntheta=20;
nphi=20;
ns = 20;
nt = 20;

thetamin=0;
thetamax=pi/2;%2*pi;

phimin=0;
phimax=pi/2;%pi;

smin = -.3;
smax =  .3;

tmin = -.3;
tmax =  .3;

dtheta=(thetamax-thetamin)/ntheta;
dphi=(phimax-phimin)/nphi;
ds=(smax-smin)/ns;
dt=(tmax-tmin)/nt;

%% ----------------------------------------------------------------
clf;         % Clear the figure in case this is a redisplay
rotate3d on; % Permit me to rotate the graph with click and drag
hold on;     % All plot3 output goes on the same graph

%% ----------------------------------------------------------------
% Draw latitude lines
for phi = phimin : dphi : phimax
	theta = [ thetamin : dtheta : thetamax ];
	x = cos(theta) * sin(phi);
	y = sin(theta) * sin(phi);
	z = cos(phi) * ones(1, ntheta+1);
	plot3(x,y,z,'g')
end

% Draw longitude lines
for theta = thetamin : dtheta : thetamax
	phi = [ phimin : dphi : phimax ];
	x = cos(theta) * sin(phi);
	y = sin(theta) * sin(phi);
	z = cos(phi);
	plot3(x,y,z,'g')
end

% Draw the equator
phi = pi/2;
theta = [ thetamin : dtheta : thetamax ];
x = cos(theta) * sin(phi);
y = sin(theta) * sin(phi);
z = cos(phi) * ones(1, ntheta+1);
plot3(x,y,z,'r')

% Draw the Greenwich meridian
theta = 0
phi = [ phimin : dphi : phimax ];
x = cos(theta) * sin(phi);
y = sin(theta) * sin(phi);
z = cos(phi);
plot3(x,y,z,'k')

% Matlab line-style codes:
%           b     blue          .     point              -     solid
%           g     green         o     circle             :     dotted
%           r     red           x     x-mark             -.    dashdot 
%           c     cyan          +     plus               --    dashed   
%           m     magenta       *     star             (none)  no line
%           y     yellow        s     square
%           k     black

%% ----------------------------------------------------------------
% Point q is (1/2, 1/2, 1/sqrt 2).  This has spherical coordinates
% theta = 45 deg and phi = 45 deg.
qx = 1/2;
qy = 1/2;
qz = 1/sqrt(2);

q = [qx; qy; qz];

qtheta = pi/4;
qphi   = pi/4;

%% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Tangent vector 1 is d/dx using graph coordinates
% Tangent vector 2 is d/dy using graph coordinates

%  q =
% [ x ]
% [ y ]
% [ (1 - x^2 - y^2)^1/2 ]

% d/dx = 
% [ 1 ]
% [ 0 ]
% [ -x/z ]

% d/dy = 
% [ 0 ]
% [ 1 ]
% [ -y/z ]

ugx = 1;
ugy = 0;
ugz = -qx/qz;

vgx = 0;
vgy = 1;
vgz = -qy/qz;

%% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Tangent vector 1 is d/dtheta using spherical coordinates
% Tangent vector 2 is d/dphi   using spherical coordinates

%  q =
% [ cos(theta) sin(phi)]
% [ sin(theta) sin(phi)]
% [ cos(phi)]

% d/dtheta =
% [-sin(theta) sin(phi)]
% [ cos(theta) sin(phi)]
% [ 0]

% d/dphi =
% [ cos(theta) cos(phi)]
% [ sin(theta) cos(phi)]
% [-sin(phi)]

% d/dtheta at q =
% [-1/2]
% [ 1/2]
% [ 0]

% d/dphi at q =
% [ 1/2 ]
% [ 1/2 ]
% [-1/sqrt(2)]

usx = -sin(qtheta) * sin(qphi);
usy =  cos(qtheta) * sin(qphi);
usz =  0;

vsx =  cos(qtheta) * cos(qphi);
vsy =  sin(qtheta) * cos(qphi);
vsz = -sin(qphi);

ug = [ugx; ugy; ugz];
vg = [vgx; vgy; vgz];
us = [usx; usy; usz];
vs = [vsx; vsy; vsz];

%q
%u
%v
%u' * v
%q' * u
%q' * v


% Draw t lines using graph coordinates
for s = smin : ds : smax
	t = [ tmin : dt : tmax ];
	x = qx + s * ugx + t * vgx;
	y = qy + s * ugy + t * vgy;
	z = qz + s * ugz + t * vgz;
	plot3(x,y,z,'b')
end

% Draw s lines using graph coordinates
for t = tmin : dt : tmax
	s = [ smin : ds : smax ];
	x = qx + s * ugx + t * vgx;
	y = qy + s * ugy + t * vgy;
	z = qz + s * ugz + t * vgz;
	plot3(x,y,z,'b')
end

% Draw t lines using spherical coordinates
for s = smin : ds : smax
	t = [ tmin : dt : tmax ];
	x = qx + s * usx + t * vsx;
	y = qy + s * usy + t * vsy;
	z = qz + s * usz + t * vsz;
	plot3(x,y,z,'k')
end

% Draw s lines using spherical coordinates
for t = tmin : dt : tmax
	s = [ smin : ds : smax ];
	x = qx + s * usx + t * vsx;
	y = qy + s * usy + t * vsy;
	z = qz + s * usz + t * vsz;
	plot3(x,y,z,'k')
end

