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

function plotMHDOrbit(dataM,elem)

Mx = dataM(:,1);
My = dataM(:,2);
Mz = dataM(:,3);

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

% Defines the 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)*[Mx';My';Mz'];

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','MHD','Position',[scrsz(3)/3 scrsz(4)/2-60 scrsz(3)/3 scrsz(4)/2])

%Plot density-----------------------------------------

Mn = dataM(:,5);

[theta,phi,r] = cart2sph(Mx,My,Mz);
Rn = r + (Mn)./100;
[Rnx,Rny,Rnz] = sph2cart(theta,phi,Rn);

[Rnx,Rny,Rnz] = ecl2orb(elem,Rnx,Rny,Rnz);

subplot(2,2,1), plot(x2,y2,'k')
hold on
subplot(2,2,1), plot(Rnx,Rny,'b')
title('Fluid density (1:100)')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')
legend('Orbit','Density')

%Plot pressure-----------------------------------------

Mp = dataM(:,12);

[theta,phi,r] = cart2sph(Mx,My,Mz);
Rp = r + (Mp).*10^15;
[Rpx,Rpy,Rpz] = sph2cart(theta,phi,Rp);

[Rpx,Rpy,Rpz] = ecl2orb(elem,Rpx,Rpy,Rpz);

subplot(2,2,2), plot(x2,y2,'k')
hold on
subplot(2,2,2), plot(Rpx,Rpy,'b')
title('Fluid pressure (1:10^{15})')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')

%Plot Velocity field--------------------------------------------------

Mux = dataM(:,6);
Muy = dataM(:,7);
Muz = dataM(:,8);
magU = sqrt(Mux.^2+Muy.^2+Muz.^2);
[theta,phi,r] = cart2sph(Mx,My,Mz);
Ru = r + magU./100;
[Rux,Ruy,Ruz] = sph2cart(theta,phi,Ru);

[Rux,Ruy,Ruz] = ecl2orb(elem,Rux,Ruy,Ruz);
[Mux,Muy,Muz] = ecl2orb(elem,Mux,Muy,Muz);

subplot(2,2,3), plot(x2,y2,'k')
hold on
subplot(2,2,3), plot(Rux,Ruy)
hold on
title('Fluid Velocity (1:100)')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')

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

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

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

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

subplot(2,2,3), quiver(Rux,Ruy,Mux,Muy,'b')

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

radius = num2str(round(mean(r))/1000);
O = num2str(omega);
i = num2str(inc);
sw = sign(dataM(1,1)/dataM(1,6));
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{MHD}','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
