%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: plotDSMCPolarOrbit.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]
%         v     - In orbit velocity
%
% Output: -
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotDSMCPolarOrbit(dataD,elem,v)

global CO;

% 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,:)];
end

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

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).*10^-6;

[theta,phi,r] = cart2sph(Dx,Dy,Dz);
[Rnx,Rny,Rnz] = sph2cart(theta,phi,Dn);
[Rnx,Rny,Rnz] = ecl2orb(elem,Rnx,Rny,Rnz);
[THn,Rn] = cart2pol(Rnx,Rny);

subplot(2,2,1), polar(THn,Rn,'b')
title('Neutral density (cm^{-3})')
xlabel('Projected x-axis')
ylabel('Projected y-axis')


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

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

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

subplot(2,2,2), polar(THu,Ru)
hold on
title('Neutral Velocity (km/s)')
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).*CO.Kb2;

[theta,phi,r] = cart2sph(Dx,Dy,Dz);
[RTx,RTy,RTz] = sph2cart(theta,phi,DT);
[RTx,RTy,RTz] = ecl2orb(elem,RTx,RTy,RTz);
[THT,RT] = cart2pol(RTx,RTy);

subplot(2,2,3), polar(THT,RT,'b')
title('Neutral temperature (eV)')
xlabel('Projected x-axis')
ylabel('Projected y-axis')


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

radius = num2str(round(mean(r))/1000);
O = num2str(omega);
i = num2str(inc);
v = num2str(v);
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)
subplot(2,2,4), text(0.1,0.1,'S/C speed:','fontsize',16)
subplot(2,2,4), text(0.6,0.1,v,'fontsize',16)
subplot(2,2,4), text(0.95,0.1,'m/s','fontsize',16)
axis off
