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

function plotMHDPlane(dataM,heli)

LP_Constants;

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

scale = span / (2 * spanS);

Mx = dataM(:,1).*10^-3;
My = dataM(:,2).*10^-3;
Mz = dataM(:,3).*10^-3;

coo = max([max(Mx) max(My) max(Mz)]);

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

MHD1 = figure('name','MHD1','position',[1 1 815 701]);
dist = heli;

if max(Mx) == 0
     
    %----Density----
    Mn = splitVector(dataM(:,5), span).*1e-6;
    
    titel = 'Density (cm^{-3})';
    maxx = max(max(Mn));
    meann = mean(mean(Mn));
    if (maxx - meann) > (3*maxx/4)
        Mn = log10(Mn);
        titel = 'Density (log_{10}(cm^{-3}))';
    end
    
    subplot(3,2,1), pcolor(x,x,Mn)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
    
    %----Particle Velocity----

    %Plot field strength
    MVx = splitVector(dataM(:,6), span);
    MVy = splitVector(dataM(:,7), span);
    MVz = splitVector(dataM(:,8), span);
    speed = (MVx.^2+MVy.^2+MVz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(3,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Particle velocity (km/s)')
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    

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

    %----Magnetic Field----

    %Plot field strength
    MBx = splitVector(dataM(:,9), span);
    MBy = splitVector(dataM(:,10), span);
    MBz = splitVector(dataM(:,11), span);
    speed = 10^9.*((MBx.^2+MBy.^2+MBz.^2).^(1/2));
    subplot(3,2,4), 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
    MBy = scaleVec(MBy, span, spanS);
    MBz = scaleVec(MBz, span, spanS);
    quiver(x,x,MBy,MBz,scale,'k')
    axis square
    hold on



    %----Temperature----
    MT = dataM(:,12)./(dataM(:,5).*CO.Kb1).*CO.Kb2;
    MT = splitVector(MT, span);
    
    titel = 'Temperature (eV)';
    maxx = max(max(MT));
    meann = mean(mean(MT));
    if (maxx - meann) > (3*maxx/4)
        MT = log10(MT);
        titel = 'Temperature (log_{10}(eV))';
    end
    
    subplot(3,2,2), pcolor(x,x,MT)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('y-axis (km)')
    ylabel('z-axis (km)')
    
end

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

    %----Particle Velocity----

    %Plot field strength
    MVx = splitVector(dataM(:,6), span);
    MVy = splitVector(dataM(:,7), span);
    MVz = splitVector(dataM(:,8), span);
    speed = (MVx.^2+MVy.^2+MVz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(3,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Particle velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')
    

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

    %----Magnetic Field----

    %Plot field strength
    MBx = splitVector(dataM(:,9), span);
    MBy = splitVector(dataM(:,10), span);
    MBz = splitVector(dataM(:,11), span);
    speed = 10^9.*((MBx.^2+MBy.^2+MBz.^2).^(1/2));
    subplot(3,2,4), 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
    MBx = scaleVec(MBx, span, spanS);
    MBz = scaleVec(MBz, span, spanS);
    quiver(x,x,MBx,MBz,scale,'k')
    axis square
    hold on



    %----Temperature----
    MT = dataM(:,12)./(dataM(:,5).*CO.Kb1).*CO.Kb2;
    MT = splitVector(MT, span);
    
    titel = 'Temperature (eV)';
    maxx = max(max(MT));
    meann = mean(mean(MT));
    if (maxx - meann) > (3*maxx/4)
        MT = log10(MT);
        titel = 'Temperature (log_{10}(eV))';
    end
    
    subplot(3,2,2), pcolor(x,x,MT)
    shading interp
    axis square
    hold on
    colorbar
    title(titel)
    xlabel('x-axis (km)')
    ylabel('z-axis (km)')
    
end

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

    %----Particle Velocity----

    %Plot field strength
    MVx = splitVector(dataM(:,6), span);
    MVy = splitVector(dataM(:,7), span);
    MVz = splitVector(dataM(:,8), span);
    speed = (MVx.^2+MVy.^2+MVz.^2).^(1/2);
    speed = speed.*10^-3;
    subplot(3,2,3), pcolor(x,x,speed)
    shading interp
    axis square
    hold on
    colorbar
    title('Particle velocity (km/s)')
    xlabel('x-axis (km)')
    ylabel('y-axis (km)')
    

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

    %----Magnetic Field----

    %Plot field strength
    MBx = splitVector(dataM(:,9), span);
    MBy = splitVector(dataM(:,10), span);
    MBz = splitVector(dataM(:,11), span);
    speed = 10^9.*((MBx.^2+MBy.^2+MBz.^2).^(1/2));
    subplot(3,2,4), 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
    MBx = scaleVec(MBx, span, spanS);
    MBy = scaleVec(MBy, span, spanS);
    quiver(x,x,MBx,MBy,scale,'k')
    axis square
    hold on



    %----Temperature----
    MT = dataM(:,12)./(dataM(:,5).*CO.Kb1).*CO.Kb2;
    MT = splitVector(MT, span);
        
    titel = 'Temperature (eV)';
    maxx = max(max(MT));
    meann = mean(mean(MT));
    if (maxx - meann) > (3*maxx/4)
        MT = log10(MT);
        titel = 'Temperature (log_{10}(eV))';
    end
    
    subplot(3,2,2), pcolor(x,x,MT)
    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'
    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
    dist = input('What distance is modelled? (1.3 AU/2.0 AU/2.7 AU/3.0 AU): ','s');
    subplot(3,2,5), text(0.1,0.9,'\bf{MHD}','fontsize',16)
    subplot(3,2,5), text(0.1,0.5,dist,'fontsize',16)    
    subplot(3,2,5), text(0.1,0.1,out,'fontsize',16)
    axis off
    saveas(MHD,'MHD','jpg')
    !lpr -P hp4600 MHD.jpg
    !rm MHD.jpg
end