function G = liugtcv_IDS(P, IDS_inputs)
%% The IDS standard has the following names
% Power supplies
% Coils are groups of windings
% Circuits are groups of coils attached to a power supply
% Circuits are typically all series connected
% Coils can be diveded into elements. This distinction is not present yet
% in TCV
% Power supplies + Coils = Circuits
% The IDS_inputs.pf_system(1).pfcircuits.connections(circuit, node, element pin attach
% flag) has the following structure
% size(connections) = [ncircuits, max( nnodes( circuit(:))), 2* nelements]
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% Active pf circits. 
%The active pf circuits considered for the reconstruction migh be a subset
%of the total active pf circuits available, and connect together active
%coils, which if in series connection has to share the same current

% Computation of the transfer matrix
% Twa = from active coils considered to windings Iw = G.Twa*Ia
% Twc = from coils to windings Iw = G.Twc * Ic. 
% Tca = from coils to active circuits considered Ic = Tca*Ia
% These matrix is used to group together windings that share the same current
% Matrix is composed by columns of 1/0

% -------------- I will assume that all the active circuits are series
% connected, there would be no reason why the should not, and how otherwise they would compute the current--------------------------%
% This assumption will be relaxed 
% ---------------------------------------------------

G.pf_active = IDS_inputs.pf_active;
G.pf_passive = IDS_inputs.pf_passive;
G.wall  = IDS_inputs.wall;

% In the G structure some information obtained by reading the connection matrix are added
G.pf_active.pfsupplies = struct();
G.pf_active.pfcircuits = struct();
G.pf_active.pfcoils = struct();

% Add extra info to the structure.
G.pf_active.pfsupplies.dim = numel(G.pf_active.supply);
G.pf_active.pfcircuits.dim = numel(G.pf_active.circuit);
G.pf_active.pfcoils.dim = numel(G.pf_active.coil);

% Get all the names in cell array
G.pf_active.pfsupplies.names = cell(G.pf_active.pfsupplies.dim,1);
for ii=1:G.pf_active.pfsupplies.dim
  G.pf_active.pfsupplies.names{ii} =  G.pf_active.supply{ii}.name;
end
G.pf_active.pfcircuits.names = cell(G.pf_active.pfcircuits.dim,1);
for ii=1:G.pf_active.pfcircuits.dim
  G.pf_active.pfcircuits.names{ii} =  G.pf_active.circuit{ii}.name;
end
G.pf_active.pfcoils.names = cell(G.pf_active.pfcoils.dim,1);
for ii=1:G.pf_active.pfcoils.dim
  G.pf_active.pfcoils.names{ii} =  G.pf_active.coil{ii}.name;
end

% Get the number of elements for each coil
G.pf_active.pfcoils.n_elements_4_each_coils = zeros(G.pf_active.pfcoils.dim,1);
for  ii=1:G.pf_active.pfcoils.dim
  G.pf_active.pfcoils.n_elements_4_each_coils(ii) = numel(G.pf_active.coil{ii}.element);
end

%-----------------------------------------------------------------
% ------TODO connection matrix could be different for every circuits
% ----------------------------------------------------------------
G.pf_active.pfcircuits.connections = G.pf_active.circuit{1}.connections;
% Interpreting connection matrix to extract relation between coils and active
% circuits
if ~isempty(G.pf_active.pfcircuits.connections)
[G.pf_active.pfsupplies,G.pf_active.pfcircuits, G.pf_active.pfcoils]= read_connection_matrix(G.pf_active.pfsupplies,G.pf_active.pfcircuits, G.pf_active.pfcoils);
end

% Construct the grouping matrix coils to windings Twc. This has to be construct dinamically beacause every coil can have a different decription.
G.rw = [];
G.zw = [];
G.Twc = [];
for ii = 1:G.pf_active.pfcoils.dim % loop on coil
  for jj = 1:G.pf_active.pfcoils.n_elements_4_each_coils(ii) % loop on element of each coil
    switch  G.pf_active.coil{ii}.element{jj}.geometry.geometry_type
     case 1 % outline
      % TODO
      error('outline not implemented')
     case 2 % rectangle
      pfrzdrdz = zeros(4,1);
      pfrzdrdz(1) = G.pf_active.coil{ii}.element{jj}.geometry.rectangle.r;
      pfrzdrdz(2) = G.pf_active.coil{ii}.element{jj}.geometry.rectangle.z;
      pfrzdrdz(3) = G.pf_active.coil{ii}.element{jj}.geometry.rectangle.width;
      pfrzdrdz(4) = G.pf_active.coil{ii}.element{jj}.geometry.rectangle.height;
      turns_with_sign = G.pf_active.coil{ii}.element{jj}.turns_with_sign;
      [rw_tmp, zw_tmp, Twc_tmp] = eval_trap(pfrzdrdz, turns_with_sign, 2); 
     case 3 % oblique
      % TODO
            error('outline not implemented')
     case 4 % arcs circle
      % TODO
            error('outline not implemented')
    end
    G.rw = [G.rw; rw_tmp];
    G.zw = [G.zw; zw_tmp];
    nw_tmp = numel(rw_tmp);
    G.Twc(end+1:end+nw_tmp,end+1) = 1;   
  end
end

% The number of nodes in a circuit are equal to the number of coils + power supply connected in case of series connection
% If the there is only one power supply attached, than nnodes = ncoils +1

%----------- LIUQE variables name convention ----------------- %
G.dima = P.active_pfcircuits_considered; % Names of the active pf circuits considered for the reconstruction%
G.na = numel(G.dima); % Number of active pf circuits used 

if ~isempty(G.pf_active.pfcircuits.connections) % If connection matrix available use a subset of coils specified in parameters
% Checks requested coils are available
check_active_coils_requested(G.na, G.pf_active.pfcircuits.dim, G.pf_active.pfcircuits.names, P.active_pfcircuits_considered)

G.Tca = zeros(G.pf_active.pfcoils.dim, G.na);
for ii=1:G.pf_active.pfcoils.dim
    tmp = strmatch(G.pf_active.pfcoils.circuit_of_each_coil{ii}, G.dima ,'exact');% Check if jj coil belongs to any active and considered circuit.
    G.Tca(ii,tmp) = 1;
end

else % If no connectin matrix is specified used all the coils
G.na = G.pf_active.pfcoils.dim; % Number of active pf circuits used 
G.Tca = eye(G.na, G.na);
end
G.Twa = G.Twc*G.Tca;

%% Flux loops
G.nf  = numel(IDS_inputs.magnetics.flux_loop); % number of flux loops -1 since the data used are the relative one
G.dimf = cell(G.nf,1); % names of the flux loops
G.rf = zeros(G.nf,1);
G.zf = zeros(G.nf,1);
%------------------------------------
%---- WARNING one could have many points to describe one single flux loop position. Need to check the FLUXLOOP position document and understand how they are defined
% ------------------------------------
for ii=1:G.nf
  G.dimf{ii} =  IDS_inputs.magnetics.flux_loop{ii}.name;
  G.rf(ii) = IDS_inputs.magnetics.flux_loop{ii}.position{1}.r;
  G.zf(ii) = IDS_inputs.magnetics.flux_loop{ii}.position{1}.z; 
end


%% Vessel
% --------------------------------------------------------------------------------------------------------------------------------------
% -----------THIS will have to be modified to take into account many different units for the vessel and different description_2d element
% ------------------------------------------------------------------------------------------------------------------------------------
G.nv = P.nv; % Number of vessel filaments
if strcmp( 'annular', strtrim(P.wall_type))
  % With the anular description I expect to have the flux loop voltages
  % in order to reconstruct the current in the vessel.
  % The resistence for each filament is computed below.
  
  % Vessel inner and outer R,Z description
  G.rv_inside  = IDS_inputs.wall.description_2d{1}.vessel.unit{1}.annular.outline_inner.r;
  G.zv_inside  = IDS_inputs.wall.description_2d{1}.vessel.unit{1}.annular.outline_inner.z;
  G.rv_outside = IDS_inputs.wall.description_2d{1}.vessel.unit{1}.annular.outline_outer.r; 
  G.zv_outside = IDS_inputs.wall.description_2d{1}.vessel.unit{1}.annular.outline_outer.z;
  
  [G.rv, G.zv] = theta_int(G.rv_inside, G.zv_inside, G.rv_outside, G.zv_outside, G.nv); % Interpolation
  G.ns = G.nf; % number of segment equal to the number of flux loops
  
  % Vessel filament resitences.
  G.eta = IDS_inputs.wall.description_2d{1}.vessel.unit{1}.annular.resistivity;
  G.Rv = eval_Rv(G.eta, G.rv, G.zv, G.rv_inside, G.zv_inside, G.rv_outside, G.zv_outside);
  
  % Grouping vessel coils into segments Is = T*Iv such that they share the same voltage
  [G.Tvs, G.T] = tcvvseg(G.Rv, G.nf, G.rf,G.zf,G.rv,G.zv);
  
  G.Rs = 1./ (G.T*(1./G.Rv)); % Rs^-1 = sum_{s \in v} Rv^-1
  %G.Rs2 = diag(G.Tvs'*diag(G.Rv)*G.Tvs); % This is equivalent to the
  %previous one

elseif strcmp( 'block',strtrim(P.wall_type))
  % With the block description I expect to have the current in each
  % block.
  % No need to compute the resistence of the filaments.
  % Need to discretized the blocks into filaments to compute the mutual
  % matrices
  
  % TODO extend the thet_int function to provide discretization for the
  % block type IDS_inputs.wall
  
  % First rude implementation
  % With a block description I would first evaluate the rzdrdz and than
  % use the exisiting tool to discretize the coils.
  % Once the rv are evaluated I would group them with the
  
  n_blocks =  numel(wall.wall2d.vessel.vessel_unit.blocks.blocks_unit);
  pfrzdrdz = zeros(n_blocks, 1);
  for ii=1:n_blocks
    [pfrzdrdz(ii, :)] = eval_pfrzdrdz(wall.wall2d.vessel.vessel_unit.blocks(ii).blocks_unit.position.r, wall.wall2d.vessel.vessel_unit.blocks(ii).blocks_unit.position.z );    
  end
  
  % Discretize as if they were rectangles
  [G.rv, G.zv, G.Tvs] = eval_trap(pfrzdrdz, G.n_turns_real, 2); 
elseif strcmp('none', strtrim(P.wall_type)) 
  % No description of the wall is available
  G.nv =  0;
  G.rv =  [];
  G.zv =  [];
  G.Tvs = [];
end
%% Tile apperture (limiter)
G.rl = G.wall.description_2d{1}.limiter.unit{1}.outline.r;
G.zl = G.wall.description_2d{1}.limiter.unit{1}.outline.z;

%% The mesh in taken in between the limiter and the vessel, if vessel is available
% -----------------------------------------------------------------
%------------ Need to get the dimension of the grid from somewhere else
% -------------------------------------------------------------
if strcmp('none', strtrim(P.wall_type))
G.rx = linspace( min(G.rl)-0.1, max(G.rl)+0.1, 28 )';
G.zx = linspace( min(G.zl)-0.1, max(G.zl)+0.1, 65 )';
else
G.rx = linspace( mean([min(G.rl),min(G.rv_inside)]), mean([max(G.rv_inside),max(G.rl)]),  28)';
G.zx = linspace( mean([min(G.zl),min(G.zv_inside)]), mean([max(G.zv_inside),max(G.zl)]),   65)';
end


%% Reinterpolation of limiter
[G.rl, G.zl] = plane_int(G.rl, G.zl, pdist(G.rx(2),1,G.rx(1),1)/4 ); % Linear interpolation, to make vessel discretization finer
G.nl = numel(G.zl);

% Poloidal magnetic probes
G.nm =  numel(IDS_inputs.magnetics.bpol_probe); 
G.dimm = cell(G.nm,1);
G.rm = zeros(G.nm,1);
G.zm = zeros(G.nm,1);
G.am = zeros(G.nm,1);
for ii=1:G.nm
  G.dimm{ii} =  IDS_inputs.magnetics.bpol_probe{ii}.name;
  G.rm(ii) = IDS_inputs.magnetics.bpol_probe{ii}.position.r;
  G.zm(ii) = IDS_inputs.magnetics.bpol_probe{ii}.position.z; 
  G.am(ii) = IDS_inputs.magnetics.bpol_probe{ii}.poloidal_angle;   
end

%% Evaluate the green function and correct the dimensions if needed
G = liug(G,P); 
return
