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

function plotHybridPlane(dataH,heli)

% Controls size of data set (occationally the data set recieved for DSMC
%                               was missing the last line)
nCompare = numel(dataH(:,1));
half = floor(nCompare/2);
if 2*half == nCompare
    dataH = [dataH;dataH(nCompare,:)];
end
   
span = round(sqrt(max(size(dataH))));
rangeL = (span - 1)/2;
if rangeL > 5
    spanS = 11;    
else
    spanS = span;
end

scale = span / (2 * spanS);

Hx = dataH(:,1).*10^-3;
Hy = dataH(:,2).*10^-3;
Hz = dataH(:,3).*10^-3;

coo = max([max(Hx) max(Hy) max(Hz)]);

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

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

if max(Hx) == 0

    %----Solar Wind Density----
        
    Hrhosw = splitVector(dataH(:,11), span).*10^-6;    
    subplot(2,2,1), pcolor(x,x,Hrhosw)
    shading interp
    axis square
    hold on
    colorbar
    title('SW-density (cm^{-3})')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    
    %----Heavy Ion Density----
    
    Hrhohi = splitVector(dataH(:,15), span).*10^-6;
    
    titel = 'HI-density (cm^{-3})';
    maxx = max(max(Hrhohi));
    meann = mean(mean(Hrhohi));
    if (maxx - meann) > (3*maxx/4)
        Hrhohi = log10(Hrhohi);
        titel = 'HI-density (log_{10}(cm^{-3}))';
    end
    
    subplot(2,2,2), pcolor(x,x,Hrhohi)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')

    
    %----Solar Wind Velocity----

    %Plot speed
    Huswx = splitVector(dataH(:,12), span);
    Huswy = splitVector(dataH(:,13), span);
    Huswz = splitVector(dataH(:,14), span);
    speed = (Huswx.^2+Huswy.^2+Huswz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('SW-velocity (km/s)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    %Scale and plot velocity arrows
    Huswy = scaleVec(Huswy, span, spanS);
    Huswz = scaleVec(Huswz, span, spanS);
    quiver(x,x,Huswy,Huswz,scale,'k')
    axis square
    hold on

    
    %----Heavy Ion Velocity----

    %Plot speed
    Huhix = splitVector(dataH(:,16), span);    
    Huhiy = splitVector(dataH(:,17), span);
    Huhiz = splitVector(dataH(:,18), span);
    speed = (Huhix.^2+Huhiy.^2+Huhiz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,4), pcolor(x,x, speed)
    shading interp
    axis square
    hold on
    colorbar
    title('HI-velocity (km/s)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')

    %Scale and plot velocity arrows
    Huhiy = scaleVec(Huhiy, span, spanS);
    Huhiz = scaleVec(Huhiz, span, spanS);
    quiver(x,x,Huhiy,Huhiz,scale,'k')
    axis square
    hold on
    

    Hybrid2 = figure('name','Hybrid2');

    %----Magnetic Field----

    %Plot field strength
    HBx = splitVector(dataH(:,5), span);
    HBy = splitVector(dataH(:,6), span);
    HBz = splitVector(dataH(:,7), span);
    speed = 10^9.*((HBx.^2+HBy.^2+HBz.^2).^(1/2));
    subplot(2,2,1), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('B-field (nT)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    %Plot field direction
    HBy = scaleVec(HBy, span, spanS);
    HBz = scaleVec(HBz, span, spanS);
    quiver(x,x,HBy,HBz,scale,'k')
    axis square
    hold on


    %----Electric field----

    %Plot field strength
    HEx = splitVector(dataH(:,8), span);
    HEy = splitVector(dataH(:,9), span);
    HEz = splitVector(dataH(:,10), span);
    speed = (HEx.^2+HEy.^2+HEz.^2).^(1/2);
    subplot(2,2,2), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('E-field (V/m)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    %Plot field direction
    HEy = scaleVec(HEy, span, spanS);
    HEz = scaleVec(HEz, span, spanS);
    quiver(x,x,HEy,HEz,scale,'k')
    axis square
    hold on


    %----Total density----
    
    HrhoTot = dataH(:,11) + dataH(:,15);
    HrhoTot = splitVector(HrhoTot, span);
    
    titel = 'Total density (cm^{-3})';
    maxx = max(max(Hrhohi));
    meann = mean(mean(Hrhohi));
    if (maxx - meann) > (3*maxx/4)
        HrhoTot = log10(HrhoTot);
        titel = 'Total density (log_{10}(cm^{-3}))';
    end
    
    subplot(2,2,3), pcolor(x,x,HrhoTot.*10^-6)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    
    %----Total Velocity----

    %Plot speed
    HuTotx = (dataH(:,11).*dataH(:,12) + dataH(:,15).*dataH(:,16))./(dataH(:,11) + dataH(:,15));
    HuTotx = splitVector(HuTotx, span);
    HuToty = (dataH(:,11).*dataH(:,13) + dataH(:,15).*dataH(:,17))./(dataH(:,11) + dataH(:,15));
    HuToty = splitVector(HuToty, span);
    HuTotz = (dataH(:,11).*dataH(:,14) + dataH(:,15).*dataH(:,18))./(dataH(:,11) + dataH(:,15));
    HuTotz = splitVector(HuTotz, span);
    speed = (HuTotx.^2+HuToty.^2+HuTotz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,4), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Total velocity (km/s)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')

    %Scale and plot velocity arrows
    HuToty = scaleVec(HuToty, span, spanS);
    HuTotz = scaleVec(HuTotz, span, spanS);
    quiver(x,x,HuToty,HuTotz,scale,'k')
    axis square
    hold on
    
end

if max(Hy) == 0
    
    %----Solar Wind Density----
        
    Hrhosw = splitVector(dataH(:,11), span).*10^-6;    
    subplot(2,2,1), pcolor(x,x,Hrhosw)
    shading interp
    axis square
    hold on
    colorbar
    title('SW-density (cm^{-3})')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')
    
    
    %----Heavy Ion Density----
    
    Hrhohi = splitVector(dataH(:,15), span).*10^-6;
    
    titel = 'HI-density (cm^{-3})';
    maxx = max(max(Hrhohi));
    meann = mean(mean(Hrhohi));
    if (maxx - meann) > (3*maxx/4)
        Hrhohi = log10(Hrhohi);
        titel = 'HI-density (log_{10}(cm^{-3}))';
    end
    
    subplot(2,2,2), pcolor(x,x,Hrhohi)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')

    
    %----Solar Wind Velocity----

    %Plot speed
    Huswx = splitVector(dataH(:,12), span);
    Huswy = splitVector(dataH(:,13), span);
    Huswz = splitVector(dataH(:,14), span);
    speed = (Huswx.^2+Huswy.^2+Huswz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('SW-velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')
    
    %Scale and plot velocity arrows
    Huswx = scaleVec(Huswx, span, spanS);
    Huswz = scaleVec(Huswz, span, spanS);
    quiver(x,x,Huswx,Huswz,scale,'k')
    axis square
    hold on

    
    %----Heavy Ion Velocity----

    %Plot speed
    Huhix = splitVector(dataH(:,16), span);    
    Huhiy = splitVector(dataH(:,17), span);
    Huhiz = splitVector(dataH(:,18), span);
    speed = (Huhix.^2+Huhiy.^2+Huhiz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,4), pcolor(x,x, speed)
    shading interp
    axis square
    hold on
    colorbar
    title('HI-velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')

    %Scale and plot velocity arrows
    Huhix = scaleVec(Huhix, span, spanS);
    Huhiz = scaleVec(Huhiz, span, spanS);
    quiver(x,x,Huhix,Huhiz,scale,'k')
    axis square
    hold on
    

    Hybrid2 = figure('name','Hybrid2');

    %----Magnetic Field----

    %Plot field strength
    HBx = splitVector(dataH(:,5), span);
    HBy = splitVector(dataH(:,6), span);
    HBz = splitVector(dataH(:,7), span);
    speed = 10^9.*((HBx.^2+HBy.^2+HBz.^2).^(1/2));
    subplot(2,2,1), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('B-field (nT)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')

    %Plot field direction
    HBx = scaleVec(HBx, span, spanS);
    HBz = scaleVec(HBz, span, spanS);
    quiver(x,x,HBx,HBz,scale,'k')
    axis square
    hold on


    %----Electric field----

    %Plot field strength
    HEx = splitVector(dataH(:,8), span);
    HEy = splitVector(dataH(:,9), span);
    HEz = splitVector(dataH(:,10), span);
    speed = (HEx.^2+HEy.^2+HEz.^2).^(1/2);
    subplot(2,2,2), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('E-field (V/m)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')

    %Plot field direction
    HEx = scaleVec(HEx, span, spanS);
    HEz = scaleVec(HEz, span, spanS);
    quiver(x,x,HEx,HEz,scale,'k')
    axis square
    hold on


    %----Total density----
    
    HrhoTot = dataH(:,11) + dataH(:,15);
    HrhoTot = splitVector(HrhoTot, span);
    
    titel = 'Total density (cm^{-3})';
    maxx = max(max(HrhoTot));
    meann = mean(mean(HrhoTot));
    if (maxx - meann) > (3*maxx/4)
        HrhoTot = log10(HrhoTot);
        titel = 'Total density (log_{10}(cm^{-3}))';
    end
    
    subplot(2,2,3), pcolor(x,x,HrhoTot.*10^-6)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')
    
    
    %----Total Velocity----

    %Plot speed
    HuTotx = (dataH(:,11).*dataH(:,12) + dataH(:,15).*dataH(:,16))./(dataH(:,11) + dataH(:,15));
    HuTotx = splitVector(HuTotx, span);
    HuToty = (dataH(:,11).*dataH(:,13) + dataH(:,15).*dataH(:,17))./(dataH(:,11) + dataH(:,15));
    HuToty = splitVector(HuToty, span);
    HuTotz = (dataH(:,11).*dataH(:,14) + dataH(:,15).*dataH(:,18))./(dataH(:,11) + dataH(:,15));
    HuTotz = splitVector(HuTotz, span);
    speed = (HuTotx.^2+HuToty.^2+HuTotz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,4), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Total velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')

    %Scale and plot velocity arrows
    HuTotx = scaleVec(HuTotx, span, spanS);
    HuTotz = scaleVec(HuTotz, span, spanS);
    quiver(x,x,HuTotx,HuTotz,scale,'k')
    axis square
    hold on
        
end

if max(Hz) == 0
    
    %----Solar Wind Density----
        
    Hrhosw = splitVector(dataH(:,11), span).*10^-6;    
    subplot(2,2,1), pcolor(x,x,Hrhosw)
    shading interp
    axis square
    hold on
    colorbar
    title('SW-density (cm^{-3})')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')
    
    
    %----Heavy Ion Density----
    
    Hrhohi = splitVector(dataH(:,15), span).*10^-6;
    
    titel = 'HI-density (cm^{-3})';
    maxx = max(max(Hrhohi));
    meann = mean(mean(Hrhohi));
    if (maxx - meann) > (3*maxx/4)
        Hrhohi = log10(Hrhohi);
        titel = 'HI-density (log_{10}(cm^{-3}))';
    end
    
    subplot(2,2,2), pcolor(x,x,Hrhohi)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')

    
    %----Solar Wind Velocity----

    %Plot speed
    Huswx = splitVector(dataH(:,12), span);
    Huswy = splitVector(dataH(:,13), span);
    Huswz = splitVector(dataH(:,14), span);
    speed = (Huswx.^2+Huswy.^2+Huswz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('SW-velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')
    
    %Scale and plot velocity arrows
    Huswx = scaleVec(Huswx, span, spanS);
    Huswy = scaleVec(Huswy, span, spanS);
    quiver(x,x,Huswx,Huswy,scale,'k')
    axis square
    hold on

    
    %----Heavy Ion Velocity----

    %Plot speed
    Huhix = splitVector(dataH(:,16), span);    
    Huhiy = splitVector(dataH(:,17), span);
    Huhiz = splitVector(dataH(:,18), span);
    speed = (Huhix.^2+Huhiy.^2+Huhiz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,4), pcolor(x,x, speed)
    shading interp
    axis square
    hold on
    colorbar
    title('HI-velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')

    %Scale and plot velocity arrows
    Huhix = scaleVec(Huhix, span, spanS);
    Huhiy = scaleVec(Huhiy, span, spanS);
    quiver(x,x,Huhix,Huhiy,scale,'k')
    axis square
    hold on
    

    Hybrid2 = figure('name','Hybrid2');

    %----Magnetic Field----

    %Plot field strength
    HBx = splitVector(dataH(:,5), span);
    HBy = splitVector(dataH(:,6), span);
    HBz = splitVector(dataH(:,7), span);
    speed = 10^9.*((HBx.^2+HBy.^2+HBz.^2).^(1/2));
    subplot(2,2,1), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('B-field (nT)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')

    %Plot field direction
    HBx = scaleVec(HBx, span, spanS);
    HBy = scaleVec(HBy, span, spanS);
    quiver(x,x,HBx,HBy,2,'k')
    axis square
    hold on


    %----Electric field----

    %Plot field strength
    HEx = splitVector(dataH(:,8), span);
    HEy = splitVector(dataH(:,9), span);
    HEz = splitVector(dataH(:,10), span);
    speed = (HEx.^2+HEy.^2+HEz.^2).^(1/2);
    subplot(2,2,2), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('E-field (V/m)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')

    %Plot field direction
    HEx = scaleVec(HEx, span, spanS);
    HEy = scaleVec(HEy, span, spanS);
    quiver(x,x,HEx,HEy,2,scale,'k')
    axis square
    hold on


    %----Total density----
    
    HrhoTot = dataH(:,11) + dataH(:,15);
    HrhoTot = splitVector(HrhoTot, span);
    
    titel = 'Total density (cm^{-3})';
    maxx = max(max(HrhoTot));
    meann = mean(mean(HrhoTot));
    if (maxx - meann) > (3*maxx/4)
        HrhoTot = log10(HrhoTot);
        titel = 'Total density (log_{10}(cm^{-3}))';
    end
    
    subplot(2,2,3), pcolor(x,x,HrhoTot.*10^-6)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')
    
    
    %----Total Velocity----

    %Plot speed
    HuTotx = (dataH(:,11).*dataH(:,12) + dataH(:,15).*dataH(:,16))./(dataH(:,11) + dataH(:,15));
    HuTotx = splitVector(HuTotx, span);
    HuToty = (dataH(:,11).*dataH(:,13) + dataH(:,15).*dataH(:,17))./(dataH(:,11) + dataH(:,15));
    HuToty = splitVector(HuToty, span);
    HuTotz = (dataH(:,11).*dataH(:,14) + dataH(:,15).*dataH(:,18))./(dataH(:,11) + dataH(:,15));
    HuTotz = splitVector(HuTotz, span);
    speed = (HuTotx.^2+HuToty.^2+HuTotz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(2,2,4), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Total velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')

    %Scale and plot velocity arrows
    HuTotx = scaleVec(HuTotx, span, spanS);
    HuToty = scaleVec(HuToty, span, spanS);
    quiver(x,x,HuTotx,HuToty,scale,'k')
    axis square
    hold on
    
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: Hybrid}','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(Hybrid1,'Hyb1','jpg')
    saveas(Hybrid2,'Hyb2','jpg')  
    !lpr -P hp4600 1st.jpg
    !lpr -P hp4600 Hyb1.jpg
    !lpr -P hp4600 Hyb2.jpg
    !rm 1st.jpg
    !rm Hyb1.jpg
    !rm Hyb2.jpg
end