(*
Solution to Problem 3 of the Bigger Problem Set - ABC flow.
At the Mathematica prompt, type
<< abc.math
*)

(* Time to integrate to, starting from t = 0 *)
tp = 10;

(* The A, B, C constants *)
a = 1;
b = 1;
c = 0;

(* Differential equations to solve, starting at x0,y0,z0 *)
deq[x0_,y0_,z0_] :=
	{x'[t] == b Cos[y[t]] + c Sin[z[t]],
	 y'[t] == a Sin[x[t]] + c Cos[z[t]],
	 z'[t] == a Cos[x[t]] + b Sin[y[t]],
	 x[0] == x0, y[0] == y0, z[0] == z0}
	
(* Solution of the diffeq *)
soln[x0_,y0_,z0_] := soln[x0,y0,z0] \
	= NDSolve[deq[x0,y0,z0], {x,y,z}, {t,0,tp}];

(* Solutions for x, y, z, starting at x0,y0,z0 *)
x[t_,x0_,y0_,z0_] := x[t] /. soln[x0,y0,z0];
y[t_,x0_,y0_,z0_] := y[t] /. soln[x0,y0,z0];
z[t_,x0_,y0_,z0_] := z[t] /. soln[x0,y0,z0];

(* Procedure to plot the evolution of x, y, z *)
plot[x0_,y0_,z0_] := Plot[{
	x[t,x0,y0,z0],
	y[t,x0,y0,z0],
	z[t,x0,y0,z0]},
	{t,0,tp}
];

(* Procedure to plot the evolution of the separation between two points *)
plot[x0_,y0_,z0_,dx0_,dy0_,dz0_] := Plot[{
	x[t,x0+dx0,y0+dy0,z0+dz0] - x[t,x0,y0,z0],
	y[t,x0+dx0,y0+dy0,z0+dz0] - y[t,x0,y0,z0],
	z[t,x0+dx0,y0+dy0,z0+dz0] - z[t,x0,y0,z0]},
	{t,0,tp}
];

(* Procedure to plot the evolution of the |separation| between two points *)
plotabs[x0_,y0_,z0_,dx0_,dy0_,dz0_] := Plot[
	 (x[t,x0+dx0,y0+dy0,z0+dz0] - x[t,x0,y0,z0])^2
	+(y[t,x0+dx0,y0+dy0,z0+dz0] - y[t,x0,y0,z0])^2
	+(z[t,x0+dx0,y0+dy0,z0+dz0] - z[t,x0,y0,z0])^2,
	{t,0,tp}
];

(* Procedure to plot flow lines in the x-y plane *)
paraplot[ic__] := ParametricPlot[ic, {t,0,tp}];

(* Table of initial conditions *)
initconds[t_] :=
	Table[Flatten[{x[t,x0,Pi/2,0], y[t,x0,Pi/2,0]}], {x0,0,2 Pi,Pi/12}];

(* Plot flow lines in the x-y plane *)
paraplot[initconds[t]]
	
(* Plot evolution of x, y, z *)
(* plot[0,2,0] *)

(* Plot the evolution of the separation between two points *)
(* plot[0,2,0,0,.2,0] *)

(* Plot the evolution of the |separation| between two points *)
(* plotabs[0,2,0,0,.2,0] *)