%  Daniel Kawano, Rose-Hulman Institute of Technology
%  Last modified:  Feb 23, 2016

clear all
close all
clc

%  Specify symbolic parameters and variables that are functions of the
%  directed curve's convected coordinate xi:

syms psi(xi) theta(xi) phi(xi) 
syms nu1(xi) nu2(xi) nu3(xi)

%  (1)  Relate the reference triad {D1,D2,D3} = {E1,E2,E3} to the triad
%       {d1,d2,d3} using a 3-2-3 set of Euler angles:

P1 = [cos(psi), sin(psi), 0; 
      -sin(psi), cos(psi), 0; 
      0, 0, 1];

P2 = [cos(theta), 0, -sin(theta);
      0, 1, 0;
      sin(theta), 0, cos(theta)];  
  
P3 = [cos(phi), sin(phi), 0; 
      -sin(phi), cos(phi), 0; 
      0, 0, 1];

P = simplify(P3*P2*P1);

%  (2)  Compute the axial vector of strains, nu, which yields a set of 
%       ODEs governing the evolution of the Euler angles for the triad
%       {d1,d2,d3} along the directed curve:

Ptensor = transpose(P);

Ktensor = simplify(transpose(Ptensor)*diff(Ptensor, xi));

Ktensor32 = [0, 0, 1]*Ktensor*[0, 1, 0]';
Ktensor13 = [1, 0, 0]*Ktensor*[0, 0, 1]';
Ktensor21 = [0, 1, 0]*Ktensor*[1, 0, 0]';

nu = [Ktensor32; Ktensor13; Ktensor21];

%  (3)  Manipulate the system of ODEs into a form suitable for numerical 
%       integration:

%  Compile the state equations and arrange the state variables:
     
StateEqns = [nu1; nu2; nu3] == nu;

StateVars = [psi; theta; phi];

%  Express the state equations in mass-matrix form, M(xi,Y)*Y'(xi) = 
%  F(xi,Y):

[Msym, Fsym] = massMatrixForm(StateEqns, StateVars);

Msym = simplify(Msym)
Fsym = simplify(Fsym)

%  Convert M(xi,Y) and F(xi,Y) to symbolic function handles with the input 
%  parameters specified:

M = odeFunction(Msym, StateVars, nu1(xi), nu2(xi), nu3(xi));  
F = odeFunction(Fsym, StateVars, nu1(xi), nu2(xi), nu3(xi));  

%  Save M(xi,Y) and F(xi,Y):

save directed_curve_ODEs.mat Msym Fsym StateVars M F