% Rosetta Earth and Mars Flyby Modeling Tool - main program
%
% Initial version: Magnus Billvik, Swedish Institute of Space Physics,
%Uppsala, June 2005
% Revised by: Niklas Edberg, Swedish Institute of Space Physics, Uppsala, Oct 2005
%
% Coordinate systems (all geocentric):
% ------------------------------------
% GEI (Geocentric Equatorial Inertial): X-aries,          Z-geo (rot) north
% GEA (Geocentric Ecliptic Aries):      X-aries,          Z-ecliptic north
% GSE (Geocentric Solar Ecliptic):      X-earth-sun line, Z-ecliptic north   <-- The most useful!
% GEO (Geographic):                     X-lat=long=0,     Z-geo (rot) north
% ROS (Rosetta Inertial System)         X-High-gain ant.  Z-the side with the ICA instrument
% MEI (Mars-centric Equat. Inertial):   X-aries           Z-geo (rot) north
% MEA (Mars-centric Ecliptic Aries):    X-aries           Z-ecliptic north
% MSE (Mars-centric Solar Ecliptic):    X-mars-sun line   Z-ecliptic north
% HEI (Heliocentric Equat. Inertial):   X-aries           Z-geo (rot) north
% HEA (Heliocentric Ecliptic Aries):    X-aries           Z-ecliptic north


if (~exist('dontread') | dontread==0)
    clear all;
    read=1;                   % read the trajectory and attitude data from files
end

% DEFINITIONS & CONSTANTS
rad_sun=695000;                              % radius of the sun (km)
rad_earth=6371.2;                            % radius of the Earth (km)
rad_mars=3397;                               % radius of Mars (km)
inclination=23.439291*pi/180;                % Earths axial inclination angle in rad
marsinclination=25.19*pi/180;                % Mars' axial inclination angle in rad
longascnode=49.57854*pi/180;                 % Longitude of ascending node of Mars
marsorbinc = 1.850*pi/180;                   % The inclination of Mars orbit to the ecliptic

disp('Rosetta Planetary Swingby Modeling tool - Magnus Billvik, May 2005');
disp('Select ROSETTA event to load trajectory and attitude data.');
disp('1 - Earth flyby 2005');
disp('2 - Earth flyby 2007');
disp('3 - Earth flyby 2009');
disp('4 - Mars flyby 2007');
disp('5 - whole mission');
event=input('?');

% To transform the Rosetta dates to julian days, Matlabs "datenum" is used:
timediff=datenum([-4713 1 1 12 0 0])+327;    % julian days between JD=0 and the year 0 ...
                                             % ..with a 327 day correction for I don't know what!
                                             % Add this value to the 'datenum' value to get Julian Days!
if read==1
    readdata;
end

% juilan date for rosetta position data
jdr=datenum(dater(1,:),dater(2,:),dater(3,:),dater(4,:),dater(5,:),dater(6,:))-timediff; 

% Here Rosetta's (GEI or MEI) coordinates are interpolated to match the time points
% for the coordinates of Sun and Moon (1 minute resolution). Also Mars' gea coords.

if event<=3
    rr_gei=zeros(3,length(jd)); % declare the variables to speed up
    vr_gei=zeros(3,length(jd)); 
    rr_gei(1,:) = interp1(jdr,RR(1,:),jd);      % interpolate the coordinates for every 1 minute
    rr_gei(2,:) = interp1(jdr,RR(2,:),jd);
    rr_gei(3,:) = interp1(jdr,RR(3,:),jd);
    vr_gei(1,:) = interp1(jdr,VR(1,:),jd);
    vr_gei(2,:) = interp1(jdr,VR(2,:),jd);
    vr_gei(3,:) = interp1(jdr,VR(3,:),jd);
    rs_gea=RS;
    if moon==1
        rm_gea=RMo;
        clear RMo;
    end
    
elseif event==4
    rr_mei=zeros(3,length(jd)); % declare the variables to speed up
    vr_mei=zeros(3,length(jd)); 
    rr_mei(1,:) = interp1(jdr,RR(1,:),jd,'spline'); %Interpolates with Splines just to get the closest approach more correct, doesn't really matter otherwise...
    rr_mei(2,:) = interp1(jdr,RR(2,:),jd,'spline');
    rr_mei(3,:) = interp1(jdr,RR(3,:),jd,'spline');
    vr_mei(1,:) = interp1(jdr,VR(1,:),jd,'spline');
    vr_mei(2,:) = interp1(jdr,VR(2,:),jd,'spline');
    vr_mei(3,:) = interp1(jdr,VR(3,:),jd,'spline');
    rs_gea=RS;
    rma_gea = RMa;
    if marsmoons==1
        rpho_gea=RPho;
        rdei_gea=RDei;
        clear RPho RDei;
    end
    
elseif event==5
    dist_earth=1.496e8;
    dist_mars=2.28e8;

    if att==0
        rr_hei=RR;
        vr_hei=VR;
        rc_hei=RC;
        jd=datenum(dater')-timediff;    % define the time from the trajectory data
        number=length(rr_hei);

    elseif att==1  % if attitude data is loaded: interpolate the trajectory

        jda_tmp=datenum(datea(1,:),datea(2,:),datea(3,:),datea(4,:),datea(5,:),datea(6,:))-timediff;  % juilan date for rosetta attitude data

        % delete attitude elements with the same timestamp (it causes error in the interpolation)
        Q_mod(:,1)=Q(:,1);
        jda(1)=jda_tmp(1);
        j=2;
        for i=2:length(jda_tmp)
            jda(j)=jda_tmp(i);
            Q_mod(:,j)=Q(:,i);
            if jda(j)==jda(j-1)
                j=j-1;
            end
            j=j+1;
        end
        clear jda_tmp;

        % interpolate the trajectory for the time steps defined by attitude data
        rr_hei(1,:)=interp1(jdr',RR(1,:),jda);
        rr_hei(2,:)=interp1(jdr',RR(2,:),jda);
        rr_hei(3,:)=interp1(jdr',RR(3,:),jda);
        vr_hei(1,:)=interp1(jdr',VR(1,:),jda);
        vr_hei(2,:)=interp1(jdr',VR(2,:),jda);
        vr_hei(3,:)=interp1(jdr',VR(3,:),jda);
        rc_hei(1,:)=interp1(jdc',RC(1,:),jda);
        rc_hei(2,:)=interp1(jdc',RC(2,:),jda);
        rc_hei(3,:)=interp1(jdc',RC(3,:),jda);
        vc_hei(1,:)=interp1(jdc',VC(1,:),jda);
        vc_hei(2,:)=interp1(jdc',VC(2,:),jda);
        vc_hei(3,:)=interp1(jdc',VC(3,:),jda);

        number=length(rr_hei);
      
        % Interpolate the quaternions for every minute
        q1 = Q_mod(1,:);
        q2 = Q_mod(2,:);
        q3 = Q_mod(3,:);
        q4 = Q_mod(4,:);
        clear Q_mod;
%        keyboard;
%        q4(find(q4>0.6))=q4(find(q4>0.6))-0.9;

        % Define rosettas attitude vectors (from the quaternions) in GEI frame of referece for all times (3d matrix),
        % the Rosetta inertial x,y,z axis being the rows of Asc
        % Formula for converting quaternions to attitude vectors are taken from the ESA document file: 
        % RO-ESC-IF-5003_2_0_DDID_Appendix_H_Data_Delivery_FD_Products_2003Oct23.pdf
        for i=1:number
            Asc_gei(1,:,i) = [q1(i)^2-q2(i)^2-q3(i)^2+q4(i)^2  2*(q1(i)*q2(i)+q3(i)*q4(i))      2*(q1(i)*q3(i)-q2(i)*q4(i))];
            Asc_gei(2,:,i) = [2*(q1(i)*q2(i)-q3(i)*q4(i))     -q1(i)^2+q2(i)^2-q3(i)^2+q4(i)^2  2*(q2(i)*q3(i)+q1(i)*q4(i))];
            Asc_gei(3,:,i) = [2*(q1(i)*q3(i)+q2(i)*q4(i))      2*(q2(i)*q3(i)-q1(i)*q4(i))     -q1(i)^2-q2(i)^2+q3(i)^2+q4(i)^2];
        end

        jd=jda;                 % define the time from the attitude data
%        clear q1 q2 q3 q4;
        clear jda jdr jdc;  % memory cleanup
    end

    istart=1;
    iend=number;
    
end
   
    
%clear RR VR RS RMa;
%clear jdmo jdr;

% ATIITUDE data - interpolate the quarternions and calculate the att. vectors (for flybys)
if att==1 & event<5

    jda_tmp=datenum(datea(1,:),datea(2,:),datea(3,:),datea(4,:),datea(5,:),datea(6,:))-timediff;  % juilan date for rosetta attitude data
    % delete attitude elements with the same timestamp (it causes error in the interpolation)
    Q_mod(:,1)=Q(:,1);
    jda(1)=jda_tmp(1);
    j=2;
    for i=2:length(jda_tmp)
        jda(j)=jda_tmp(i);
        Q_mod(:,j)=Q(:,i);
        if jda(j)==jda(j-1)
            j=j-1;
        end
        j=j+1;
    end
    clear jda_tmp;
    
    % Interpolate the quaternions for every minute
    q1 = interp1(jda,Q_mod(1,:),jd);
    q2 = interp1(jda,Q_mod(2,:),jd);
    q3 = interp1(jda,Q_mod(3,:),jd);
    q4 = interp1(jda,Q_mod(4,:),jd);
    clear Q_mod;

    % Define rosettas attitude vectors (from the quaternions) in GEI frame of referece for all times (3d matrix),
    % the Rosetta inertial x,y,z axis being the rows of Asc
    % Formula for converting quaternions to attitude vectors are taken from the ESA document file: 
    % RO-ESC-IF-5003_2_0_DDID_Appendix_H_Data_Delivery_FD_Products_2003Oct23.pdf
    for i=1:number
        Asc_gei(1,:,i) = [q1(i)^2-q2(i)^2-q3(i)^2+q4(i)^2  2*(q1(i)*q2(i)+q3(i)*q4(i))      2*(q1(i)*q3(i)-q2(i)*q4(i))];
        Asc_gei(2,:,i) = [2*(q1(i)*q2(i)-q3(i)*q4(i))     -q1(i)^2+q2(i)^2-q3(i)^2+q4(i)^2  2*(q2(i)*q3(i)+q1(i)*q4(i))];
        Asc_gei(3,:,i) = [2*(q1(i)*q3(i)+q2(i)*q4(i))      2*(q2(i)*q3(i)-q1(i)*q4(i))     -q1(i)^2-q2(i)^2+q3(i)^2+q4(i)^2];
    end

    clear q1 q2 q3 q4; % memory cleanup
    clear jda;

end

coordtransform;


% define index variables for the interesting range of the flyby
% find the first and last index with traj data that is not NaN
if event<4 % for the earth flybys
    i=1;
    while isnan(rr_gea(1,i))
        i=i+1;
    end
    istart=i;
    i=number;
    while isnan(rr_gea(1,i))
        i=i-1;
    end
    iend=i;
    mindist=min(distance);
    iclosest=find(mindist==distance);
elseif event==4 % for the mars flyby
    i=1;
    while isnan(rr_mea(1,i))
        i=i+1;
    end
    istart=i;
    i=number;
    while isnan(rr_mea(1,i))
        i=i-1;
    end
    iend=i;
    mindist=min(distance);
    iclosest=find(mindist==distance);
    mindistpho=min(disttopho);
    iclosepho=find(mindistpho==disttopho);
    mindistdei=min(disttodei);
    iclosedei=find(mindistdei==disttodei);
end

disp('Calculations done!');


vis;
