clear
clc
close all

% Read in the kinetic energy and position of ions and the magnetic field
% strength at their location
[KE0 B0 x0 y0 z0]=reads('Run 5 Magnetic0e-4.csv');
[KE1 B1 x1 y1 z1]=reads('Run 5 Magnetic1e-4.csv');
[KE2 B2 x2 y2 z2]=reads('Run 5 Magnetic2e-4.csv');
[KE3 B3 x3 y3 z3]=reads('Run 5 Magnetic3e-4.csv');
[KE4 B4 x4 y4 z4]=reads('Run 5 Magnetic4e-4.csv');
[KE5 B5 x5 y5 z5]=reads('Run 5 Magnetic5e-4.csv');
[KE6 B6 x6 y6 z6]=reads('Run 5 Magnetic6e-4.csv');
[KE7 B7 x7 y7 z7]=reads('Run 5 Magnetic7e-4.csv');
[KE8 B8 x8 y8 z8]=reads('Run 5 Magnetic8e-4.csv');
[KE9 B9 x9 y9 z9]=reads('Run 5 Magnetic9e-4.csv');
%[KE10 B10 x10 y10 z10]=reads('Run 5 Magnetic10e-4.csv');

% Add all variables into a long vector to simulate a continuous flow
% instead of one burst.
Total_KE = [KE0',KE0',KE0',KE0',KE0',KE0',KE0',KE0',KE0',KE0',KE1',KE1',KE1',KE1',KE1',KE1',KE1',KE1',KE1',KE2',KE2',KE2',KE2',KE2',KE2',KE2',KE2',KE3',KE3',KE3',KE3',KE3',KE3',KE3',KE4',KE4',KE4',KE4',KE4',KE4',KE5',KE5',KE5',KE5',KE5',KE6',KE6',KE6',KE6',KE7',KE7',KE7',KE8',KE8',KE9'];
Total_KE = Total_KE';
Total_B = [B0',B0',B0',B0',B0',B0',B0',B0',B0',B0',B1',B1',B1',B1',B1',B1',B1',B1',B1',B2',B2',B2',B2',B2',B2',B2',B2',B3',B3',B3',B3',B3',B3',B3',B4',B4',B4',B4',B4',B4',B5',B5',B5',B5',B5',B6',B6',B6',B6',B7',B7',B7',B8',B8',B9'];
Total_B = Total_B';
Total_X = [x0',x0',x0',x0',x0',x0',x0',x0',x0',x0',x1',x1',x1',x1',x1',x1',x1',x1',x1',x2',x2',x2',x2',x2',x2',x2',x2',x3',x3',x3',x3',x3',x3',x3',x4',x4',x4',x4',x4',x4',x5',x5',x5',x5',x5',x6',x6',x6',x6',x7',x7',x7',x8',x8',x9'];
Total_X = Total_X';
Total_Y = [y0',y0',y0',y0',y0',y0',y0',y0',y0',y0',y1',y1',y1',y1',y1',y1',y1',y1',y1',y2',y2',y2',y2',y2',y2',y2',y2',y3',y3',y3',y3',y3',y3',y3',y4',y4',y4',y4',y4',y4',y5',y5',y5',y5',y5',y6',y6',y6',y6',y7',y7',y7',y8',y8',y9'];
Total_Y = Total_Y';
Total_Z = [z0',z0',z0',z0',z0',z0',z0',z0',z0',z0',z1',z1',z1',z1',z1',z1',z1',z1',z1',z2',z2',z2',z2',z2',z2',z2',z2',z3',z3',z3',z3',z3',z3',z3',z4',z4',z4',z4',z4',z4',z5',z5',z5',z5',z5',z6',z6',z6',z6',z7',z7',z7',z8',z8',z9'];
Total_Z = Total_Z';

% Declare the needed variables.
counter = 0;
total = size(Total_X,1);
Zspacing = 0.5; % spacing for the Z direction, m
Rspacing = 0.5; % spacing for the Radial direction, m
halfZ = Zspacing/2;
halfR = Rspacing/2;
R_total = 60/Rspacing+2;
[Z,R] = meshgrid(-30.5:Zspacing:30,-30.5:Rspacing:30);
B_avg = zeros(R_total,R_total);
V_avg = zeros(R_total,R_total);
V_Alfven = zeros(R_total,R_total);
Beta = zeros(R_total,R_total);
Density = zeros(R_total,R_total);
Density_kg = zeros(R_total,R_total);
Volume =(Rspacing*Zspacing*0.1);
Z2 = linspace(-30,30,R_total);
scale = 2.2612e21/size(Total_X,1);

% Step 1: Cut a slice out of the 30 m sphere.
for i=1:total
    if Total_X(i)<=0.05 && Total_X(i)>=-0.05 && Total_Y(i)<29.9% slice is -0.05m < x < 0.05m
        counter = counter+1;
        SliceR(counter)  = Total_Y(i);
        SliceZ(counter)  = Total_Z(i);
        SliceB(counter)  = Total_B(i);
        SliceKE(counter) = Total_KE(i);
    end
end

SliceV = sqrt(2*SliceKE/6.633e-26); % Determine velocity magnitude from KE.

% Step 2: Bin the sliced ions into Z x R bins
for i=1:R_total
    Rindex = -31 + i*Rspacing; 
    for j=1:R_total
       Zindex = -31 + j*Zspacing;
       number = 0;
       Bsum = 0;
       Vsum = 0;
        for k=1:counter
            if SliceR(k)<=(Rindex+halfR) && SliceZ(k)<=(Zindex+halfZ) && SliceR(k)>(Rindex-halfR) && SliceZ(k)>(Zindex-halfZ)
                number = number + 1;
                Bsum = Bsum + SliceB(k);
                Vsum = Vsum + SliceV(k);
            end
        end
       Density(i,j)= number*scale/Volume; % # of ions/m^3 x a scaling factor
       Density_kg(i,j)=Density(i,j)*(150e-6/2.2612e21);% kg/m^3
       
       if number > 0
            B_avg(i,j)= Bsum/number; % Average B in this square
            V_avg(i,j)= Vsum/number; % Average V in this square
       end
       Total(i,j)=number;
    end
end

Count = sum(sum(Total,1)); % check to make sure total number of ions was counted.

for l=1:R_total
    for m=1:R_total
        if Density_kg(l,m)>0
            V_Alfven(l,m) = (B_avg(l,m))/sqrt((1.25664e-6)*Density_kg(l,m));
            Beta(l,m) = (V_avg(l,m)/V_Alfven(l,m))^2;
        end
    end
end

for l=1:R_total
    for m=1:R_total
        if Beta(l,m)>0
          Log_Beta(l,m)=log10(Beta(l,m));
        else Log_Beta(l,m)=0;
        end
    end
end

% Interpolate to get a smoother graph.
INT_Log_Beta = INTERP2(Log_Beta,2);
INT_Density = INTERP2(Density,2);
INT_VAlfven = INTERP2(V_Alfven,2);
[Zi,Ri] = meshgrid(-30.5:Zspacing/4:30,-30.5:Rspacing/4:30);

figure
contourf(Ri,Zi,INT_Density',50,'LineColor','none')
axis([2.8 5 -0.7 0])
xlabel('Axial distance, m')
ylabel('Radius, m')
colorbar
%caxis([0 5])

figure
contourf(Ri,Zi,INT_Log_Beta',50,'LineColor','none')
axis([-1 29 -30 30])
%axis([.5 5 -0.7 0])
xlabel('Axial distance, m')
ylabel('Radius, m')
caxis([-2 2])
colorbar

figure
h = bar3(Z2,Log_Beta);
colormap jet
axis([0 R_total -30 30 -Inf Inf])
ylabel('Axial distance, m')
xlabel('Radius, m')
title('Beta')
set(gca,'XTickLabel',{'-30';'-20';'-10';'0';'10';'20';'30'}) 
colorbar
caxis([-2 5])
for i = 1:length(h)
    zdata = get(h(i),'ZData');
    set(h(i),'CData',zdata) % Add back edge color removed by interpolating shading
    set(h,'EdgeColor','k') 
end


