function [time, phi, theta, u, v, frequency] = ExactInertialSolution( varargin )
%EXACTINERTIALSOLUTION Exact solution for an inertial oscillation on the Earth.
%Jeffrey J. Early, 2011. http://jeffreyearly.com
%
%   [time,phi,theta,u,v,frequency]=ExactInertialSolution(latitude,u)
%   returns the time, the longitude (phi), the latitude (theta), the
%   velocity (u, v) and the frequency of an inertial particle constrained
%   to the surface of the earth, given initial conditions of latitude
%   (latitude) and the eastward velocity (u0). The duration is
%   automatically set to exactly one inertial period (inclusive), or some
%   reasonable maximum if the inertial period is too long. The resolution
%   is then set to exactly duration/100.
%   
%   [time,phi,theta,u,v,frequency]=ExactInertialSolution(latitude,u,
%   duration) returns the exact inertial solution to the specified
%   duration. The resolution is then set to exactly duration/100.
%   
%   [time,phi,theta,u,v,frequency]=ExactInertialSolution(latitude,u,
%   duration,resolution) returns the exact inertial solution to the
%   specified duration with a given resolution (in seconds).
%   
%   ---------------------------------------------------------------------
%   
%   Angles are in degrees, time is in seconds, speed is in meters per
%   second, and the frequency is in cycles per second.
%   
%   The exactly velocity (u,v)=(R*\dot{phi}*cos(theta_0),R*\dot{theta})
%   is computed directly from analytical functions.
%   
%   ---------------------------------------------------------------------
%   
%   Example usage:
%   Case m < 1
%     [time,phi,theta,u,v,f]=ExactInertialSolution(45,25);
%   
%   Case m = 1
%     [time,phi,theta,u,v,f]=ExactInertialSolution(45,136.2242);
%   
%   Case m > 1
%     [time,phi,theta,u,v,f]=ExactInertialSolution(45,200);

if strcmp(varargin{1}, '--f')
  ExactInertialSolution_fig,return
end

if strcmp(varargin{1}, '--test')
  ExactInertialSolution_test,return
end

if nargin < 2
	fprintf('At least two inputs are required (latitude and initial eastward velocity)\n'), return
else
	latitude = varargin{1};
	if ( latitude < -90.0 || latitude > 90.0 )
		fprintf('Latitude must be between -90 and 90.'), return
	end
	
	u0 = varargin{2};
	if ( u0 < 0.0 )
		fprintf('The velocity must be greater than or equal to zero.'), return
	end
	
	if nargin >= 3
		duration = varargin{3};
		
		if nargin >= 4
			resolution = varargin{4};
		else
			resolution = duration/100;
		end
	else
		% Constrain our automated duration to twice the equatorial inertial period.
		duration = 1/ExactInertialFrequency(latitude,u0);
		maxduration = 2/ExactInertialFrequency(0.001,u0);
		if (duration > maxduration)
			duration = maxduration;
		end
		
		% Lets also set the resolution so that our time intervals are exact fractions of the inertial period.
		resolution = duration/100;
	end;
	
end


time=(0:resolution:duration);

% Constants of the problem
sidereal_day = 86164.1;
earth_radius = 6378100.0;
omega = 2 * pi / sidereal_day;
theta_0 = latitude * pi / 180;
f = 2*omega*sin(theta_0);
beta = 2*omega*cos(theta_0)/earth_radius;

% Modulus and squared modulus of the elliptic functions
m = 4*( u0*beta/(f*f) + u0*u0/(earth_radius*earth_radius*f*f) );
k = sqrt(m);

if abs(1-m) < 1e-2
	sch = sech( omega * sin(theta_0) * time);
	tnh = tanh( omega * sin(theta_0) * time);

	% Path
	theta = asin( sin(theta_0) * sch );
	phi = omega * cos(theta_0) *( time + (1/(omega*cos(theta_0))) *atan( tan(theta_0) * tnh ) ) - omega*time;
	
	% Velocity
	coeff1 = 1 ./ sqrt( 1 + tan(theta_0)*tan(theta_0)* (tnh.*tnh) );
	coeff2 = f*f/(2*beta);
	u = u0 - coeff2 * ( tnh.*tnh );
	v = -coeff2 * ( sch.*tnh );
	u = u .* coeff1;
	v = v .* coeff1;
	
elseif  m < 1
	u_phi = u0 / (earth_radius * cos(theta_0) * omega);
	[sn, cn, dn, am] = EllipticJ( omega * sin(theta_0) * time, m);
	
	% Path
	theta = asin( sin(theta_0) * dn );
	phi = (1 + u_phi) * EllipticPi( am, m, -m*tan(theta_0)*tan(theta_0)) / sin(theta_0) - omega*time;
	
	% Velocity
	coeff1 = 1 ./ sqrt( 1 + m*tan(theta_0)*tan(theta_0)*(sn.*sn) );
	coeff2 = 2*(u0 + u0*u0/(earth_radius*earth_radius*beta));
	u = u0 - coeff2 * ( sn.*sn );
	v = -coeff2 * ( sn.*cn );
	u = u .* coeff1;
	v = v .* coeff1;
elseif m > 1
	u_phi = u0 / (earth_radius * cos(theta_0) * omega);
	[sn, cn, dn, am] = EllipticJ( k * omega * sin(theta_0) * time, 1/m);
	
	% Path
	theta = asin( sin(theta_0) * cn );
	phi = (1 + u_phi) * EllipticPi( am, 1/m, -tan(theta_0)*tan(theta_0)) / (k * sin(theta_0) ) - omega*time;
	
	% Velocity
	coeff1 = 1 ./ sqrt( 1 + tan(theta_0)*tan(theta_0)*(sn.*sn) );
	coeff2 = 2*(u0 + u0*u0/(earth_radius*earth_radius*beta));
	coeff3 = f*f/(2*beta);
	u = u0 - coeff3 * ( sn.*sn );
	v = -coeff3 * k * ( sn.*dn );
	u = u .* coeff1;
	v = v .* coeff1;
end

phi = phi*180/pi;
theta = theta*180/pi;
frequency = ExactInertialFrequency( latitude, u0);

end

function[]=ExactInertialSolution_fig
	frequency1 = ExactInertialFrequency( 45, 75 );
	frequency2 = ExactInertialFrequency( 45, 136.2242 );
	frequency3 = ExactInertialFrequency( 20, 50 );
	
	duration1 = 10/frequency1;
	duration2 = 1/frequency2;
	duration3 = 5/frequency3;
	
	[time1, phi1, theta1, u1, v1, f1] = ExactInertialSolution( 45, 75, duration1, duration1/200 );
	phi1 = phi1 + 220;
	[time2, phi2, theta2, u2, v2, f2] = ExactInertialSolution( 45, 136.2242 );
	phi2 = phi2 + 250;
	[time3, phi3, theta3, u3, v3, f3] = ExactInertialSolution( 20, 50, duration3, duration3/200 );
	phi3 = phi3+120;
	
	
	radius = 1.001;
	
	load('topo.mat','topo','topomap1');
	
	figure;
	[x,y,z] = sphere(50);

	cla reset
	axis square off
	props.AmbientStrength = 0.1;
	props.DiffuseStrength = 1;
	props.SpecularColorReflectance = .5;
	props.SpecularExponent = 20;
	props.SpecularStrength = 1;
	props.FaceColor= 'texture';
	props.EdgeColor = 'none';
	props.FaceLighting = 'phong';
	props.Cdata = topo;
	surface(x,y,z,props);
	
	colormap(topomap1);
	light('position',[-1 0 1]);
	light('position',[-1.5 0.5 -0.5], 'color', [.6 .2 .2]);
	
	
	hold on;
	
	plot3( radius*cosd(phi1).*cosd(theta1), radius*sind(phi1).*cosd(theta1), radius*sind(theta1), 'Color','red','LineWidth',2 )
	
	plot3( radius*cosd(phi2).*cosd(theta2), radius*sind(phi2).*cosd(theta2), radius*sind(theta2), 'Color','red','LineWidth',2 )
	
	plot3( radius*cosd(phi3).*cosd(theta3), radius*sind(phi3).*cosd(theta3), radius*sind(theta3), 'Color','red','LineWidth',2 )
	
	hold off;
	
	axis equal
	view(-72,12)
	zoom(1.5)
end

function[]=ExactInertialSolution_test

% Constants of the problem
sidereal_day = 86164.1;
earth_radius = 6378100.0;
omega = 2 * pi / sidereal_day;

frequency1 = ExactInertialFrequency( 45, 75 );
frequency2 = ExactInertialFrequency( 45, 136.2242 );
frequency3 = ExactInertialFrequency( 20, 50 );

duration1 = 10/frequency1;
duration2 = 1/frequency2;
duration3 = 5/frequency3;

% The three possible solutions
[time1, phi1, theta1, u1, v1, f1] = ExactInertialSolution( 45, 75, duration1, duration1/200 );
[time2, phi2, theta2, u2, v2, f2] = ExactInertialSolution( 45, 136.2242 );
[time3, phi3, theta3, u3, v3, f3] = ExactInertialSolution( 20, 50, duration3, duration3/200 );
	
% Is the energy conserved? (yeah, yeah, we're looking at the speed, but that's okay).
energy1 = sqrt( u1.*u1 + v1.*v1 );
energy2 = sqrt( u2.*u2 + v2.*v2 );
energy3 = sqrt( u3.*u3 + v3.*v3 );

% How about the angular momentum?
L1 = earth_radius * cosd( theta1 ) .* ( u1 + omega*earth_radius*cosd(theta1));
L2 = earth_radius * cosd( theta2 ) .* ( u2 + omega*earth_radius*cosd(theta2));
L3 = earth_radius * cosd( theta3 ) .* ( u3 + omega*earth_radius*cosd(theta3));

figure

subplot(2,1,1)
plot( time1, energy1, 'b');
hold on
plot( time2, energy2, 'r');
plot( time3, energy3, 'k');
xlabel('Time (s)')
ylabel('Speed (m/s)')

subplot(2,1,2)
plot( time1, L1, 'b');
hold on
plot( time2, L2, 'r');
plot( time3, L3, 'k');
xlabel('Time (s)')
ylabel('Angular momentum (m^2/s)')

end


function q = EllipticPi( phi, m, c )

q = zeros(1, length(phi) );

fHandle = @(u)1./((1 - c.*sin(u).^2).*sqrt(1 - m.*sin(u).^2));

for n=1:length(phi)
	q(n) = quadgk(fHandle, 0, phi(n),'RelTol',1e-8,'AbsTol',1e-12);
end

end

function [sn,cn,dn,am] = EllipticJ(u,m,tol)
%ELLIPJ Jacobi elliptic functions and Jacobi's amplitude.
%   [Sn,Cn,Dn,Am] = ELLIPJ(U,M) returns the values of the Jacobi
%   elliptic functions SN, CN, DN and AM evaluated for corresponding
%   elements of argument U and parameter M.  The arrays U and M must
%   be the same size (or either can be scalar).  As currently
%   implemented, M is limited to 0 <= M <= 1. 
%
%   [Sn,Cn,Dn,Am] = ELLIPJ(U,M,TOL) computes the elliptic functions to
%   the accuracy TOL instead of the default TOL = EPS.  
%
%   Some definitions of the Jacobi elliptic functions use the modulus
%   k instead of the parameter m.  They are related by m = k^2.
%
%   See also ELLIPKE.

%   L. Shure 1-9-88
%   Copyright 1984-2001 The MathWorks, Inc. 
%   $Revision: 5.14 $  $Date: 2001/04/15 12:01:40 $
%
%   Modified by Moiseev Igor, 
%   moiseev[at]sissa.it
%   34106, SISSA, via Beirut n. 2-4,  Trieste, Italy
%   Date: 2005/10/04
%
%   The modification of orginal script fixes the problem of slow convergence
%   of the AGM algorithm for the value of parameter M=1.


%   ELLIPJ uses the method of the arithmetic-geometric mean
%   described in [1].
%
%   References:
%   [1] M. Abramowitz and I.A. Stegun, "Handbook of Mathematical
%       Functions" Dover Publications", 1965, Ch. 16-17.6.


if nargin<3, tol = eps; end
if nargin<2, error('Not enough input arguments.'); end

if ~isreal(u) | ~isreal(m)
    error('Input arguments must be real.')
end

if length(m)==1, m = m(ones(size(u))); end
if length(u)==1, u = u(ones(size(m))); end
if ~isequal(size(m),size(u)), error('U and M must be the same size.'); end

am = zeros(size(u));
cn = zeros(size(u));
sn = cn;
dn = sn;
m = m(:).';    % make a row vector
u = u(:).';

if any(m < 0) | any(m > 1), 
  error('M must be in the range 0 <= M <= 1.');
end

I = uint32( find(m ~= 1 & m ~= 0) );
if ~isempty(I)
    [mu,J,K] = unique(m(I));   % extracts unique values from m
    K = uint32(K);
    mumax = length(mu);

    % pre-allocate space and augment if needed
	chunk = 7;
	a = zeros(chunk,mumax);
	c = a; 
	b = a;
	a(1,:) = ones(1,mumax);
	c(1,:) = sqrt(mu);
	b(1,:) = sqrt(1-mu);
	n = uint32( zeros(1,mumax) );
	i = 1;
	while any(abs(c(i,:)) > tol)                                    % Arithmetic-Geometric Mean of A, B and C
        i = i + 1;
        if i > size(a,1)
          a = [a; zeros(2,mumax)];
          b = [b; zeros(2,mumax)];
          c = [c; zeros(2,mumax)];
        end
        a(i,:) = 0.5 * (a(i-1,:) + b(i-1,:));
        b(i,:) = sqrt(a(i-1,:) .* b(i-1,:));
        c(i,:) = 0.5 * (a(i-1,:) - b(i-1,:));
        in = uint32( find((abs(c(i,:)) <= tol) & (abs(c(i-1,:)) > tol)) );
        if ~isempty(in)
          [mi,ni] = size(in);
          n(in) = ones(mi,ni)*(i-1);
        end
	end

    mmax = length(I);
	phin = zeros(1,mmax);
	phin(:) = (2 .^ double(n(K))).*a(i,K).*u(I);
	while i > 1
        i = i - 1;
        in = uint32( find(n(K) >= i) );
        if ~isempty(in)
          phin(in) = 0.5*(asin(c(i+1,K(in)).*sin(phin(in))./a(i+1,K(in))) + phin(in));
        end
	end
	am(I) = phin;
	sn(I) = sin(phin);
	cn(I) = cos(phin);
	dn(I) = sqrt(1 - m(I).*sin(phin).^2);
end

% Special cases: m = {0, 1} 
m0 = find(m == 0);
am(m0) = u(m0);
sn(m0) = sin(u(m0));
cn(m0) = cos(u(m0));
dn(m0) = ones(size(m0));

m1 = find(m == 1);
am(m1) = asin(tanh(u(m1)));
sn(m1) = tanh(u(m1));
cn(m1) = sech(u(m1));
dn(m1) = sech(u(m1));
end