function [e,inc,w_p,capom,Lsp,obl,capm,omega,azo,f,a,jd,r_out,v_out] = orbelm_jpl(t,yr,idx);
%% HEADER: orbelm_jpl.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%Description: Calculates orbital elements from vector of julan dates 
%              (or years) using JPL ephemeris data files
%
%%Inputs:
%       t: vector of desired julian dates OR (IF YR FLAG SET):
%           vector of yrs since current epoch 
%          (valid range: -1e3 <= t <= 100)
%      yr: flag to indicate that t vector is in years, just has to exist or
%          be nonzero in length
%     idx: index corresponding to desired planet for ephmeris (default=6)
%         Possible Values:
%                         1: Mercury
%                         2: Venus
%                         3: Earth
%                         4: Asteroid Belt Barycenter
%                         5: Jupiter
%                         6: Saturn (default)
%                         7: Neptune
%                         8: Uranus
%                         9: Neptune
%                        10: Pluto
%                        11: Sun (will be 0) {NOT ALLOWED]
%                        12: Solar System BaryCenter
%                        13: Earth Moon BaryCenter
%
%%Outputs:
%       e: eccentricty
%       inc: incidence angle [radians]
%       w_p: longitude of perihelion [radians]
%       capom: longitude of ascending node [radians]
%       Lsp: Ls of perihelion (Saturn Only!) [radians]
%       obl: Obliquity  (Saturn Only!) [radians] (calculated from secular theory)
%       omega: argument of periapse [radians]
%       f: true anomaly [radians]
%       a: mean anomaly [radians]
%       jd: julian dates corresponding to input t vector
%       r_out: Cartesian position of saturn in JPL Horizons sun body centered
%              frame
%       v_out: Cartesian velocities of saturn in JPL Horizons sun body centered
%              frame
%
%%Calls: orbelm_xv2orbelm.m
%        init_eph.m (Matlab Ephemeris Package)
%        eph_global.m (Matlab Ephemeris Package)
%        pleph.m (Matlab Ephemeris Package)
%
%%Required Files: ASC_MERGE_405.dat
%                 ASC_MERGE_406.dat
%                 ASC_MERGE_414.dat
%
%%Source:
%      JPL Ephermis Data and Matlab Ephemeris Package
%
%%Author: Alexander G. Hayes (hayes@gps.caltech.edu)
%
%%Change Log: 
%            09/13/08: Draft Code Escalated to V1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% BEGIN CODE: orbelm_jpl.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%give zero length values to outputs to forestall warning messages
e = [];
inc = [];
w_p = [];
capom = [];
omega   = [];
f   = [];
capm   = [];
a   = [];
Lsp = [];
obl = [];
r_out = [];
v_out = [];
jd = [];

%idx defaults to saturn
if ~exist('idx','var');
    idx = 6;
elseif length(idx)==0;
    idx = 6;
end;
if length(idx) > 1;
    disp(['Only One IDX Value Allowed Per Run. Using position 1 (' ...
        num2str(idx(1)) ') ...']);
    idx = idx(1);
end;   
idx = idx( (idx > 1) & (idx < 14) );
if length(idx) == 0;
    disp(['IDX Values ' num2str(idx) ' Out Of Range (1 <= idx <= 13)']);
    return;
end;

%if we are inputing years, convert to julian date
if exist('yr','var');
    J2000 = 2451545; 
    %convert t to julian date from fractional year
    t = t.*365.25+J2000;
    jd = t;
end;

%define the julian date limit files
type = {'102','200','202','403','405','406','413','414'};
type_min = [2414992.5 2305424.5 2433264.5 2433264.5 2305424.5 1684496.5 2414992.5 2414992.5];
type_max = [2488080.5 2513392.5 2469808.5 2469808.5 2525008.5 2415056.5 2469808.5 2469808.5];
type_step= [32    32    32    32    32    64   32    32];
type_max(6) = 2414990;

%%initalize the correct data file
%first find the location of the minimum julian data
%ind1 = find( (type_min <= min(t)) & (type_max >= min(t)) );
ind1 = find( t < min(type_min) );
if length(ind1) == numel(t);
    disp(['Problem: Julian Dates of Input Time are All'  ...
        ' Too Old for Dataset (Min T = ' num2str(min(type_min)) ' ).' ...
        ' Use Alternative Method (Secular Theory or Swift Integration);' ...
        ' Returning...']);
    return;
end;
%define the output vectors
r_out = zeros(numel(t),3);
v_out = zeros(numel(t),3);
r_out(:) = nan;
v_out(:) = nan;
data_mask = zeros(size(t));
ind1 = find( (t >= min(type_min(6))) & (t <= max(type_max(6))) );


init_eph('ASC_MERGE_406.dat');
%initialize global variables
eph_global;
%get the ephemeris position in cartesian
if length(ind1) > 0;
    for i=ind1;      
        [r,v] = pleph(t(i),idx,11,0);
        if (sum(r) + sum(v)) ~= 0;
            r_out(i,:) = r;
            v_out(i,:) = v;
        end;
    end;
end;
%fill in the mask array
data_mask(ind1) = 1;
if length(ind1) ~= numel(t);
    ind2 = find( (t >= min(type_min(8))) & (t <= max(type_max(8))) );
    %initialize the ephemeris file
    init_eph(['ASC_MERGE_' type{8} '.dat']);
    %initialize global
    eph_global
    %get the ephemeris position in cartesian
    if length(ind2) > 0;
        for i=ind2
            [r,v] = pleph(t(i),idx,11,0);
            if (sum(r) + sum(v)) ~= 0;
                r_out(i,:) = r;
                v_out(i,:) = v;
            end;
        end;
    end;   
    %fill in the mask array
    data_mask(ind2) = 1;
end;
if length(find(data_mask)) ~= numel(t);
    ind3 = find( (t >= min(type_min(5))) & (t <= max(type_max(5))) );
    %initialize the ephemeris file
    init_eph(['ASC_MERGE_' type{5} '.dat']);
    %initialize global
    eph_global
    %get the ephemeris position in cartesian
    if length(ind3) > 0;
        for i=ind3;
            [r,v] = pleph(t(i),idx,11,0);
            if (sum(r) + sum(v)) ~= 0;
                r_out(i,:) = r;
                v_out(i,:) = v;
            end;
        end;
    end;   
    %fill in the mask array
    data_mask(ind3) = 1;
end;

if length(find(data_mask)) ~= numel(t);
    disp(['Warning: Not All Values of t Were In Integration Data Range']);
end;

%convert r,v into the orbital elements
ind = find( ~isnan(r_out(:,1)) & ~isnan(v_out(:,1)) );
if length(ind) > 0;
    %move r,v into a elciptic mean equinox coord system from the default
    %mean equinox earth equator system
    earth_obl = 23.45.*pi./180; %Earth Obliquity at J2000 (JPL Horizons)
    offset_x = -4.336703746110615e-9; %JPL Horizons
    offset_y = 1.68732754435296799262e-8; %JPL Horizons
    offset_z = -1.8896072528521624395e-4; %JPL Horizons
    %first negate the z-axis (orbit normal verse pole position)
    r_out(:,3) = -r_out(:,3);
    v_out(:,3) = -v_out(:,3);
    %rotate about the x-axis by earths obliquity
    Rx = [1 0               0; ...
          0 cos(earth_obl) -sin(earth_obl); ...
          0 sin(earth_obl)  cos(earth_obl); ];
    v_out = v_out';
    r_out = r_out';    
    v_out = [Rx * v_out]';
    r_out = [Rx * r_out]';
    r_out(:,1) = r_out(:,1)+offset_x;
    r_out(:,2) = r_out(:,2)+offset_y;
    r_out(:,3) = r_out(:,3)+offset_z;
    r_out(:,3) = -r_out(:,3);
    v_out(:,3) = -v_out(:,3);
    [a,e,inc,capom,omega,capm,f] = orbelm_xv2orbelm(r_out,v_out);
    w_p = omega+capom;
    
    %calculate the obliquity and pole precession values
    %if we are inputting julian date covert to years
    if 1;%~exist('yr','var');
        jd = t;
        J2000 = 2451545;
        %convert t to julian date from fractional year
        t = (t-J2000)./365.25;
    end;
    if idx == 6;
        %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 1973)
        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 (TITAN ONLY!)
        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 (TITAN ONLY)
    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'.*180./pi + omega_o + sum(obl_secular3,1);
    pole_lon = mod(pole_lon.*pi./180,2.*pi);
    azo = pole_lon;

    %calculate Ls of perihelion (Lsp) using the position of Titans spin pole
    Lsp = mod(omega'-azo-pi/2,2.*pi);                   
    end;
end;
%% END CODE: orbelm_jpl.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%