% Displays a number of plots over the flyby
%
% Initial version: Magnus Billvik, Swedish Institute of Space Physics,
%Uppsala, June 2005

%global number nlstart date jd mjd date datea dater timediff inclination radiuse theta;
%global RR VR RM VM RS VS;
%global rr_gei vr_gei rs_gei vs_gei rm_gei vm_gei;
%global rr_gea vr_gea rs_gea vs_gea rm_gea vm_gea;
%global rr_gse vr_gse rs_gse vs_gse rm_gse vm_gse;
%global rr_geo vr_geo rs_geo vs_geo rm_geo vm_geo;
%%global vr_ros re_ros rs_ros rm_ros;
%global longitude latitude altitude;
%global Q q1 q2 q3 q4;
%global Asc_gea Asc_gse;
%global srho sphi stheta;
%global mrho mphi mtheta;
%global vrho vphi vtheta;
%global ux uy uz distx disty sistz dist_ur dotp;
%global anghor;

%fb=2882;    % index for closest approach (earth)
%ua=2868;    % index for closest approach (uppsala)
%nstart=ua-60;
%nend=ua+20;
%nstart=ua-2*60;
%nend=ua+60;

% uppsla's coordinates
lat_u=(90-59.86)*pi/180;
long_u=17.64*pi/180;
ux=radiuse*sin(lat_u)*cos(long_u);
uy=radiuse*sin(lat_u)*sin(long_u);
uz=radiuse*cos(lat_u);
dist_u=sqrt(ux.^2+uy.^2+uz.^2);

rx=rr_geo(1,:)*radiuse;
ry=rr_geo(2,:)*radiuse;
rz=rr_geo(3,:)*radiuse;

% the vector from uppsala to rosetta
distx=rx-ux;
disty=ry-uy;
distz=rz-uz;
dist_ur=sqrt(distx.^2+disty.^2+distz.^2);

% the vector pointing east from Uppsala
long_u_east=long_u+pi/2;   % +90 degrees longitude
east(1)=sin(lat_u)*cos(long_u_east);
east(2)=sin(lat_u)*sin(long_u_east);
east(3)=cos(lat_u);


% calc the dot product of the uppsala xyz-vector and the uppsala-rosetta xyz-vector
for i=1:length(distx)
    dotp(i)=dot([ux uy uz], [distx(i) disty(i) distz(i)]);
end

% calculate the angle above the horizon from the dot product above
anghor=90-(180/pi)*acos(dotp./(dist_u*dist_ur));
disp(size(anghor));

for i=1:length(distx)
    % calc the cross product of the uppsala vector and the uppsala-rosetta vector
    crossprod=cross([ux uy uz], [distx(i) disty(i) distz(i)]);

    % calculate the azimuth 
    % = the angle between the south vector and rosettas projection on the horizontal plane 
    % = the angle between the east vector and the cross prod above
    azimuth(i)=(180/pi)*acos(dot(east, crossprod)/(sqrt(east*east')*sqrt(crossprod*crossprod')))-90;
end

fb=iclosest;
ua=find(dist_ur==min(dist_ur));
nstart=ua-40;
nend=ua+40;
radiuse=rad_earth;

figure(1);
set(gcf,'Position',[0 0 1275 945]); % maximizes the window size

subplot(3,1,1);
hold on;
plot(jd([nstart:nend])+timediff, anghor([nstart:nend]));
plot(jd([nstart nend])+timediff, [0 0], 'g-');
plot(jd(fb)+timediff, anghor(fb), 'ro');
text(jd(fb)+timediff, anghor(fb), 'Closest approach');
plot(jd(ua)+timediff, anghor(ua), 'ro');
text(jd(ua)+timediff, anghor(ua), 'Closest Uppsala');
grid on;
title('Rosetta Earth flyby 4/3-05');
xlabel('time (UT)');
ylabel('Rosettas angle over the horizon as seen from Uppsala (deg)');
datetick('x',15);

subplot(3,1,2)
hold on;
plot(jd([nstart:nend])+timediff,dist_ur([nstart:nend]));
grid on;
title('Rosetta Earth flyby 4/3-05');
xlabel('time (UT)');
ylabel('distance Uppsala - Rosetta (km)');
datetick('x',15);

%subplot(3,1,3)
%hold on;
%plot(jd([nstart:nend])+timediff,altitude([nstart:nend]));
%grid on;
%title('Rosetta Earth flyby 4/3-05');
%xlabel('time (UT)');
%ylabel('altitude from earth surface');
%datetick('x',15);

subplot(3,1,3)
hold on;
plot(jd([nstart:nend])+timediff,azimuth([nstart:nend]));
grid on;
title('Rosetta Earth flyby 4/3-05');
xlabel('time (UT)');
ylabel('azimuth (degrees from south towards west)');
datetick('x',15);


mindist=min(dist_ur);
minindex=find(dist_ur==mindist);
format short g;
disp('Minimum distance (km) Uppsala-Rosetta:');
disp(mindist);
disp('Index for min. distance:');
disp(minindex);
disp('Date and time (GMT) for minimum distance:');
date(:,minindex)';
