%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: plotMHDPolarOrbit.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 plotMHDPolarOrbit(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,:);

kJ = 1.38*10^-23;   % J/K
keV =8.62*10^-5;    % eV/K

% 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 B-field------------------------------------------
MBx = dataM(:,9).*1e9;
MBy = dataM(:,10).*1e9;
MBz = dataM(:,11).*1e9;
magB = sqrt(MBx.^2+MBy.^2+MBz.^2);
[theta,phi,r] = cart2sph(Mx,My,Mz);
[RBx,RBy,RBz] = sph2cart(theta,phi,magB);
[RBx,RBy,RBz] = ecl2orb(elem,RBx,RBy,RBz);
[THB,RB] = cart2pol(RBx,RBy);

[MBx,MBy,MBz] = ecl2orb(elem,MBx,MBy,MBz);

subplot(2,2,2), polar(THB,RB)
hold on
title('B-field (nT)')

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

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

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

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

%quiver(RBx,RBy,MBx,MBy,'b')




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

Mn = dataM(:,5).*10^-6;

[theta,phi,r] = cart2sph(Mx,My,Mz);
[Rnx,Rny,Rnz] = sph2cart(theta,phi,Mn);
[Rnx,Rny,Rnz] = ecl2orb(elem,Rnx,Rny,Rnz);
[THn,Rn] = cart2pol(Rnx,Rny);

subplot(2,2,1), polar(THn,Rn,'b')
title('Fluid density (cm^{-3})')


%Plot temperature-----------------------------------------

MT = dataM(:,12)./(kJ.*dataM(:,5)).*keV;

[theta,phi,r] = cart2sph(Mx,My,Mz);
[Rpx,Rpy,Rpz] = sph2cart(theta,phi,MT);
[Rpx,Rpy,Rpz] = ecl2orb(elem,Rpx,Rpy,Rpz);
[THp,Rp] = cart2pol(Rpx,Rpy);

%subplot(2,2,2), polar(THp,Rp,'b')
%title('Fluid temperature (eV)')

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

Mux = dataM(:,6);
Muy = dataM(:,7);
Muz = dataM(:,8);
magU = sqrt(Mux.^2+Muy.^2+Muz.^2).*10^-3;
[theta,phi,r] = cart2sph(Mx,My,Mz);
[Rux,Ruy,Ruz] = sph2cart(theta,phi,magU);
[Rux,Ruy,Ruz] = ecl2orb(elem,Rux,Ruy,Ruz);
[THu,Ru] = cart2pol(Rux,Ruy);

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

subplot(2,2,3), polar(THu,Ru)
hold on
title('Fluid Velocity (km/s)')

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
