% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Octave/MATLAB script computes the radial magnetic field Br (denoted as "By" here) from a 2D map of the axial field Bz. 
% Rename this with a ".m" extension, so that Octave/MATLAB recognize it.
% 
% Enjoy!
% Stefano Boccelli, 2020. 
% boccelliengineering.altervista.org
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear
clc

try
  page_screen_output(0);
end

% ++++++++++++++ WRITE HERE DATA ++++++++++++++

z = [41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67]; % [mm] Axial positions from the support
hh = z - z(1) + 3; % [mm] scale them to distance from the electrode surface

B_mV = [-40.15 -51.86 -57.35 -65.82 -74.02 -81.97 -88.69 -93.60 -97.27 -100.71 -103.56 -105.78 -108.12 -107.53 -102.97 -97.36 -92.49 -85.32 -79.89 -72.41 -62.49 -54.78 -44.15 -34.32 -24.43 -14.39 -4.25 6.72 14.91 24.67 38.79 49.34 66.13 77.71 92.98 104.98 109.61 122.12 135.08 148.58 153.59 160.82 172.07 171.79 176.01 188.41 189.21 183.61 184.36 182.50;
-42.59 -42.41 -45.87 -50.42 -54.65 -58.99 -62.72 -65.07 -66.64 -68.30 -69.58 -70.52 -71.71 -71.24 -68.41 -64.95 -61.97 -57.62 -54.86 -49.55 -44.14 -39.10 -32.60 -27.02 -20.37 -14.18 -7.44 0.15 5.64 11.84 21.27 27.81 36.87 45.55 54.55 62.63 64.71 72.69 80.50 88.82 92.02 95.96 103.47 103.39 105.75 111.95 113.27 109.39 109.77 105.97;
-34.88 -34.61 -36.38 -39.16 -41.69 -44.25 -46.29 -47.76 -48.49 -49.33 -49.98 -50.54 -51.11 -50.63 -48.90 -46.62 -44.64 -41.88 -39.86 -37.25 -32.23 -29.96 -25.72 -22.27 -18.44 -13.83 -9.01 -4.37 -0.91 3.32 9.35 13.72 19.25 25.28 31.28 36.13 37.60 42.64 47.54 52.33 54.90 57.28 61.88 61.74 62.93 67.06 67.56 65.05 65.02 63.44;
-24.51 -27.63 -29.02 -30.74 -32.27 -33.86 -35.18 -36.09 -36.60 -37.03 -37.41 -37.81 -38.18 -37.91 -36.73 -35.21 -34.15 -32.45 -31.19 -29.35 -27.21 -24.89 -22.26 -19.99 -16.85 -14.16 -11.61 -8.35 -6.03 -3.23 0.56 3.29 7.84 10.99 14.92 17.86 18.96 22.09 25.64 28.44 30.13 31.49 34.40 34.26 35.19 38.19 37.98 36.20 36.09 35.00;
-22.37 -21.88 -22.70 -23.91 -24.93 -26.02 -26.79 -27.36 -27.71 -27.97 -28.24 -28.48 -28.83 -28.68 -27.94 -27.14 -26.36 -25.33 -24.62 -23.55 -22.21 -20.69 -19.20 -17.80 -16.38 -14.49 -12.71 -10.73 -9.39 -7.77 -5.29 -3.71 -1.03 1.06 3.44 5.38 5.96 7.93 9.80 12.01 12.75 13.70 15.46 15.29 15.72 17.41 17.50 16.35 16.29 15.79;
-18.19 -17.96 -18.57 -19.33 -20.03 -20.78 -21.34 -21.73 -21.94 -22.11 -22.34 -22.59 -22.76 -22.74 -22.27 -21.66 -21.26 -20.63 -20.16 -19.53 -18.62 -17.79 -16.80 -15.96 -14.87 -13.79 -12.68 -11.44 -10.57 -9.52 -7.96 -6.84 -5.23 -3.94 -2.45 -1.20 -0.78 0.46 1.66 3.04 3.48 4.04 5.16 5.11 5.33 6.31 6.34 5.57 5.54 5.19;
-14.09 -14.76 -15.13 -15.72 -16.23 -16.73 -17.14 -17.46 -17.61 -17.72 -17.87 -18.04 -18.24 -18.25 -17.91 -17.63 -17.34 -16.99 -16.71 -16.36 -15.83 -15.28 -14.66 -14.19 -13.56 -12.80 -12.17 -11.31 -10.89 -10.16 -9.29 -8.66 -7.67 -6.84 -5.94 -5.18 -4.90 -4.16 -3.29 -2.63 -2.36 -2.02 -1.38 -1.43 -1.29 -0.64 -0.71 -1.18 -1.24 -1.45;
-12.40 -12.52 -12.86 -13.24 -13.61 -14.02 -14.33 -14.56 -14.65 -14.75 -14.90 -15.09 -15.24 -15.24 -15.01 -14.82 -14.64 -14.41 -14.28 -14.08 -13.75 -13.39 -13.01 -12.72 -12.29 -11.87 -11.39 -10.96 -10.64 -10.16 -9.59 -9.21 -8.66 -8.10 -7.53 -7.02 -6.79 -6.37 -5.92 -5.48 -5.25 -4.99 -4.64 -4.66 -4.59 -4.27 -4.33 -4.61 -4.64 -4.86;
-10.89 -10.64 -10.87 -11.23 -11.49 -11.81 -12.02 -12.23 -12.27 -12.39 -12.51 -12.64 -12.81 -12.88 -12.76 -12.63 -12.51 -12.38 -12.34 -12.24 -12.03 -11.81 -11.56 -11.42 -11.26 -10.92 -10.63 -10.39 -10.21 -9.89 -9.54 -9.40 -9.02 -8.66 -8.35 -8.03 -7.88 -7.63 -7.38 -7.12 -6.98 -6.86 -6.63 -6.63 -6.62 -6.48 -6.47 -6.62 -6.71 -6.86;
-8.79 -8.79 -8.99 -9.24 -9.49 -9.73 -9.87 -9.99 -10.12 -10.22 -10.34 -10.48 -10.62 -10.65 -10.62 -10.53 -10.50 -10.50 -10.50 -10.49 -10.39 -10.27 -10.22 -10.13 -10.02 -9.88 -9.76 -9.64 -9.53 -9.39 -9.23 -9.13 -8.94 -8.76 -8.63 -8.42 -8.36 -8.18 -8.11 -7.92 -7.87 -7.79 -7.73 -7.68 -7.66 -7.63 -7.65 -7.74 -7.76 -7.87;
-7.24 -7.62 -7.76 -7.99 -8.14 -8.36 -8.49 -8.62 -8.75 -8.85 -8.94 -9.01 -9.21 -9.27 -9.25 -9.25 -9.25 -9.24 -9.25 -9.26 -9.26 -9.21 -9.14 -9.16 -9.13 -9.01 -8.99 -8.89 -8.87 -8.76 -8.66 -8.63 -8.56 -8.49 -8.37 -8.27 -8.20 -8.13 -8.09 -8.00 -7.97 -7.89 -7.88 -7.88 -7.87 -7.88 -7.88 -7.89 -7.96 -8.00 ;
-6.49 -6.63 -6.75 -6.97 -7.12 -7.27 -7.39 -7.50 -7.61 -7.68 -7.79 -7.89 -8.01 -8.12 -8.11 -8.11 -8.13 -8.16 -8.20 -8.25 -8.26 -8.24 -8.24 -8.26 -8.25 -8.21 -8.14 -8.13 -8.14 -8.10 -8.04 -8.10 -8.01 -7.99 -7.96 -7.88 -7.79 -7.76 -7.76 -7.75 -7.74 -7.71 -7.68 -7.66 -7.65 -7.71 -7.72 -7.75 -7.75 -7.79;
-5.64 -5.77 -5.89 -6.02 -6.18 -6.35 -6.44 -6.51 -6.63 -6.65 -6.78 -6.87 -6.99 -7.06 -7.04 -7.11 -7.11 -7.16 -7.24 -7.25 -7.31 -7.27 -7.33 -7.38 -7.38 -7.38 -7.37 -7.37 -7.38 -7.37 -7.37 -7.38 -7.38 -7.38 -7.36 -7.33 -7.25 -7.24 -7.26 -7.26 -7.25 -7.25 -7.25 -7.24 -7.25 -7.27 -7.29 -7.27 -7.29 -7.35;
-4.90 -5.02 -5.12 -5.24 -5.38 -5.49 -5.63 -5.69 -5.74 -5.83 -5.90 -5.99 -6.10 -6.19 -6.17 -6.25 -6.25 -6.32 -6.37 -6.44 -6.46 -6.49 -6.51 -6.58 -6.62 -6.61 -6.63 -6.61 -6.63 -6.63 -6.63 -6.68 -6.73 -6.70 -6.74 -6.72 -6.63 -6.66 -6.71 -6.74 -6.68 -6.71 -6.73 -6.67 -6.71 -6.75 -6.75 -6.75 -6.75 -6.77;
];

% ++++++++++++++ START POSTPROCESSING +++++++++++++++++++

B_mV = B_mV(:,2:end); % Remove first point (it's always noisy)

rr = linspace(23, 0, size(B_mV,2)) - 1.5;

% Convert to Gauss
B_mT = B_mV*1800/2100;
Gs   = B_mT*10;

% ++++++++++++ Create a matrix for contourplot ++++++++++++++++
GsMAT = Gs';
[hhMAT, rrMAT] = meshgrid(hh, rr);

% Remove the last four elements inside rr, that go over the axis.
% (only use data from the axis (r = 0), outwards)
hhMAT = hhMAT(1:end-4,:);
rrMAT = rrMAT(1:end-4,:);
Bz    = Gs(:,1:end-4);

rr = rr(1:end-4);
hh = hh;

% Compute z derivative for every radial position (rr)
for ii = 1:numel(rr) % Loop on height coordinate (axis z)
  dz_Bz(:, ii) = diff(Bz(:,ii))./diff(hh)';
end

% Now, compute By
By = zeros(size(Bz));

dr = rr(2) - rr(1);

for ii = 1:numel(hh)-1

  By_now = 0 + 1./(rr(end:-1:1)).*cumsum(rr(end:-1:1).*dz_Bz(ii,:))*dr;
  By(ii+1, :) = By_now(end:-1:1);

end

By = - By;

% Smooth a bit (kind of Laplacian)
for ii = 1:5

  for kk = 1:size(By,1)
    By(kk,2:end-1) = (By(kk,3:end) +2*By(kk, 2:end-1) + By(kk, 1:end-2))/4;
  end

end

% Plot contour of Bz and By fields
figure
subplot(1,2,1)
contourf(rrMAT, hhMAT, Bz', 20)
title('Axial field Bz [Gauss]')
xlabel('Radial position "r" [mm]')
ylabel('Axial position "z" [mm]')
colorbar
axis equal

subplot(1,2,2)
contourf(rrMAT, hhMAT, By', 20)
colorbar
axis equal
title('Radial field Br [Gauss]')
xlabel('Radial position "r" [mm]')
ylabel('Axial position "z" [mm]')

% ++++++ Plot field lines ++++++
% QUIVER
figure
quiver(rrMAT, hhMAT, By', Bz', 'linewidth',2)
axis equal

% NOW PLOT SOME FIELD LINES
dx = 0.1; % Step for tracing field lines

% r_start = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]; % Seed starting points
% r_start = [10:0.5:20];

r_start = [min(rr),1,2,3,4,5,6,7,8]; % Seed starting points

for kkk = 1:numel(r_start)

  fprintf('Computing field line %d of %d\n', kkk, numel(r_start));
 
  r(1) = r_start(kkk);
  z(1) = min(hh);
  
  for ii = 2:400 % Loop for a while
  
    r_now = r(ii-1);
    z_now = z(ii-1);
  
    Br_now = interp2(rrMAT', hhMAT', By, r_now, z_now);
    Bz_now = interp2(rrMAT', hhMAT', Bz, r_now, z_now);
  
    Babs = sqrt(Br_now^2 + Bz_now^2);
  
    r(ii) = r(ii-1) + Br_now/Babs*dx;
    z(ii) = z(ii-1) + Bz_now/Babs*dx;

    if (z(ii) < 2)
      break
    end

  end
  
  hold on
  plot(r, z, 'r', 'linewidth',2)

end