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

function plotDSMCOrbit(dataD,elem)

% Controls size of data set (occationally the data set recieved for DSMC
%                               was missing the last line)
nCompare = numel(dataD(:,1));
half = floor(nCompare/2);
if 2*half ~= nCompare
    dataD = [dataD;dataD(nCompare,:)];  % Duplicates last row
end

Dx = dataD(:,1);
Dy = dataD(:,2);
Dz = dataD(:,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);];

Result = inv(R)*[Dx';Dy';Dz'];      % Performs the rotation

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


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

Dn = dataD(:,9);

[theta,phi,r] = cart2sph(Dx,Dy,Dz);
Rn = r + (Dn).*10^-8;
[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('Neutral density (1:10^{-8})')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')
legend('Orbit','Density')


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

Dux = dataD(:,5);
Duy = dataD(:,6);
Duz = dataD(:,7);
magU = sqrt(Dux.^2+Duy.^2+Duz.^2);
[theta,phi,r] = cart2sph(Dx,Dy,Dz);
Ru = r + magU./1;
[Rux,Ruy,Ruz] = sph2cart(theta,phi,Ru);

[Rux,Ruy,Ruz] = ecl2orb(elem,Rux,Ruy,Ruz);
[Dux,Duy,Duz] = ecl2orb(elem,Dux,Duy,Duz);

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

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

new = zeros(1,max(size(Duy))/10);
m = 1;
while m<max(size(Duy))/10
    new(m) = Duy(10*(m-1)+1);
    m = m+1;
end
Duy = 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,2), quiver(Rux,Ruy,Dux,Duy,'b')


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

DT = dataD(:,8);

[theta,phi,r] = cart2sph(Dx,Dy,Dz);
RT = r + (DT).*100;
[RTx,RTy,RTz] = sph2cart(theta,phi,RT);

[RTx,RTy,RTz] = ecl2orb(elem,RTx,RTy,RTz);

subplot(2,2,3), plot(x2,y2,'k')
hold on
subplot(2,2,3), plot(RTx,RTy,'b')
title('Neutral temperature (1:100)')
axis equal
xlabel('Projected x-axis')
ylabel('Projected y-axis')
legend('Orbit','Temperature')


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

radius = num2str(round(mean(r))/1000);
O = num2str(omega);
i = num2str(inc);
subplot(2,2,4), text(0.1,0.9,'\bf{DSMC}','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)
axis off
