%% Load geometry information
clear;
fromfile = true; % set true to load from file
fname = 'ex1_data.mat'; 
if fromfile && ~isempty(which(fname))
  load('ex1_data.mat'); % if file exists, load it
else
  % Loading geometry data structure |G| of TCV from |fge:|
  % Vessel model with 256 filaments.
  L = fge('tcv',61400,[],'selu','v','izgrid',1);
  L.G = meqg(L.G,L.P,'Brxu','Bzxu','Brxa','Bzxa'); 
  G=L.G;
  % uncomment this to save file
  % save(fname,'G');
end

% The G structures now contains the full filament model of the vessel, since above we used `selu='v'`.
% Note that G.Tvu is the identity matrix and G.Muu is not diagonal.
% We'll copy these variables to Mvv, Rv etc for clarity
G.Mvv = G.Muu; G.Mav = G.Mau; G.Rv = G.Ru; G.Mxv=G.Mxu; G.Brxv=G.Brxu; G.Bzxv=G.Bzxu; G.Rvv = diag(G.Rv); % since selu='v'
fprintf('Mvv diagonal? '); if isdiag(G.Mvv); fprintf('yes\n'); else, fprintf('no\n'); end
fprintf('nv = %d\n',G.nv)
% copy some more variables for brevity
na = G.na; nv = G.nv;

%% Eigenmode decomposition
[V,D] = eig(G.Mvv\G.Rvv); % eigenvector and eigenvalues

% esort(P)  sorts the complex eigenvalues in the vector P in
% descending order by real part.  The unstable eigenvalues
% (in the continuous-time sense) will appear first.
[~,isort] = esort(diag(-D)); D = D(isort,isort); V = V(:,isort); % reorder eigenvalues
% keep only neig eigenvalues
neig = 30; ieig = 1:neig; 

% Define new Tvu and Tuv to diagonalize Muu and Ruu
Tvu = V; Tuv = V'; % this choice of Tuv maintains Mua=Mau'
Tuv = Tuv(ieig,:); Tvu = Tvu(:,ieig); % select only subset of eigenmodes
Ruu = Tuv*diag(G.Rv)*Tvu; Muu = Tuv*G.Mvv*Tvu; Mau = G.Mav*Tvu; Mua = Tuv*G.Mav';

% check if new Ruu, Muu are diagonal
% seems of diagonal values are nonzero, due to numerical round off, rounding to 10 digits shows diagonal matrices though
fprintf('Muu diagonal? '); if isdiag(round(Muu,10)); fprintf('yes\n'); else, fprintf('no\n'); end
fprintf('Ruu diagonal? '); if isdiag(round(Ruu,10)); fprintf('yes\n'); else, fprintf('no\n'); end



%% Some helpful MEQ fields for later on in the exercise

% The names of the TCV coils can be found in G.dima. Their order
% corresponds to their use elsewhere in MEQ (e.g. in G.Ra and G.Maa)
% The location the coils can be obtained from G.rc and G.zc

% The position of the current elements x in the plasma are collected in
% G.rx and G.zx 
% The mutual inductances between the current elements x and the active
% coils and the vessel are in the matrices G.Mxa and G.Mxu
% The radial and vertical magnetic fields are related to the currents
% through the matrices G.Brxa and G.Brxu for the radial field and G.Bzxa
% and G.Bzxu for the vertical field.

% These long 1D array storing data corresponding to R and Z locations can
% be converted to a 2D structure as follows:
Mxcoil1_2D = reshape(G.Mxa(:,1),G.nz,G.nr); % mutual inductance between the first active coil and all current elements in a 2D matrix structure
% The same can be done for other quantities of course






