function [ecc,inc,w_p,capom,Lsp,obl,index] = orbelm_secular(t,index);
%% HEADER: orbelm_secular.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%Description: Returns orbital elements for planet defined by index
%              using secular theory outlined in Murray and Dermot and Ward
%              2004/1974.
%
%%Inputs:
%       t: vector of yrs since current epoch 
%       index: vector of index positions for desired orbital elements
%              valid range: (1 <= index <= 8)
%                         1: Mercury
%                         2: Venus
%                         3: Earth
%                         4: Asteroid Belt Barycenter
%                         5: Jupiter
%                         6: Saturn (default)
%                         7: Neptune
%                         8: Uranus
%
%%Outputs:
%       ecc: eccentricty
%       inc: incidence angle [radians]
%       w_p: longitude of perihelion [radians]
%       capom: longitude of ascending node [radians]
%       Lsp: Ls of perihelion [radians] (Only Valid For Saturn)
%       obl: Obliquity of Saturn [radians] (ONly Valid For Saturn)
%       index: input index array
%
%%Calls: None
%
%%Required Files: None
%
%%Source:
%        Murray and Dermot, Capter 7, Ward 2004 (Tilting Saturn I and II)
%
%%Author: Alexander G. Hayes (hayes@gps.caltech.edu)
%
%%Notes: Should be updated with newer secular constants thant published in 
%        Murray and Dermot!
%
%%Change Log: 
%            09/13/08: Draft Code Escalated to V1
%            01/10/09: Modified Oblquity Calculations and updated secular 
%                      constants for obliquity ONLY!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% BEGIN CODE: orbelm_secular.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%inputs
%t: yrs from J2000
%index: index of desired values (defaults to Saturn [6]);

%%outputs
%ecc: eccentricy
%inc: orbital inclination to ecliptic
%w_p: longitude of perihelion
%capom: longitude of acsending node
%Lsp: Ls of perihelion (Saturn Only)
%obl: Obliquity (Saturn Only)
%index: vector of valid output index values
%
%NOTE: Lsp and obl only valid for Saturn!!!!

%define the variables to avoid any upleasent warnings
ecc = [];
inc = [];
w_p = [];
capom = [];
Lsp = [];
obl = [];

%make sure t is in correct array format
if numel(t) ~= size(t,2); t = reshape(t,1,numel(t)); end;

%default for index is to return Saturn's values
if ~exist('index','var');
    index = 6;
elseif length(index)==0;
    index = 6;
else;
    index = index( (index >=1) & (index <=8) );
end;
if length(index)==0;
    disp(['All Values of Index are Out of Range! (1<=INDEX<=8)']);
    return;
end;

%define eij vector (pg. 303, Murray and Dermot)
eij = [17454   610   392    64    -1   -1     0    0;...%1
       -2550  2094  1634   290    -1   -1     0    0;...%2
         154 -1255  1043  2972     0   -1     0    0;...%3
        -169  1483 -1483  7237     0   -1     0    0;...%4
        3571  1909  1834  1861  4482 3291 -3188   63;...%5
          96     6   283  1502 -1535 4863  -248  -11;...%6
          43    42    44    59   147  148  3043 -323;...%7
           1     1     1     2     5    6   143  938;...%8
         -21   -11   -12   -20   -16  -49    20    0;...%9
          10     4    15    64   -54  202   -16    0];  %10
eij = eij./1e5;

%define Iij vector (pg. 304, Murray and Dermot)
Iij = [12449  1178  849   180   -2   -3    2    0;...%1
       -3545  1004  810   180   -1   -2    1    0;...%2
         409 -2680 2448 -3589    0    0    0    0;...%3
         116  -685  453  5025    0   -2    0    0;...%4
        2575  2575 2757  2757 2757 2757 2757 2757;...%5
          28    12  281   965 -631 1572  -69   -8;...%6
        -332  -192 -173  -126  -96  -78 1760 -207;...%7
        -145  -132 -130  -123 -117 -113  110 1175];  %8
Iij = Iij./1e5;
    
%define secular eigenfrequencies
%gi fi Bi Gi
eigen = [ 5.46326  -5.20154  92.182  19.433;...
          7.34474  -6.57080 196.881 318.057;...
         17.32832 -18.74359 335.224 255.031;...
         18.00233 -17.63331 317.948 296.541;...
          4.29591   0.00000  29.550 107.102;...
         27.77406 -25.73355 125.120 127.367;...
          2.71931  -2.90266 131.944 315.063;...
          0.63332  -0.67752  69.021 202.293;...
        -19.18225   0.00000 293.979 0.00000;...
         51.25222   0.00000 220.691 0.00000];
eigen(:,1:2) = eigen(:,1:2)./3600.*pi./180; %deg
eigen(:,3:4) = eigen(:,3:4).*pi./180;%deg

%calculate equations at motions
hj = (eij'*sin(repmat(eigen(:,1),1,length(t)).*repmat(t,size(eigen,1),1) + ...
    repmat(eigen(:,3),1,length(t))) );
kj = (eij'*cos(repmat(eigen(:,1),1,length(t)).*repmat(t,size(eigen,1),1) + ...
    repmat(eigen(:,3),1,length(t))) );
pj = (Iij'*sin(repmat(eigen(1:8,2),1,length(t)).*repmat(t,8,1) + ...
    repmat(eigen(1:8,4),1,length(t))) );
qj = (Iij'*cos(repmat(eigen(1:8,2),1,length(t)).*repmat(t,8,1) + ...
    repmat(eigen(1:8,4),1,length(t))) );

%calculate the eccentricity
ecc = sqrt( hj.*hj + kj.*kj );

%calculate the inclination
inc = sqrt( pj.*pj + qj.*qj );

%calculate the longitude of perihelion
w_p = atan2(hj,kj);

%calculate the longitude of the ascending node
capom = atan2(pj,qj);

pos_saturn = find( index == 6);
if length(pos_saturn) > 0;
   %calculate the obliquity (WARD 2004, HAMILTON, 2004, WARD 1974, BRETAGNON 1973)
    obl_J2000 = 26.73;
    %define required constants
    n =  [];%Heliocentric mean motion
    w=  [] ;%Spin Frequency of Saturn
    J2 = 1.6297e-2; %coefficient of quadrapole moment of gravity field
    q = 0.05164; %(French 1993)
    %Effective Quadrupole Moment of satellite sytem
    %q/J2 = ratio of soalr torque on the satellites to that directly
    %exterted on planet
    l = 0.00278; %(French 1993)
    lambda = 0.2199; %(Hubbard and Marley (1989)
    %Moment of Inertia of Saturn normalized to MsR^2
    %define precessional constant (WARD 2004 Eq (1))
    %alpha = 3/2.*n.^2./w.*(J2+q)./(lambda+l);
    alpha = 0.8306/3600; %deg/year
    %projected alpha onto mean orbit plane
    alpha_cosobj = alpha.*cos(obl_J2000.*pi./180);%0.7380./3600;%0.7427./3600; %deg/year
    %large amplitude terms for saturns orbit (WARD 2004, TABLE 1 from BRETAGNON 1974)
    gj = [-26.34,-2.99,-0.692]./3600; %deg/yr
    Ij = [0.910,0.045,0.064]; %deg
    %define phase offset (BRETAGNON 1973)
    delta_j = [125.383423,316.173594,201.171559];%deg
    %define initial pole azimuth
    az_o = 59.9392-90; %59.120 from Ward, 2004, 0.8192 from JPL SPICE analysis
    %Projection of Spin Vector onto orbit plane
    %Clockwise angle from vector perpendicular
    %   to line of ascending node
    % Instantaneous Value at J2000
    %define initial longitude of ascending node
    omega_o = 113.6669;
    %obl (WARD 2004 EQ.(5) / WARD 1974 EQ. (32)
    for i=1:3;
        obl_secular(i,:) = ( gj(i).*Ij(i)./(alpha_cosobj+gj(i)).*...
            (sin( (alpha_cosobj.*t+gj(i).*t+delta_j(i) - az_o - omega_o).*pi./180) - ...
            (sin( (delta_j(i) - az_o - omega_o).*pi./180) ) ) );
    end;
    obl = obl_J2000 - sum(obl_secular,1);
    obl = obl.*pi./180;

    %calculate the azimuth change in Titan's Spine Pole (WARD 1974)
    %%First WARD 1974 EQ 34
    for i=1:3;
        obl_secular2(i,:) = ( gj(i).*Ij(i)./(alpha_cosobj+gj(i)).*...
            (sin( (delta_j(i) - az_o - omega_o).*pi./180) ) );
    end;
    obl_o_star = obl_J2000 + sum(obl_secular2,1);
    %now define spin axis axzimuth (Ward 1974 EQ 36 and 37)
    obl_o = obl_J2000.*pi./180;
    for i=1:3;
        obl_secular3(i,:) = ( gj(i).*Ij(i)./(alpha_cosobj+gj(i))./sin(obl_o).*...
            (gj(i).*cos(obl_o) + alpha)./(gj(i)+alpha_cosobj).*...
            (cos( (alpha_cosobj.*t+gj(i).*t+delta_j(i) - az_o - omega_o).*pi./180) - ...
            (cos( (delta_j(i) - az_o - omega_o).*pi./180) ) ) );
    end;
    pole_lon = az_o - alpha.*t.*cos(obl_o_star.*pi./180) - ...
        capom(6,:).*180./pi + omega_o + sum(obl_secular3,1);
    pole_lon = mod(pole_lon.*pi./180,2.*pi);
    az = pole_lon;

    %calculate Ls of perihelion (Lsp) using the position of Titans spin pole
    Lsp = mod(w_p(6,:)-capom(6,:)-az-pi/2,2.*pi);

end;

%now organize outputs to match index
ecc_out = zeros(numel(t),length(index));
inc_out = ecc_out;
w_p_out = ecc_out;
capom_out = ecc_out;
if pos_saturn > 0;
    Lsp_out = ecc_out;
    obl_out = ecc_out;
    Lsp_out(:,pos_saturn) = Lsp;
    obl_out(:,pos_saturn) = obl;
end;
for i=1:length(index);
    ecc_out(:,i) = ecc(index(i),:);
    inc_out(:,i) = inc(index(i),:);
    w_p_out(:,i) = w_p(index(i),:);
    capom_out(:,i) = capom(index(i),:);
end;
%rename output files for convienience
ecc = ecc_out;
inc = inc_out;
w_p = w_p_out;
capom = capom_out;
if pos_saturn > 0;
    Lsp = Lsp_out;
    obl = obl_out;
end;
%% END CODE: orbelm_secular.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%