%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: plotDSMCPlane.m
% Author: J. Kristoffer Hultgren
%
% Description:
%
% Plots the plasma environment in a given plane at a given heliocentric
% distance
%
% Input:  dataD - DSMC model data
%         heli  - Heliocentric distance
%
% Output: -
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotDSMCPlane(dataD,heli)

% 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


span = round(sqrt(max(size(dataD))));
rangeL = (span - 1)/2;
if rangeL > 5
    spanS = 11;    
else
    spanS = span;
end

scale = span / (2 * spanS);

Dx = dataD(:,1).*10^-3;
Dy = dataD(:,2).*10^-3;
Dz = dataD(:,3).*10^-3;

coo = max([max(Dx) max(Dy) max(Dz)]);

x = linspace(-coo,coo,span);

DSMC1 = figure('name','DSMC1');
dist = heli;


if max(Dx) == 0
           
    %----Velocity----

    %Plot speed    
    DVx = splitVector(dataD(:,5), span);
    DVy = splitVector(dataD(:,6), span);
    DVz = splitVector(dataD(:,7), span);
    speed = (DVx.^2+DVy.^2+DVz.^2).^(1/2);
    subplot(2,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Neutral gas velocity (km/s)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    %Plot field direction
    DVy = scaleVec(DVy, span, spanS);
    DVz = scaleVec(DVz, span, spanS);
    quiver(x,x,DVy,DVz,scale,'k')
    axis square
    hold on
    

    %----Temperature----
    
    DT = splitVector(dataD(:,8), span);    
    subplot(1,2,2), pcolor(x,x,DT)
    shading interp
    axis square
    hold on
    colorbar
    title('Neutral temperature (K)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    
    %----Density----
    Dn = splitVector(dataD(:,9), span).*1e-6;
    
    titel = 'Neutral density (cm^{-3})';
    maxx = max(max(Dn));
    meann = mean(mean(Dn));
    if (maxx - meann) > (3*maxx/4)
        Dn = log10(Dn);
        titel = 'Neutral density (log_{10}(cm^{-3}))';
    end
    
    subplot(1,2,1), pcolor(x,x,Dn)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
end

if max(Dy) == 0
    
    %----Velocity----

    %Plot speed    
    DVx = splitVector(dataD(:,5), span);
    DVy = splitVector(dataD(:,6), span);
    DVz = splitVector(dataD(:,7), span);
    speed = (DVx.^2+DVy.^2+DVz.^2).^(1/2);
    subplot(2,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Neutral gas velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')
    
    %Plot field direction
    DVx = scaleVec(DVx, span, spanS);
    DVz = scaleVec(DVz, span, spanS);
    quiver(x,x,DVx,DVz,scale,'k')
    axis square
    hold on
       
    
    %----Temperature----
    DT = splitVector(dataD(:,8), span);    
    subplot(1,2,2), pcolor(x,x,DT)
    shading interp
    axis square
    hold on
    colorbar
    title('Neutral temperature (K)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')



    %----Density----
    Dn = splitVector(dataD(:,9), span).*1e-6;
    
    titel = 'Neutral density (cm^{-3})';
    maxx = max(max(Dn));
    meann = mean(mean(Dn));
    if (maxx - meann) > (3*maxx/4)
        Dn = log10(Dn);
        titel = 'Neutral density (log_{10}(cm^{-3}))';
    end
    
    subplot(1,2,1), pcolor(x,x,Dn)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')
end

if max(Dz) == 0
    %----Velocity----

    %Plot speed    
    DVx = splitVector(dataD(:,5), span);
    DVy = splitVector(dataD(:,6), span);
    DVz = splitVector(dataD(:,7), span);
    speed = (DVx.^2+DVy.^2+DVz.^2).^(1/2);
    subplot(2,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Neutral gas velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')
    
    %Plot field direction
    DVx = scaleVec(DVx, span, spanS);
    DVy = scaleVec(DVy, span, spanS);
    quiver(x,x,DVx,DVy,scale,'k')
    axis square
    hold on
       
    
    %----Temperature----
    DT = splitVector(dataD(:,8), span);    
    subplot(1,2,2), pcolor(x,x,DT)
    shading interp
    axis square
    hold on
    colorbar
    title('Neutral temperature (K)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')



    %----Density----
    Dn = splitVector(dataD(:,9), span).*1e-6;
    
    titel = 'Neutral density (cm^{-3})';
    maxx = max(max(Dn));
    meann = mean(mean(Dn));
    if (maxx - meann) > (3*maxx/4)
        Dn = log10(Dn);
        titel = 'Neutral density (log_{10}(cm^{-3}))';
    end
    
    subplot(1,2,1), pcolor(x,x,Dn)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')
end

reply = 'n';
reply = input('Do you want to print? (y/n): ','s');

if isempty(reply)
    reply = 'n';
end

if reply == 'y'
    firstpage = figure('name','firstpage');
    text(0.1,0.9,'\bf{Model: DSMC}','fontsize',16)
    text(0.1,0.7,'\bf{Distance: }','fontsize',16)
    text(0.3,0.7,dist,'fontsize',16)    
    axis off
    saveas(firstpage,'1st','jpg')
    saveas(DSMC1,'DSMC1','jpg')
    saveas(DSMC2,'DSMC2','jpg')  
    !lpr -P hp4600 1st.jpg
    !lpr -P hp4600 DSMC1.jpg
    !lpr -P hp4600 DSMC2.jpg
    !rm 1st.jpg
    !rm DSMC1.jpg
    !rm DSMC2.jpg
end