%% Exercise 3.1b ------------------------------------------------------------
% Solution to exercise 3.1b
%
% Originally prepared by Artur Palha, 2015/20/02
% Updated and adapted by Pedro Molina, 04/02/2023
%--------------------------------------------------------------------------

%% Load TCV geoemtry data -------------------------------------------------
clear;
load('ex3.mat'); % loads L variable

%% Unit current magnetic field contribution -------------------------------
% for example: how to plot the contribution from a given coil to 
% a total magnetic field (Br,Bz). Apply a unit current (1A) on each coil for
% instance and see what field you get. 

figure('position', [0, 0, 1200, 800])        % the Br plot
Br_figure = gcf;
subplot(4,4,1) 
figure('position', [0, 0, 1200, 800])        % the Bz plot
Bz_figure = gcf;
subplot(4,4,1)
figure('position', [0, 0, 1200, 800])        % the (Br,Bz) vector plot
B_figure = gcf;
subplot(4,4,1)
% prepare R/Z grids for plotting nicely
[meshR, meshZ]  = meshgrid(L.G.rx, L.G.zx);
% collated R/Z grid
colR            = reshape(meshR, size(G.Brxa(:,1),1), 1);
colZ            = reshape(meshZ, size(G.Brxa(:,1),1), 1); 

% plot effect of coils E-F: 1-16. See G.dima
for n_coil=1:16
    figure(Br_figure)
    subplot(4,4,n_coil)
    % make the plots

    % plot Br for coil number n_coil 
    pcolor( meshR, meshZ, reshape(G.Bzxa(:,n_coil),L.nzx,L.nrx))
    shading interp
    colorbar
    %caxis([-0.02,0.02])
    axis([0.5,1.3,-0.9,0.9])
    box on
    title('$B_{r}$','Interpreter','latex')
    xlabel('$r$','Interpreter','latex')
    ylabel('$z$','Interpreter','latex')

    % plot Bz for coil number n_coil
    figure(Bz_figure)    
    subplot(4,4,n_coil) 
    pcolor( meshR, meshZ, reshape(G.Brxa(:,n_coil),L.nzx,L.nrx))
    shading interp
    colorbar
    %caxis([-0.02,0.02])
    axis([0.5,1.3,-0.9,0.9])
    box on
    title('$B_{z}$','Interpreter','latex')
    xlabel('$r$','Interpreter','latex')
    ylabel('$z$','Interpreter','latex')

    % plot the vector plot of the magnetic field for coil number n_coil
    figure(B_figure)    
    subplot(4,4,n_coil) 
    quiver( colR(1:10:end),colZ(1:10:end),G.Brxa(1:10:end,n_coil),...
            G.Bzxa(1:10:end,n_coil), 2)
    axis([0.5,1.3,-0.9,0.9])
    box on
    title('$\vec{B}$','Interpreter','latex')
    xlabel('$r$','Interpreter','latex')
    ylabel('$z$','Interpreter','latex')
end


%% Arbitrary current total magnetic field ---------------------------------

% plot the total magnetic field for a combination of currents in the coils
Ia = rand(16,1); % just generate random currents between -1 and 1 A on the
                 % coils

% compute the total magnetic field (just add the contribution from each
% E and F coil 1:16 in Ia. See G.dima)
B_r_total = G.Brxa(:,1:end-3)*Ia; 
B_z_total = G.Bzxa(:,1:end-3)*Ia;

% make the plots
figure('position', [0, 0, 1200, 600])
subplot(1,3,1) % plot Br total
pcolor(meshR, meshZ, reshape(B_r_total,L.nzx,L.nrx))
shading interp 
colorbar
%caxis([-0.02,0.02])
axis([0.5,1.3,-0.9,0.9])
box on
title('$B_{r}$','Interpreter','latex')
xlabel('$r$','Interpreter','latex')
ylabel('$z$','Interpreter','latex')
subplot(1,3,2) % plot Bz total
pcolor(meshR, meshZ, reshape(B_z_total,L.nzx,L.nrx))
shading interp
colorbar
%caxis([-0.02,0.02])
axis([0.5,1.3,-0.9,0.9])
box on
title('$B_{z}$','Interpreter','latex')
xlabel('$r$','Interpreter','latex')
ylabel('$z$','Interpreter','latex')
subplot(1,3,3) % plot the vector plot of the total magnetic field
quiver(colR, colZ, B_r_total, B_z_total,2)
axis([0.5,1.3,-0.9,0.9])
box on
title('$\vec{B}$','Interpreter','latex')
xlabel('$r$','Interpreter','latex')
ylabel('$z$','Interpreter','latex')


%% Least-Squares fit ------------------------------------------------------

% the Least-Squares problem tries to solve a non-quadratic system of
% equations:
%   Ax = b
% by minimizing the error.
% This results in the following algebraic system of equations:
% A^{t}Ax = A^{t}b


% the required field components
Br_required = 0.0;
Bz_required = -0.1183;

% compute the number of equations, which is the same a the number of points
% where where we want the magnetic field to be the required one
% Use only E and F coils. Thus n_equations = length(G.dima)-3
n_equations = size(G.Brxa(:,1:end-3),1);

% construct b (the right hand side)
br = Br_required*ones(n_equations,1); % the r component
bz = Bz_required*ones(n_equations,1); % the z component
b = [br;bz]; % the two components together

% construct the A matrix
A = [G.Brxa(:,1:end-3);G.Bzxa(:,1:end-3)];

% solve the system
Ia = (A'*A)\(A'*b);

% plot the solution
% compute the total magnetic field (just add the contribution from each
% coil)
B_r_total = G.Brxa(:,1:end-3)*Ia; 
B_z_total = G.Bzxa(:,1:end-3)*Ia;

% make the plots
figure('position', [0, 0, 1200, 600])
subplot(1,3,1) % plot Br total
pcolor(meshR, meshZ,reshape(B_r_total,L.nzx,L.nrx))
shading interp
colorbar
caxis([-0.02,0.02])
axis([0.5,1.3,-0.9,0.9])
box on
title('$B_{r}-B_{r,{expected}}$','Interpreter','latex')
xlabel('$r$','Interpreter','latex')
ylabel('$z$','Interpreter','latex')
subplot(1,3,2) % plot Bz total
pcolor(meshR, meshZ, reshape(B_z_total,L.nzx,L.nrx)-Bz_required)
shading interp
colorbar
caxis([-0.02,0.02])
axis([0.5,1.3,-0.9,0.9])
box on
title('$B_{z}-B_{z,\mathrm{expected}}$','Interpreter','latex')
xlabel('$r$','Interpreter','latex')
ylabel('$z$','Interpreter','latex')
subplot(1,3,3) % plot the vector plot of the total magnetic field
quiver(colR, colZ, B_r_total, B_z_total,2)
axis([0.5,1.3,-0.9,0.9])
box on
title('$\vec{B}$','Interpreter','latex')
xlabel('$r$','Interpreter','latex')
ylabel('$z$','Interpreter','latex')

%if save_fig %save file externally
%    export_fig(['./figures/', 'B_Ia_LS.png'],'-png','-transparent','-r300','-zbuffer')
%end

% compute the residual
R2 = (sum((B_z_total - Bz_required).^2) + sum(B_r_total.^2));

fprintf('\nResidual of LS approximation: %.4f\n\n',R2);
fprintf('\nAverage error of LS approximation: %.4f\n\n',sqrt(R2/(2*n_equations)));


