%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: plotHybridOrbit.m
% Author: J. Kristoffer Hultgren
%
% Description:
%
% Plots the plasma environment in a given orbit at a given heliocentric
% distance
%
% Input:  dataH - Hybrid model data
%         elem  - Orbital elements [r inc omega]
%
% Output: -
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotHybridOrbit(dataH,elem)

Hx = dataH(:,1);
Hy = dataH(:,2);
Hz = dataH(:,3);

omega = elem(3);
inc = elem(2);

% Rotation matrix
R = [cos(omega) -sin(omega)*cos(inc) sin(omega)*sin(inc); 
        sin(omega) cos(omega)*cos(inc) -cos(omega)*sin(inc); 
        0 sin(inc) cos(inc)];

% Performs the rotation
Result = inv(R)*[Hx';Hy';Hz'];

x2 = Result(1,:);
y2 = Result(2,:);
z2 = Result(3,:);

% Checks for transformation errors
if max(z2)>1
    disp('Transformation failure')
    return
end

scrsz = get(0,'ScreenSize');
figure('name','Hybrid','Position',[1 scrsz(4)/2-60 scrsz(3)/3 scrsz(4)/2])


%Plot SW & HI densities-----------------------------------------

Hrhosw = dataH(:,11);
Hrhohi = dataH(:,15);

[theta,phi,r] = cart2sph(Hx,Hy,Hz);
Rsw = r + (Hrhosw)./100;
[Rswx,Rswy,Rswz] = sph2cart(theta,phi,Rsw);

[Rswx,Rswy,Rswz] = ecl2orb(elem,Rswx,Rswy,Rswz);

[theta,phi,r] = cart2sph(Hx,Hy,Hz);
Rhi = r + (Hrhohi)./500;
[Rhix,Rhiy,Rhiz] = sph2cart(theta,phi,Rhi);

[Rhix,Rhiy,Rhiz] = ecl2orb(elem,Rhix,Rhiy,Rhiz);

subplot(2,2,1), plot(x2,y2,'k')
hold on
subplot(2,2,1), plot(Rswx,Rswy,'b')
subplot(2,2,1), plot(Rhix,Rhiy,'g')
title('Solar wind & ion density (1:100 & 1:500)')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')
legend('Orbit','SW-density','HI-density')


%Plot HI-field--------------------------------------------------

Huhix = dataH(:,16);
Huhiy = dataH(:,17);
Huhiz = dataH(:,18);
magUHI = sqrt(Huhix.^2+Huhiy.^2+Huhiz.^2);
[theta,phi,r] = cart2sph(Hx,Hy,Hz);
Ruhi = r + magUHI;
[Ruhix,Ruhiy,Ruhiz] = sph2cart(theta,phi,Ruhi);

[Ruhix,Ruhiy,Ruhiz] = ecl2orb(elem,Ruhix,Ruhiy,Ruhiz);
[Huhix,Huhiy,Huhiz] = ecl2orb(elem,Huhix,Huhiy,Huhiz);

subplot(2,2,2), plot(x2,y2,'k')
hold on
subplot(2,2,2), plot(Ruhix,Ruhiy)
hold on
title('Heavy Ion Velocity')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')

new = zeros(1,max(size(Huhix))/10);
m = 1;
while m<max(size(Huhix))/10
    new(m) = Huhix(10*(m-1)+1);
    m = m+1;
end
Huhix = new;

new = zeros(1,max(size(Huhiy))/10);
m = 1;
while m<max(size(Huhiy))/10
    new(m) = Huhiy(10*(m-1)+1);
    m = m+1;
end
Huhiy = new;

new = zeros(1,max(size(Ruhix))/10);
m = 1;
while m<max(size(Ruhix))/10
    new(m) = Ruhix(10*(m-1)+1);
    m = m+1;
end
Ruhix = new;

new = zeros(1,max(size(Ruhiy))/10);
m = 1;
while m<max(size(Ruhiy))/10
    new(m) = Ruhiy(10*(m-1)+1);
    m = m+1;
end
Ruhiy = new;

subplot(2,2,2), quiver(Ruhix,Ruhiy,Huhix,Huhiy,'b')


%Plot SW-field--------------------------------------------------

Huswx = dataH(:,12);
Huswy = dataH(:,13);
Huswz = dataH(:,14);
magUsw = sqrt(Huswx.^2+Huswy.^2+Huswz.^2)/20;
[theta,phi,r] = cart2sph(Hx,Hy,Hz);
Rusw = r + magUsw;
[Ruswx,Ruswy,Ruswz] = sph2cart(theta,phi,Rusw);
[Ruswx,Ruswy,Ruswz] = ecl2orb(elem,Ruswx,Ruswy,Ruswz);
[Huswx,Huswy,Huswz] = ecl2orb(elem,Huswx,Huswy,Huswz);

subplot(2,2,3), plot(x2,y2,'k')
hold on
subplot(2,2,3), plot(Ruswx,Ruswy)
hold on
title('Solar Wind Velocity (1:20)')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')

new = zeros(1,max(size(Huswx))/10);
m = 1;
while m<max(size(Huswx))/10
    new(m) = Huswx(10*(m-1)+1);
    m = m+1;
end
Huswx = new;

new = zeros(1,max(size(Huswy))/10);
m = 1;
while m<max(size(Huswy))/10
    new(m) = Huswy(10*(m-1)+1);
    m = m+1;
end
Huswy = new;

new = zeros(1,max(size(Ruswx))/10);
m = 1;
while m<max(size(Ruswx))/10
    new(m) = Ruswx(10*(m-1)+1);
    m = m+1;
end
Ruswx = new;

new = zeros(1,max(size(Ruswy))/10);
m = 1;
while m<max(size(Ruswy))/10
    new(m) = Ruswy(10*(m-1)+1);
    m = m+1;
end
Ruswy = new;

subplot(2,2,3), quiver(Ruswx,Ruswy,Huswx,Huswy,'b')


%Plot info----------------------------------------

radius = num2str(round(mean(r))/1000);
O = num2str(omega);
i = num2str(inc);
sw = sign(dataH(1,1)/dataH(1,12));
if sw < 0
    out = 'Sun is in positive x-direction';
end
if sw > 0
    out = 'Sun is in negative x-direction';
end
subplot(2,2,4), text(0.1,0.9,'\bf{Hybrid}','fontsize',16)
subplot(2,2,4), text(0.1,0.7,'\it{r}  =','fontsize',16)
subplot(2,2,4), text(0.7,0.7,'(km)','fontsize',16)
subplot(2,2,4), text(0.25,0.7,radius,'fontsize',16)
subplot(2,2,4), text(0.1,0.5,'\Omega  \it{=}','fontsize',16)
subplot(2,2,4), text(0.7,0.5,'(rad)','fontsize',16)
subplot(2,2,4), text(0.3,0.5,O,'fontsize',16)
subplot(2,2,4), text(0.1,0.3,'\it{i}  =','fontsize',16)
subplot(2,2,4), text(0.7,0.3,'(rad)','fontsize',16)
subplot(2,2,4), text(0.25,0.3,i,'fontsize',16)
subplot(2,2,4), text(0.1,0.1,out,'fontsize',16)
axis off
