%Visulatization program for the Rosetta trajectory
%
% Initial version: Magnus Billvik, Swedish Institute of Space Physics, Uppsala, June 2005
% Revised by: Niklas Edberg, Swedish Institute of Space Physics, Uppsala, Oct 2005

disp('Select visualization type/option:');
disp('a - distance plot');
disp('b - 3D GSE (or MSO) plot');
disp('c - 2D GSE (or MSO) plots');
disp('d - 3D GSE plot with s/c attitude');
disp('e - 2D GSE plots with s/c attitude');
disp('f - GEO plot with world map (only earth flybys)');
disp('g - LAP: sun and velocity vector');
disp('h - LAP: make image file');
disp('i - LAP: Calculate "sunlit/eclipse and wake/non wake" from image file');
disp('j - LAP: Plot "sunlit/eclipse and wake/non wake" vs. time');
disp('k - Whole mission plot');
disp('l - set current time span');
disp('m - dump trajectory data file');
disp('n - dump eclipse/wake data file');
disp('o - LAP: plot measured photo-current data in lap plot');
disp('p - MEO plot with Mars map (only Mars flyby)');
disp('q - 2D MSO with Mars mangetosphere (only Mars flyby)');
disp('r - 2D MRX plots with Mars rotation axis (only Mars flyby)');

item=input('?','s');

if isempty(get(0, 'CurrentFigure'))
    fign=1;
else
    fign=get(0, 'CurrentFigure')+1;
end

windowsize=[0 0 1275 945]; % for maximizing the window size (probably different for different computers...)
scol=0; %Sets the color on the sphere that is beeing drawn, 0 for Earth (default) and 1 for Mars

if event<4
    x=rr_gse(1,:)/radius; %Divide by radius to get the scale in planet radii
    y=rr_gse(2,:)/radius;
    z=rr_gse(3,:)/radius;
    mx=rm_gse(1,:)/radius;
    my=rm_gse(2,:)/radius;
    mz=rm_gse(3,:)/radius;
    index=[istart:iend];
    indexx=[istart:60:iend];
    coordsyst=' GSE';

    
elseif event==4
    
    x=rr_mso(1,:)/radius;
    y=rr_mso(2,:)/radius;
    z=rr_mso(3,:)/radius;
    
    xpho=rpho_mso(1,:)/radius;
    ypho=rpho_mso(2,:)/radius;
    zpho=rpho_mso(3,:)/radius;
    
    xdei=rdei_mso(1,:)/radius;
    ydei=rdei_mso(2,:)/radius;
    zdei=rdei_mso(3,:)/radius;
    
    index=[istart:iend];
    indexx=[istart:60:iend];
    coordsyst=' MSO';
  
    marsax_mrx(1,:)     =   zeros(1,length(rr_mso)); %These lines sets up the martian rotation axis. It should have components (0,0,1) in MRX.
    marsax_mrx(2,:)     =   zeros(1,length(rr_mso));
    marsax_mrx(3,:)     =   ones(1,length(rr_mso)).*3;
    marsax_mei          =   mrx2mei(marsax_mrx);
    marsax_mea          =   mei2mea(marsax_mei,inclination); %It is then transformed to mea -> moc -> mso-coordinates
    marsax_moc          =   mea2moc(marsax_mea);
    marsax_mso          =   moc2mso(marsax_moc,rs_moc);
    
    marsver_mrx(1,:)    =   ones(1,length(rr_mso)).*3; %These lines sets up the direction to the mars vernal equinox. It should have components (1,0,0) in MRX.
    marsver_mrx(2,:)    =   zeros(1,length(rr_mso));
    marsver_mrx(3,:)    =   zeros(1,length(rr_mso));
    marsver_mei         =   mrx2mei(marsver_mrx);
    marsver_mea         =   mei2mea(marsver_mei,inclination); %It is then transformed to mea -> moc -> mso-coordinates
    marsver_moc         =   mea2moc(marsver_mea);
    marsver_mso         =   moc2mso(marsver_moc,rs_moc);
    
    marszxx_mrx(1,:)    =   zeros(1,length(rr_mso)); %These lines sets up the direction to the y-axis. It should have components (0,1,0) in MRX.
    marszxx_mrx(2,:)    =   ones(1,length(rr_mso)).*3;
    marszxx_mrx(3,:)    =   zeros(1,length(rr_mso));
    marszxx_mei         =   mrx2mei(marszxx_mrx);
    marszxx_mea         =   mei2mea(marszxx_mei,inclination); %It is then transformed to mea -> moc -> mso-coordinates
    marszxx_moc         =   mea2moc(marszxx_mea);
    marszxx_mso         =   moc2mso(marszxx_moc,rs_moc);   
    
    
elseif event==5
    x=rr_hea(1,:)/radius;
    y=rr_hea(2,:)/radius;
    z=rr_hea(3,:)/radius;
    index=[istart:iend];
%    indexx=[istart:60:iend];
    coordsyst=' HEA';
end



switch item

% Altitude plot
case 'a'
    figure(fign)
    set(gcf,'Position', windowsize);   
    sta=25;
    sto=25;
    plot(jd((iclosest-sta:iclosest+sta))+timediff, (distance((iclosest-sta:iclosest+sta))-radius));
    hold on;
    plot([jd(iclosest-sta) jd(iclosest+sto)]+timediff, [1 1]*(distance(iclosest)-radius), 'k--');
    ax=axis;
    axis([ax(1) ax(2) 0 ax(4)]);
    ylabel('altitude from planetary surface (km)');
    titlestr=strcat(tit, ' - Altitude plot (',datestr(jd(iclosest-sta)+timediff,20), '@',...
        datestr(jd(iclosest-sta)+timediff,15), '-',datestr(jd(iclosest+sto)+timediff,20),'@',...
        datestr(jd(iclosest+sto)+timediff,15),')');
    title(titlestr);
    datetick('x', 15);


% Plot 3D path
case 'b'
    figure(fign)
    set(gcf,'Position', windowsize);
    
    %Nytt
    %rr_mse=gea2gse(rr_mea,rs_mea)./radius;
    %plot3(rr_mse(1,:),rr_mse(2,:),rr_mse(3,:))
    %hold on
    %plot3(rr_mse(1,885),rr_mse(2,885),rr_mse(3,885),'*')
    %plot3(rr_mso(1,iend-895),rr_mso(2,iend-895),rr_mso(3,iend-895),'*')
    %axis([-10 10 -30 30 -2 2])
    %Slut nytt
    
    plot3(x(index),y(index),z(index));
    hold on;
    plot3(x(istart),y(istart),z(istart),'*');
    if moon==1
        plot3(mx(index),my(index),mz(index),'y');    % moon
        plot3(mx(istart), my(istart), mz(istart), 'yo'); % plot a ring at the start position
    end
    if marsmoons == 1
        scol=1;
        plot3(xpho(index),ypho(index),zpho(index),'r');
        plot3(xdei(index),ydei(index),zdei(index),'g');
        plot3(xpho(istart), ypho(istart), zpho(istart), 'r*');
        plot3(xdei(istart), ydei(istart), zdei(istart), 'g*');
       
        plot3(xpho(iclosepho), ypho(iclosepho), zpho(iclosepho), 'rs'); %Plot the closest approach to the moon
        plot3(x(iclosepho), y(iclosepho), z(iclosepho), 'rs');
       
        plot3(xdei(iclosedei), ydei(iclosedei), zdei(iclosedei), 'gs'); %Plot the closest approach to the moon
        plot3(x(iclosedei), y(iclosedei), z(iclosedei), 'gs');
        %plot the martian rotation axis and the direction to mars vernal equinox and the y-axis 
        quiver3([zeros(1,length(rr_mso))],[zeros(1,length(rr_mso))],[zeros(1,length(rr_mso))],marsax_mso(1,:),marsax_mso(2,:),marsax_mso(3,:)); 
        quiver3([zeros(1,length(rr_mso))],[zeros(1,length(rr_mso))],[zeros(1,length(rr_mso))],marsver_mso(1,:),marsver_mso(2,:),marsver_mso(3,:));
        quiver3([zeros(1,length(rr_mso))],[zeros(1,length(rr_mso))],[zeros(1,length(rr_mso))],marszxx_mso(1,:),marszxx_mso(2,:),marszxx_mso(3,:));
        
        %Plot the start position for axis([-30 30 -10 10 -2 2])
        plot3(rr_mso(1,885)/radius,rr_mso(2,885)/radius,rr_mso(3,885)/radius,'*')
    end
    sphere(1,scol);
    axis('equal');
    grid on;
    %axis([-30 30 -10 10 -2 2])
    splot3(x(index),y(index),z(index),axis)
    if moon==1
        splot3(mx(index),my(index),mz(index),axis);
    end
    ssphere(1,axis);

    xlabel('x (planetary radii) - towards sun');
    ylabel('y (planetary radii)');
    zlabel('z (planetary radii)');
    %titlestr=strcat(tit, ' - 3D ', 'MSE',' plot (',datestr(jd(885)+timediff,20),...
     %   '@', datestr(jd(885)+timediff,15), '-',datestr(jd(iend-895)+timediff,20),'@',...
      %  datestr(jd(iend-895)+timediff,15),')'); 
    
    titlestr=strcat(tit, ' - 3D ', coordsyst,' plot (',datestr(jd(istart)+timediff,20),...
        '@', datestr(jd(istart)+timediff,15), '-',datestr(jd(iend)+timediff,20),'@',...
        datestr(jd(iend)+timediff,15),')');
    title(titlestr);
    if moon==1
        legend('Rosetta path', 'start', 'Moon path', 'start');
    elseif marsmoons ==1
         legend('Rosetta path', 'Rosetta start', 'Phobos path','Deimos path', 'Phobos start', 'Deimos start','Point closest to Phobos','Point closest to Phobos','Point closest to Deimos','Point closest to Deimos','Mars spin axis','To Mars vernal equinox');
    end
    %axis([-10 10 -10 10 -2 2]);

% Plot 2D projections
case 'c'
    windowsize=[0 0 1000 1000];
    figure(fign)
    set(gcf,'Position', windowsize); 

    subplot(2,2,1)
    plot(x(index),y(index));            % rosetta
    hold on;
    plot(x(istart), y(istart), '*');    % plot a star at the start position
    if moon==1
        plot(mx(index),my(index), 'y');     % moon
        plot(mx(istart), my(istart), 'y*'); % plot a ring at the start position
    end
    if marsmoons==1
        plot(xpho(index),ypho(index),'r');
        plot(xdei(index),ydei(index),'g');
    end
    circle(0,0,1,'k');                  % earth/mars
    axis('equal');
    grid on;
    xlabel('x (planet radii)');
    ylabel('y (planet radii)');
    title('X/Y');
    axis([-12 12 -12 12]);
    
    subplot(2,2,2)
    plot(x(index),z(index));            % rosetta
    hold on;
    plot(x(istart), z(istart), '*');    % plot a ring at the start position
    if moon==1
        plot(mx(index),mz(index),'y');      % moon
        plot(mx(istart), mz(istart), 'y*'); % plot a ring at the start position
    end 
    if marsmoons==1
        plot(xpho(index),zpho(index),'r');
        plot(xdei(index),zdei(index),'g');
    end
    circle(0,0,1,'k');                  % earth
    axis('equal');
    grid on;
    xlabel('x (planet radii)');
    ylabel('z (planet radii)');
    title('X/Z');
    axis([-12 12 -12 12]);
    
    subplot(2,2,3)
    plot(y(index),z(index));            % rosetta
    hold on;
    plot(y(istart), z(istart), '*');    % plot a ring at the start position
    if moon==1
        plot(my(index),mz(index),'y');      % moon
        plot(my(istart), mz(istart), 'y*'); % plot a ring at the start position
    end
    if marsmoons==1
        plot(ypho(index),zpho(index),'r');
        plot(ydei(index),zdei(index),'g');
    end
    circle(0,0,1,'k');                  % earth
    axis('equal');
    grid on;
    xlabel('y (planet radii)');
    ylabel('z (planet radii)');
    title('Y/Z');

%    subplot(2,2,4)
%    plot(0,0);
%    plot(0,0,'*');
%    legend('path', 'start');
    
    title(strcat(tit, ' 2D', coordsyst,' plots (',datestr(jd(istart)+timediff,20), '@',...
        datestr(jd(istart)+timediff,15), '-',datestr(jd(iend)+timediff,20),'@',...
        datestr(jd(iend)+timediff,15),') [printed ', datestr(now), ']'));
    
    axis([-12 12 -12 12]);
    
% Plot 3D projections and attitude vectors
case 'd'
    figure(fign)
    set(gcf,'Position', windowsize); 

    % following is the coordinates for a square representing rosetta (only to make the viewing of the xyz axis easier)
    ros3dmodela_sc=0.5*[ 1 0 1 ; 1 1 1 ; 1 1 0 ; 1 0 0 ; 0 0 1 ; 0 1 1 ; 0 1 0 ; 0 0 0 ; 1 0 1 ; 1 1 1 ; 1 1 0 ; 1 0 0 ]';
    ros3dmodelb_sc=0.5*[ 1 1 1 ; 1 1 0 ; 1 0 0 ; 1 0 1 ; 0 1 1 ; 0 1 0 ; 0 0 0 ; 0 0 1 ; 0 0 1 ; 0 1 1 ; 0 1 0 ; 0 0 0 ]';
    x1=Asc_gse(1,1,:); x2=Asc_gse(1,2,:); x3=Asc_gse(1,3,:); % xyz gse components of the rosetta's x-axis
    y1=Asc_gse(2,1,:); y2=Asc_gse(2,2,:); y3=Asc_gse(2,3,:); % xyz gse components of the rosetta's y-axis
    z1=Asc_gse(3,1,:); z2=Asc_gse(3,2,:); z3=Asc_gse(3,3,:); % xyz gse components of the rosetta's z-axis

    plot3(x(index),y(index),z(index));       % rosetta
    hold on;
    if moon==1
        plot3(mx(index),my(index),mz(index),'y');  % moon
        plot3(mx(istart), my(istart), mz(istart), 'yo');            % plot a ring at the start position
        splot3(mx(index),my(index),mz(index),axis);
    end
    sphere(1,0);            % earth
    splot3(x(index),y(index),z(index),axis);                    % gray 2d-projections on the sides
    ssphere(1,axis);                                            % earth
    plot3(x(istart), y(istart), z(istart), 'o');                % a ring at the start position

    % plot the 3d cube
    for i=1:length(indexx)
        a=inv(Asc_gse(:,:,indexx(i)))*ros3dmodela_sc; % transform to GSE ref frame
        b=inv(Asc_gse(:,:,indexx(i)))*ros3dmodelb_sc;
        for l=1:12
            a(1,l)=a(1,l)+x(indexx(i)); % translation to earth centered origo
            a(2,l)=a(2,l)+y(indexx(i));
            a(3,l)=a(3,l)+z(indexx(i));
            b(1,l)=b(1,l)+x(indexx(i));
            b(2,l)=b(2,l)+y(indexx(i));
            b(3,l)=b(3,l)+z(indexx(i));
            plot3([a(1,l) b(1,l)], [a(2,l) b(2,l)], [a(3,l) b(3,l)], 'color', [0.7 0.7 0.7]);
        end
    end
    % plot rosettas x,y and z axis in rgb colours
    point3d(x(indexx), y(indexx), z(indexx), x1(indexx), x2(indexx), x3(indexx), 3, 'r');
    point3d(x(indexx), y(indexx), z(indexx), y1(indexx), y2(indexx), y3(indexx), 3, 'g');
    point3d(x(indexx), y(indexx), z(indexx), z1(indexx), z2(indexx), z3(indexx), 3, 'b');

    axis('equal');
    grid on;
    xlabel('x (planet radii)');
    ylabel('y (planet radii)');
    zlabel('z (planet radii)');
    titlestr=strcat(tit, ' - 3D ', coordsyst,' plot (',datestr(jd(istart)+timediff,20),...
        '@', datestr(jd(istart)+timediff,15), '-',datestr(jd(iend)+timediff,20),'@',...
        datestr(jd(iend)+timediff,15),')');
    title(titlestr);
    legend('rosetta trajectory', 'moon trajectory');
    
    clear ros3dmodela_sc ros3dmodelb_sc x1 x2 x3 y1 y2 y3 z1 z2 z3 a b;
    
    
% Plot 2D path and attitude vectors
case 'e'

    x1=Asc_gse(1,1,:); x2=Asc_gse(1,2,:); x3=Asc_gse(1,3,:); % xyz gse components of the rosetta's x-axis
    y1=Asc_gse(2,1,:); y2=Asc_gse(2,2,:); y3=Asc_gse(2,3,:); % xyz gse components of the rosetta's y-axis
    z1=Asc_gse(3,1,:); z2=Asc_gse(3,2,:); z3=Asc_gse(3,3,:); % xyz gse components of the rosetta's z-axis

    figure(fign)
    set(gcf,'Position', windowsize); 

    subplot(2,2,1)
    plot(x(index),y(index));        % rosetta
    hold on;
    plot(x(istart), y(istart), 'o');    % plot a ring at the start position
    if moon==1
        plot(mx(index),my(index),'y');      % moon
        plot(mx(istart), my(istart), 'yo'); % plot a ring at the start position
    end
    circle(0,0,1,'k');                  % earth
    % plot rosettas x,y and z axis in rgb colours
    point(x(indexx), y(indexx), x1(indexx), x2(indexx), 4, 'r');
    point(x(indexx), y(indexx), y1(indexx), y2(indexx), 4, 'g');
    point(x(indexx), y(indexx), z1(indexx), z2(indexx), 4, 'b');
    axis('equal');
    grid on;
    xlabel('x (planet radii)');
    ylabel('y (planet radii)');
    title('X/Y');
  
    subplot(2,2,2)
    plot(x(index),z(index));        % rosetta
    hold on;
    plot(x(istart), z(istart), 'o');    % plot a ring at the start position
    if moon==1
        plot(mx(index),mz(index),'y');      % moon
        plot(mx(istart), mz(istart), 'yo'); % plot a ring at the start position
    end
    circle(0,0,1,'k');                  % earth
    % plot rosettas x,y and z axis in rgb colours
    point(x(indexx), z(indexx), x1(indexx), x3(indexx), 4, 'r');
    point(x(indexx), z(indexx), y1(indexx), y3(indexx), 4, 'g');
    point(x(indexx), z(indexx), z1(indexx), z3(indexx), 4, 'b');
    axis('equal');
    grid on;
    xlabel('x (planet radii)');
    ylabel('z (planet radii)');
    title('X/Z');


    subplot(2,2,3)
    plot(y(index),z(index));        % rosetta
    hold on;
    plot(y(istart), z(istart), 'o');    % plot a ring at the start position
    if moon==1
        plot(my(index),mz(index),'y');      % moon
        plot(my(istart), mz(istart), 'yo'); % plot a ring at the start position
    end
    circle(0,0,1,'k');                  % earth
    % plot rosettas x,y and z axis in rgb colours
    point(y(indexx), z(indexx), x2(indexx), x3(indexx), 4, 'r');
    point(y(indexx), z(indexx), y2(indexx), y3(indexx), 4, 'g');
    point(y(indexx), z(indexx), z2(indexx), z3(indexx), 4, 'b');
    axis('equal');
    grid on;
    xlabel('y (planet radii)');
    ylabel('z (planet radii)');
    title('Y/Z');

    suptitle(strcat(tit, ' 2D', coordsyst,' plots with attitude (',datestr(jd(istart)+timediff,20),...
        '@', datestr(jd(istart)+timediff,15), '-',datestr(jd(iend)+timediff,20),'@',...
        datestr(jd(iend)+timediff,15),')'));
    clear x1 x2 x3 y1 y2 y3 z1 z2 z3;

    
% Plot Rosetta's path on a World map    
case 'f'
    figure(fign)
    set(gcf,'Position', windowsize); 
    [world, col]=imread('world_blueblack.png','png');
    newplot;
    set(gca,'YDir','normal');
    %image([-10 350], [85 -60], world, 'YData', [85 -60]), colormap(col);
    image([0 360], [90 -90], world, 'YData', [90 -90]), colormap(col);
    set(gca,'YDir','normal');
    hold on;
    plot(longitude(istart:iend),latitude(istart:iend),'r.');    % the s/c path
    plot(longitude(istart),latitude(istart),'b.');              % the start of the data

    temp{1,1}=datestr(jd(iclosest)+timediff);
    temp{2,1}=strcat('lat=', num2str(latitude(iclosest),'%3.1f'), '\circ');
    temp{3,1}=strcat('long=', num2str(longitude(iclosest),'%3.1f'), '\circ');
    temp{4,1}=strcat('alt=', num2str(round(distance(iclosest)-radius),'%4.0f'), ' km');
    h=text(longitude(iclosest)-0.4, latitude(iclosest)+0.4, temp, 'FontWeight', 'bold');
    set(h, 'Color', 'w');
    h=text(longitude(iclosest), latitude(iclosest), temp, 'FontWeight', 'bold');
    set(h, 'Color', 'b');

    temp{1,1}=datestr(jd(istart)+timediff);
    temp{2,1}=strcat('lat=', num2str(latitude(istart),'%3.1f'), '\circ');
    temp{3,1}=strcat('long=', num2str(longitude(istart),'%3.1f'), '\circ');
    temp{4,1}=strcat('alt=', num2str(round(distance(istart)-radius),'%4.0f'), ' km');
    h=text(longitude(istart)-0.4, latitude(istart)+0.4,temp, 'FontWeight', 'bold');
    set(h, 'Color', 'w');
    h=text(longitude(istart), latitude(istart),temp, 'FontWeight', 'bold');
    set(h, 'Color', 'b');

    plot(longitude(istart), latitude(istart),'bo');

    xlabel('longitude');
    ylabel('latitude');
    titlestr=strcat(tit,' - GEO long/lat plot (',datestr(jd(istart)+timediff,20), ...
        '@', datestr(jd(istart)+timediff,15), '-',datestr(jd(iend)+timediff,20),'@',...
        datestr(jd(iend)+timediff,15),')');
    title(titlestr);
    axis([0 360 -90 90]);

    clear temp h;

% Plot sun vector or velocity vector as seen from LAP
case 'g'
    proben=input('LAP probe number: ');
    temp=input('Sun vector (1), velocity vector (2) or none (3)');
    relradius=input('High-gain antenna relative radius (1=100%): ');
    dontplot=0;
    
    switch temp
    case 1
        phi=sphi;
        theta=stheta;
        titlestr=strcat(tit,' - Movement of the Sun as seen from LAP probe #', num2str(proben),...
            ' (',datestr(jd(istart)+timediff,20), '@', datestr(jd(istart)+timediff,15), '-',...
            datestr(jd(iend)+timediff,20),'@',datestr(jd(iend)+timediff,15),')');
    case 2
        phi=vphi;
        theta=vtheta;
        titlestr=strcat(tit,' - Movement of Rosettas velocity vector as seen from LAP probe #', ...
            num2str(proben), ' (',datestr(jd(istart)+timediff,20), '@', datestr(jd(istart)+timediff,15),...
            '-',datestr(jd(iend)+timediff,20),'@',datestr(jd(iend)+timediff,15),')');
    case 3
        dontplot=1;
        titlestr=strcat('The Rosetta spacecraft as seen from LAP probe #', num2str(proben));
    end        

    % to make the plots look fine, take away the point before the "jump" (but save the orig. values)
    figure(fign)
    set(gcf,'Position', windowsize);

    if dontplot==0
        phi_save=phi;
        theta_save=theta;
        itmp=find(abs(diff(phi))>300);
        phi(itmp)=NaN;
        theta(itmp)=NaN;

        plot(phi(index), theta(index), 'r');
        hold on;
        plot(phi(istart), theta(istart), 'r*');
    end
    drawrosetta(proben, relradius);
    if dontplot==0
        plot(phi(istart), theta(istart), 'r*');
        plot(phi(index), theta(index), 'r');
        grid on;
   end
    axis([-180 180 -90 90]);
    set(gca,'XDir','reverse');
    xlabel('azimuth in zx plane (x -> -z)');
    ylabel('elevation above the zx plane');
    title(titlestr);
    if dontplot==0
        legend('path', 'start');
    end

    clear phi theta phi_save theta_save;


% Save LAP image file
case 'h'
    proben=input('LAP probe number: ');

    figure(fign)
    set(gcf,'Position', windowsize);
    set(gca,'XDir','reverse');          % reverses the phi-axis
    set(gcf,'Position',[0 0 1275 945]); % maximizes the window size
    axis([-180 180 -90 90]);
    
    drawrosetta(proben,0);              % draw the spacecraft without the antenna
    set(gca,'XDir','reverse');          % reverses the phi-axis
    axis([-180 180 -90 90]);
    filename=strcat('image_p', num2str(proben), '_1');
    print(gcf, '-dpng', '-r400', filename)

    hold off;
    drawrosetta(proben,1);              % draw the spacecraft with the antenna
    set(gca,'XDir','reverse');          % reverses the phi-axis
    axis([-180 180 -90 90]);
    filename=strcat('image_p', num2str(proben), '_2');
    print(gcf, '-dpng', '-r400', filename)

    hold off;
    drawrosetta(proben,1.1);            % draw the spacecraft with an extra big antenna (110% radius)
    set(gca,'XDir','reverse');          % reverses the phi-axis
    axis([-180 180 -90 90]);
    filename=strcat('image_p', num2str(proben), '_3');
    print(gcf, '-dpng', '-r400', filename)


    disp('Images saved in file "image_pX_X.png". Manually convert the files to INDEXED GIF before calculating wake/eclipse data.');

    clear filename;


% Calculate wake- and eclipse data from image file
case 'i'
    proben=input('LAP probe number: ');
    clear scindex vcindex;
    scindex=zeros(3,length(sphi));
    vcindex=zeros(3,length(sphi));
    for n=1:3               % loop through the three image files
        filename=strcat('image_p',num2str(proben),'_', num2str(n));
        disp(strcat('Reading image file "',filename , '.gif"...'));
        [I, colmap] = imread(filename,'gif');
        disp('Cropping image...');
        img = imcrop(I,[418 181 2477 1954]);
        colbg=img(1,1);     % index number for the (white) background color
        disp('Calculating...');
        disp(strcat('Size of image matrix: ',num2str(size(img))));

        % this loop checks whether the coordinate is "behind" the s/c or visible from the probe
        for i=1:length(sphi)
            if isnan(sphi(i))       % if data values are missing..
                scindex(n,i)=NaN;
                vcindex(n,i)=NaN;
            else
                [sx,sy]=imagexy(sphi(i),stheta(i));     % calculate the "pixel" corresponding to a 
                [vx,vy]=imagexy(vphi(i),vtheta(i));     % certain phi-theta pair
                scindex(n,i)=img(sy,sx);                % (color) value for the image coordiante/pixel
                vcindex(n,i)=img(vy,vx);
                if scindex(n,i)==colbg                  % if visible (sunlit) set the value to 0
                    scindex(n,i)=0;
                else                                    % if in the shadow of the s/c set the value to 1
                    scindex(n,i)=1;
                end
                if vcindex(n,i)==colbg                  % if not in wake
                    vcindex(n,i)=0;
                else                                    % if in wake
                    vcindex(n,i)=1;
                end
            end
        end
    end
    


    switch proben
    % add together the values calculated from the three images 
    % to get non-binary variables for the eclipse and wake data
    case 1
        wake1=(vcindex(1,:)+vcindex(2,:)+vcindex(3,:))/3;
        eclipse1=(scindex(1,:)+scindex(2,:)+scindex(3,:))/3;
        disp('Wake-data for the current time span saved in variable "wake1"');
        disp('Eclipse-data for the current time span saved in variable "eclipse1"');
    case 2
        wake2=(vcindex(1,:)+vcindex(2,:)+vcindex(3,:))/3;
        eclipse2=(scindex(1,:)+scindex(2,:)+scindex(3,:))/3;
        disp('Wake-data for the current time span saved in variable "wake2"');
        disp('Sunlit-data for the current time span saved in variable "Eclipse2"');
    end

    clear vcindex scindex sx sy vx vy I colmap bgcol filename;


% Plot wake- and sunlit data vs. time
case 'j'
    proben=input('LAP probe number: ');
    switch proben
    case 1
        wake=wake1;
        eclipse=eclipse1;
    case 2
        wake=wake2;
        eclipse=eclipse2;
    end

    figure(fign)
    set(gcf,'Position', windowsize);

    subplot(2,1,1);
    plot(jd(index)+timediff,eclipse(index));
    grid on;
    axis([jd(istart)+timediff jd(iend)+timediff -0.2 1.2]);
    datetick('x',19);
    xlabel('time');
    ylabel('eclipse index');
    ax=axis;
    title('Eclipse information');

    subplot(2,1,2);
    plot(jd(index)+timediff,wake(index));
    grid on;
    axis([jd(istart)+timediff jd(iend)+timediff -0.2 1.2]);
    datetick('x',19);
    xlabel('time');
    ylabel('wake index');
    ax=axis;
    title('Wake information');

%    suptitle(strcat(tit, ' LAP probe #',num2str(proben), ' [printed ', datestr(now), ']'));
    suptitle(strcat(tit, ' LAP probe ',num2str(proben), ' (',datestr(jd(istart)+timediff,20),...
        '@', datestr(jd(istart)+timediff,15),'-',datestr(jd(iend)+timediff,20),'@',...
        datestr(jd(iend)+timediff,15),')'));
    
    clear wake Eclipse ax titlestr;
   
case 'k'
    % plot the whole mission
    figure(fign)
    set(gcf,'Position', windowsize); 
    AU=149597870.691;                   % define the astronomical unit
    plot3(rr_hea(1,:), rr_hea(2,:), rr_hea(3,:));
    hold on;
    plot3(rc_hea(1,:), rc_hea(2,:), rc_hea(3,:), 'k--');
    plot3(0,0,0,'y*');
    % plot mars and earth orbit
    theta=(0:1:360);
    [xx,yy] = pol2cart(theta*pi/180,1);
    plot3(dist_earth*xx, dist_earth*yy, zeros(1,length(xx)),'g--');
    plot3(dist_mars*xx, dist_mars*yy, zeros(1,length(xx)),'r--');
    axis equal;
    grid on;
    xlabel('x');
    ylabel('y');
    zlabel('z');
    title(strcat('ROSETTA mission 3D plot (J2000) [printed ', datestr(now), ']'));

    d(1,:)=[2004 02 04];    % launch
    d(2,:)=[2005 03 04];    % earth flyby 1
    d(3,:)=[2007 02 25];    % mars flyby
    d(4,:)=[2007 11 13];    % earth flyby 2
    d(5,:)=[2009 11 13];    % earth flyby 3
    d(6,:)=[2014 05 10];    % comet approach
    d(7,:)=[2015 12 31];    % end of mission
    r1=find(min(abs(datenum(d(1,:))-datenum(dater(1:3,:)')))==abs(datenum(d(1,:))-datenum(dater(1:3,:)'))); % rosetta
    r2=find(min(abs(datenum(d(2,:))-datenum(dater(1:3,:)')))==abs(datenum(d(2,:))-datenum(dater(1:3,:)'))); % indecies
    r3=find(min(abs(datenum(d(3,:))-datenum(dater(1:3,:)')))==abs(datenum(d(3,:))-datenum(dater(1:3,:)'))); % for the 
    r4=find(min(abs(datenum(d(4,:))-datenum(dater(1:3,:)')))==abs(datenum(d(4,:))-datenum(dater(1:3,:)'))); % dates
    r5=find(min(abs(datenum(d(5,:))-datenum(dater(1:3,:)')))==abs(datenum(d(5,:))-datenum(dater(1:3,:)')));
    r6=find(min(abs(datenum(d(6,:))-datenum(dater(1:3,:)')))==abs(datenum(d(6,:))-datenum(dater(1:3,:)')));
    c1=find(min(abs(datenum(d(1,:))-datenum(datec(1:3,:)')))==abs(datenum(d(1,:))-datenum(datec(1:3,:)'))); % comet
    c2=find(min(abs(datenum(d(2,:))-datenum(datec(1:3,:)')))==abs(datenum(d(2,:))-datenum(datec(1:3,:)'))); % indecies
    c3=find(min(abs(datenum(d(3,:))-datenum(datec(1:3,:)')))==abs(datenum(d(3,:))-datenum(datec(1:3,:)'))); % for the
    c4=find(min(abs(datenum(d(4,:))-datenum(datec(1:3,:)')))==abs(datenum(d(4,:))-datenum(datec(1:3,:)'))); % dates
    c5=find(min(abs(datenum(d(5,:))-datenum(datec(1:3,:)')))==abs(datenum(d(5,:))-datenum(datec(1:3,:)')));
    c6=find(min(abs(datenum(d(6,:))-datenum(datec(1:3,:)')))==abs(datenum(d(6,:))-datenum(datec(1:3,:)')));
    c7=find(min(abs(datenum(d(7,:))-datenum(datec(1:3,:)')))==abs(datenum(d(7,:))-datenum(datec(1:3,:)')));

    irs=[r1(1) r2(1) r3(1) r4(1) r5(1)];  % make an array with (the first of) the resp. dates
    ire=[r2(1) r3(1) r4(1) r5(1) r6(1)];
    ics=[c1(1) c2(1) c3(1) c4(1) c5(1) c6(1)];
    ice=[c2(1) c3(1) c4(1) c5(1) c6(1) c7(1)];
    
    figure(fign+1)
    set(gcf,'Position', windowsize); 
    xmin=min(rc_hea(1,:));
    xmax=max(rc_hea(1,:));
    ymin=min(rc_hea(2,:));
    ymax=max(dist_mars);
    subname='abcdef';

    for i=1:6
%        figure(fign+i)
        hold on;
        subplot(2,3,i);
        if i<6        
            plot(rr_hea(1,irs(i):ire(i))/AU, rr_hea(2,irs(i):ire(i))/AU);
        end
        hold on;
        plot(rc_hea(1,ics(i):ice(i))/AU, rc_hea(2,ics(i):ice(i))/AU, 'k');
        hold on;
        if i==6
            plot(rc_hea(1,ice(i))/AU, rc_hea(2,ice(i))/AU, 'ko');
        end
        plot(0,0,'y*');
        % plot mars and earth orbit
        plot(dist_earth*xx/AU, dist_earth*yy/AU,'g--');
        plot(dist_mars*xx/AU, dist_mars*yy/AU, 'r--');
        if i<6
 %           plot(rr_gea(1,ire(i))/AU, rr_gea(2,ire(i))/AU, 'bo');
        end
        axis([xmin xmax ymin ymax]/AU);
        axis equal;
%        grid on;
        xlabel('x');
        ylabel('y');
        title(strcat(subname(i), ') ', datestr(datenum(d(i,:))),' - ', datestr(datenum(d(i+1,:)))));
        grid on;    
    end

    suptitle(strcat('ROSETTA mission plot'));
    
    clear xx yy;   
%    subplot(2,3,6)
%    plot(0,0,'b',0,0,'k',0,0,'g',0,0,'r');
%    legend('path of Rosetta', 'path of the comet', 'earth orbit', 'mars orbit');
%    axis off;
    
case 'l'
    % set the current time span    
    disp('Current time span:');
    disp(strcat('Minimum traj data index :1 :', datestr(jd(1)+timediff)));
    disp(strcat('Start index             :', num2str(istart),' :', datestr(jd(istart)+timediff)));
    disp(strcat('Closest approach index  :', num2str(iclosest), ' :', datestr(jd(iclosest)+timediff)));
    disp(strcat('End index               :', num2str(iend), ' :', datestr(jd(iend)+timediff)));
    disp(strcat('Maximum dtraj ata index :', num2str(length(jd)), ' :', datestr(jd(length(jd))+timediff)));
    istart=input('Enter new value for start index:');
    iend=input('Enter new value for end index:');
    
case 'm'
    % save coordinates
    [yy,mm,temp1,temp2,temp3,temp4]=datevec(jd+timediff);
    dd=time(1,:);
    HH=time(2,:);
    MM=time(3,:);
    X=rr_gse(1,:);
    Y=rr_gse(2,:);
    Z=rr_gse(3,:);
    if att==0
        savedata=[yy; mm; dd; HH; MM; X; Y; Z];
    else
        savedata=[yy; mm; dd; HH; MM; X; Y; Z; ];
    end

    fid=fopen(savename, 'w');

    fprintf(fid, '%04u %02u %02u %02u %02u %13e %13e %13e\n', savedata);
    
    fclose(fid);

    clear yy mm temp1 temp2 temp3 temp4 dd HH MM X Y Z savedata;
    disp(strcat('Data saved in "',savename,'".'));


case 'n'
    % save sunlit and wake data to use with "getswp"
    format long;
    savedata = [(jd+timediff)' eclipse1' eclipse2' wake1' wake2'];
    save sunwake.txt savedata -ascii -tabs -double
    disp('Sunlit/wake for the two probes saved in "sunwake.txt" in the format');
    disp('MATLAB day number, eclipse probe 1, eclipse probe 2, wake probe 1, wake probe 2');
    disp('1 represent eclipse resp. not wake, 2 represent eclipse resp. wake.');
    
   
case 'o'
    % plot "measured data color coded" dots in the long-lat plot as seen from one of the probes
    proben=input('LAP probe number: ');
    relradius=input('High-gain antenna relative radius (1=100%): ');

    load fb1photo;

    switch proben
    case 1
        it=t1;
        iphoto=if1*1e9;
        titlestr=strcat(tit,' - Sun vector as seen from LAP probe #1 (',datestr(jd(istart)+timediff,20),...
            '@', datestr(jd(istart)+timediff,15), '-',datestr(jd(iend)+timediff,20),'@', ...
            datestr(jd(iend)+timediff,15),') [printed ', datestr(now), ']');
    case 2
        it=t2;
        iphoto=if2*1e9;
        titlestr=strcat(tit,' - Sun vector as seen from LAP probe #2 (',datestr(jd(istart)+timediff,20),...
            '@', datestr(jd(istart)+timediff,15), '-',datestr(jd(iend)+timediff,20),'@', ...
            datestr(jd(iend)+timediff,15),') [printed ', datestr(now), ']');
    end

    sphiint=interp1(jd+timediff,sphi,it);
    sthetaint=interp1(jd+timediff,stheta,it);
    ik=1/abs(max(iphoto)-min(iphoto));
    
    figure(fign)
    set(gcf,'Position', windowsize);
    drawrosetta(proben, relradius);

    % just for testing...
%    angleshift=-15;
    angleshift=0;
    
    for i=1:length(sphiint)
        if ~isnan(iphoto(i))
            plot(sphiint(i)+angleshift, sthetaint(i), '.', 'color', ...
                [1-abs(iphoto(i)-min(iphoto))*ik 0 abs(iphoto(i)-min(iphoto))*ik]);
        end
    end
    axis([-180 180 -90 90]);
    grid on;
    set(gca,'XDir','reverse');
    xlabel('azimuth in zx plane (x -> -z)');
    ylabel('elevation above the zx plane');
    title(titlestr);
    disp(strcat('RED=', int2str(min(iphoto)), ' nA, BLUE=', int2str(max(iphoto)), ' nA.'));
%    clear it iphoto sphiint sthetaint ik;

    
    
case 'p'
    figure(fign)
    set(gcf,'Position', [0 0 1275 1275/2]); 
    [mars, col]=imread('totmagfield.png','png');
    newplot;
    set(gca,'YDir','normal');
    image([0 360], [90 -90], mars, 'YData', [90 -90]), colormap(col);
    hold on;
    set(gca,'YDir','normal');
    
    
    plot(longitude,latitude,'r.');
    
    
    %Point out the closest approach, interpolate and plot, needs a better
    %interpolation method!!!!!!!!!!!!
  
    %distinter=interp1(jd(istart:iend),distance(istart:iend), [jd(iclosest-1):0.0000002:jd(iclosest+1)]);
    %latinter=interp1(distance(istart:iend), latitude(istart:iend), min(distinter));
    %longinter=interp1(distance(istart:iend), longitude(istart:iend), min(distinter));
    %plot(longinter,latinter,'r*');
    
    plot((longitude(iclosest)+longitude(iclosest+1))/2,(latitude(iclosest)+latitude(iclosest+1))/2,'rs'); %Here is an ugly adjustment!!!!!!
    tempo=strcat('Closest approach: ', num2str(round(distance(iclosest)-radius),'%4.0f'), ' km ',' @ ', datestr(jd(iclosest)+timediff,15));
    text((longitude(iclosest)+longitude(iclosest-1))/2-12,(latitude(iclosest)+latitude(iclosest-1))/2+6, tempo,'FontWeight','bold'); %Here too!!!!!!!!!
    
    %Plot start and stop, which is hard to see otherwise...
    h=text(longitude(istart),latitude(istart)+3,'Start','Fontweight', 'bold');
    h=text(longitude(iend),latitude(iend)+3,'End','Fontweight','bold');
    set(h,'color','k');
     
    %Plot out 10000km and time on the way in
    for e=istart:iend-1
        if distance(e) >= (10000+radius) & distance(e+1) <= (10000+radius)
        ten=e;
        plot((longitude(e)+longitude(e+1))/2,(latitude(e)+latitude(e+1))/2,'r*');
        tempo=strcat('10000 km@ ', datestr(jd(ten)+timediff,15));
        text((longitude(e)+longitude(e+1))/2,(latitude(e)+latitude(e+1))/2+3,tempo ,'FontWeight','bold');
        
        end
    end
    
    %Interpolate some nice distances on the way in...
    for j=[300 500 1000 3000 5000];
        timesinter=interp1(distance(istart:iclosest),jd(istart:iclosest),j+radius);
        latinter(j)=interp1(jd(istart:iclosest),latitude(istart:iclosest),timesinter);
        longinter(j)=interp1(jd(istart:iclosest),longitude(istart:iclosest),timesinter);
        plot(longinter(j),latinter(j),'r*');
        avst=strcat(num2str(j));
        text(longinter(j),latinter(j)+3,avst,'FontWeight','bold');  
    end
    
    
        %Interpolate some nice distances on the way out...
    for j=[300 500 1000 3000 5000];
        timesinter=interp1(distance(iclosest:iend),jd(iclosest:iend),j+radius);
        latinter(j)=interp1(jd(iclosest:iend),latitude(iclosest:iend),timesinter);
        longinter(j)=interp1(jd(iclosest:iend),longitude(iclosest:iend),timesinter);
        plot(longinter(j),latinter(j),'r*');
        avst=strcat(num2str(j));
        text(longinter(j),latinter(j)+3,avst,'FontWeight','bold');  
    end
    

    %Plot out 10000km and the time on the way out
    for e=istart:iend-1
        if distance(e+1) >= (10000+radius) & distance(e) <= (10000+radius)
        ten=e
        plot((longitude(e)+longitude(e+1))/2,(latitude(e)+latitude(e+1))/2,'r*');
        tempo=strcat('10000 km@ ', datestr(jd(ten)+timediff,15));
        text((longitude(e)+longitude(e+1))/2,(latitude(e)+latitude(e+1))/2-3,tempo ,'FontWeight','bold');
        
        end
    end
    
    title('The Rosetta trajectory projected on a map of the Martian total magnetic field (nT), 24-feb-2007 - 25-feb-2007, Magnetic map acquired using Jafar Arkani-Hameds model')
    xlabel('East longitude');
    ylabel('Latitude');

case 'q'
    figure(fign)
    set(gcf,'Position', [0 0 655 472]); 
    [phopos2, col]=imread('phobos2.png','png'); %Mars magnetosphere picture from ASPERA on Phobos-2, from Lundin and Barabash 2004 
    newplot;
    set(gca,'YDir','normal');
    image([-25.5 5], [10 -10], phopos2, 'YData', [10 -10]), colormap(col);
    hold on;
    set(gca,'YDir','normal');
    
    %Plotta ut banan
    plot(x,y,'r');
    
    %Plotta ut closest approach
    plot(x(iclosest),y(iclosest),'bs');
    temp{1,1}=strcat('Closest approach: ', num2str(round(distance(iclosest)-radius),'%3.1f'),'km');
    temp{2,1}=strcat('@', datestr(jd(iclosest)+timediff,15));
    h=text( x(iclosest), y(iclosest)+1,temp, 'FontWeight', 'bold');
    set(h, 'Color', 'b');
    
    for w=1050:15:1215 %Prints out the times
        plot(x(w),y(w),'bs');
        tid=strcat(datestr(jd(w)+timediff,15));
        text(x(w)-1.8,y(w),tid,'FontWeight','bold');
    end
    
    
    title('The Rosetta trajectory projected on a map of the Martian magnetosphere, 24-feb-2007 - 25-feb-2007')
    xlabel('x, positive x towards sun, (Mars radii)');
    ylabel('y (Mars radii)');
    
    case 'r'
        figure(fign)
        
        subplot(2,2,1)
        plot(rpho_mrx(1,:)./radius,rpho_mrx(2,:)./radius,'r')
        hold on
        plot(rdei_mrx(1,:)./radius,rdei_mrx(2,:)./radius,'g')
        plot(rr_mrx(1,:)./radius,rr_mrx(2,:)./radius);
        %quiver(0,0,0,0)
        
        axis([-10 10 -10 10])
        title('x-y plane')
        
        subplot(2,2,2)
        plot(rpho_mrx(1,:)./radius,rpho_mrx(3,:)./radius,'r')
        hold on
        plot(rdei_mrx(1,:)./radius,rdei_mrx(3,:)./radius,'g')
        plot(rr_mrx(1,:)./radius,rr_mrx(3,:)./radius);
        quiver(0,0,0,3);
        axis([-10 10 -10 10])
        title('x-z plane')
        
        subplot(2,2,3)
        plot(rpho_mrx(2,:)./radius,rpho_mrx(3,:)./radius,'r')
        hold on
        plot(rdei_mrx(2,:)./radius,rdei_mrx(3,:)./radius,'g')
        plot(rr_mrx(2,:)./radius,rr_mrx(3,:)./radius)
        quiver(0,0,0,3);
        axis([-10 10 -10 10])
        title('y-z plane')   
end