function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,varargin);
%
%  [ids_equilibrium,ids_equilibrium_description,varargout] = get_ids_equilibrium_fixed_boundary(shot,varargin);
% 
%

ids_equilibrium = struct([]);
ids_equilibrium_description = struct([]);

% initialize input parser
p = inputParser;
% no required inputs here so one can call with empty args and get defaults parameters
% effectively required inputs should have defaults as empty
p.addOptional('shot', [], @(x) isempty(x) || (isnumeric(x) && isscalar(x) && (x == round(x)))); % integer
p.addOptional('empty_struct', -1, @(x) isstruct(x));
p.addOptional('comment', 'default comment', @(x) isempty(x) || ischar(x)); % char

p.parse;
defaults_equilibrium = p.Results; % to keep track of defaults

if nargin==1
  p.parse(shot);
  params = p.Results;
elseif nargin>=2
  p.parse(shot,varargin{:});
  params = p.Results;
else
  p.parse;
  if nargout>=3; varargout{1} = rmfield(p.Results,'shot'); end
  return
end
% replace empty inputs with defaults
names = fieldnames(params);
mask = structfun(@isempty,params);
if any(mask),
  params = rmfield(params,unique([names(mask); p.UsingDefaults.']));
  if ~isfield(params, 'shot') || isnan(params.shot)
    warning('No shot entered');
    return
  end
  p.parse(params.shot,rmfield(params,'shot'));
  params = p.Results;
end
if nargout>=3; varargout{1} = rmfield(p.Results,'shot'); end

params_equilibrium = params;

ids_equilibrium = params_equilibrium.empty_struct;
clear ids_equilibrium_description % so can add new subfields

%
% ids_properties
%
ids_equilibrium.ids_properties.comment = params_equilibrium.comment;
ids_equilibrium.ids_properties.homogeneous_time = 1;
ids_equilibrium.ids_properties.source = ['TCV shot #' num2str(params_equilibrium.shot)];
ids_equilibrium.ids_properties.provider = 'get_ids_equilibrium_fixed_boundary';
ids_equilibrium.ids_properties.creation_date = date;

% As a general rule, for a new substructure under the main ids, construct a local structure like:
% "global_quantities" with subfields being the relevant data to get and a local structure:
% "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used
%
% vacuum_toroidal_field and time, using homogeneous
%
vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','source','liuqe'); % to get on liuqe time array
vacuum_toroidal_field_desc.b0 = '''b0'',''source'',''liuqe''';
vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'',''source'',''liuqe'']';
ids_equilibrium.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
ids_equilibrium.vacuum_toroidal_field.b0 = vacuum_toroidal_field.b0.data;
ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;

ids_equilibrium.time = vacuum_toroidal_field.b0.t;
ids_equilibrium_description.time = '.t subfield from: [''b0'',''source'',''liuqe'']';

ids_equilibrium.time_slice(1:length(ids_equilibrium.time)) = ids_equilibrium.time_slice(1);

% load time array data to copy to time_slices

% global_quantities data into local global_quantities.* structure with correct end names and global_quantities_desc.* with description. Use temp.* and temp_desc.* structures for temporary data

% brute force solution load all eqdsks
% $$$ for it=1:length(ids_equilibrium.time)
% $$$   ids_equilibrium.time(it)
% $$$   temp.eqdsks{it}=gdat(params_equilibrium.shot,'eqdsk','time',ids_equilibrium.time(it),'write',0);
% $$$ end
% $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_equilibrium.time(it)';

global_quantities.area = gdat(params_equilibrium.shot,'area_edge');
global_quantities_desc.area = 'area_edge';
global_quantities.beta_normal = gdat(params_equilibrium.shot,'betan');
global_quantities_desc.beta_normal = 'betan';
global_quantities.beta_pol = gdat(params_equilibrium.shot,'betap');
global_quantities_desc.beta_pol = 'betap';
global_quantities.beta_tor = gdat(params_equilibrium.shot,'beta');
global_quantities_desc.beta_tor = 'beta';
global_quantities.energy_mhd = gdat(params_equilibrium.shot,'w_mhd');
global_quantities_desc.energy_mhd = 'w_mhd';
global_quantities.ip = gdat(params_equilibrium.shot,'ip');
global_quantities_desc.ip = 'ip';
% length_pol = gdat(params_equilibrium.shot,'length_pol'); % to be added
global_quantities.li_3 = gdat(params_equilibrium.shot,'li');
global_quantities_desc.li_3 = 'li';
temp.r_magnetic_axis = gdat(params_equilibrium.shot,'r_axis');
temp_desc.r_magnetic_axis = 'r_axis';
temp.z_magnetic_axis = gdat(params_equilibrium.shot,'z_axis');
temp_desc.z_magnetic_axis = 'z_axis';
temp.psi_axis = gdat(params_equilibrium.shot,'psi_axis'); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value
temp_desc.psi_axis = 'psi_axis';
global_quantities.psi_boundary = gdat(params_equilibrium.shot,'psi_edge');
global_quantities_desc.psi_boundary = 'psi_edge';
global_quantities.q_95 = gdat(params_equilibrium.shot,'q95');
global_quantities_desc.q_95 = 'q95';
global_quantities.q_axis = gdat(params_equilibrium.shot,'q0'); % will be checked with q_rho?
global_quantities_desc.q_axis = 'q0';
temp.q_rho = gdat(params_equilibrium.shot,'q_rho');
temp_desc.q_rho = 'q_rho';
% surface = gdat(params_equilibrium.shot,'surface'); % to be added
global_quantities.volume = gdat(params_equilibrium.shot,'volume');
global_quantities_desc.volume = 'volume';
global_quantities.w_mhd = gdat(params_equilibrium.shot,'w_mhd');
global_quantities_desc.w_mhd = 'w_mhd';

global_quantities_fieldnames = fieldnames(global_quantities);
special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
for it=1:length(ids_equilibrium.time)
  for i=1:length(global_quantities_fieldnames)
    if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
      if ~isstruct(ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}))
        ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = ...
            global_quantities.(global_quantities_fieldnames{i}).data(it);
      else
        special_fields{end+1} = global_quantities_fieldnames{i};
      end
    end
  end
end

% special case
for it=1:length(ids_equilibrium.time)
  ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r = temp.r_magnetic_axis.data(it);
  ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.z = temp.z_magnetic_axis.data(it);
  ids_equilibrium.time_slice{it}.global_quantities.psi_axis = temp.psi_axis.data(it) + ...
      ids_equilibrium.time_slice{it}.global_quantities.psi_boundary;
  [zz,izz]=min(temp.q_rho.data(:,it));
  ids_equilibrium.time_slice{it}.global_quantities.q_min.value = zz;
  ids_equilibrium.time_slice{it}.global_quantities.q_min.rho_tor_norm = temp.q_rho.grids_1d.rhotornorm(izz);
end

% for boundary in addition to lcfs
% active_limiter_point = gdat(params_equilibrium.shot,'active_limiter_point');
boundary.elongation = gdat(params_equilibrium.shot,'kappa');
boundary_desc.elongation = 'kappa';
% elongation_lower = gdat(params_equilibrium.shot,'elongation_lower');
% elongation_upper = gdat(params_equilibrium.shot,'elongation_upper');
boundary.minor_radius = gdat(params_equilibrium.shot,'a_minor');
boundary_desc.minor_radius = 'a_minor';
% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner');
% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer');
% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner');
% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer');
% strike_point = gdat(params_equilibrium.shot,'strike_point');
boundary.triangularity = gdat(params_equilibrium.shot,'delta');
boundary_desc.triangularity = 'delta';
boundary.triangularity_lower = gdat(params_equilibrium.shot,'delta_bottom');
boundary_desc.triangularity_lower = 'delta_bottom';
boundary.triangularity_upper = gdat(params_equilibrium.shot,'delta_top');
boundary_desc.triangularity_upper = 'delta_top';
temp.n_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''n_xpts'''',''''liuqe.m'''')');
temp_desc.n_x_point = '''tcv_eq(''''n_xpts'''',''''liuqe.m'''')''';
temp.r_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''r_xpts'''',''''liuqe.m'''')');
temp_desc.r_x_point = '''tcv_eq(''''r_xpts'''',''''liuqe.m'''')''';
temp.z_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''z_xpts'''',''''liuqe.m'''')');
temp_desc.z_x_point = '''tcv_eq(''''z_xpts'''',''''liuqe.m'''')''';
temp.rgeom = gdat(params_equilibrium.shot,'rgeom');
temp_desc.rgeom = 'rgeom';
temp.zgeom = gdat(params_equilibrium.shot,'zgeom');
temp_desc.zgeom = 'zgeom';
temp.r_lcfs = gdat(params_equilibrium.shot,'r_contour_edge');
temp_desc.r_lcfs = 'r_contour_edge';
temp.z_lcfs = gdat(params_equilibrium.shot,'z_contour_edge');
temp_desc.z_lcfs = 'z_contour_edge';

boundary_fieldnames = fieldnames(boundary);
special_fields = {'lcfs', 'geometric_axis', 'x_point'}; % fields needing non-automatic treatments
for it=1:length(ids_equilibrium.time)
  for i=1:length(boundary_fieldnames)
    if ~any(strcmp(boundary_fieldnames{i},special_fields))
      if ~isstruct(ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}))
        ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}) = ...
            boundary.(boundary_fieldnames{i}).data(it);
      else
        special_fields{end+1} = boundary_fieldnames{i};
      end
    end
  end
end

% special cases
for it=1:length(ids_equilibrium.time)
  ids_equilibrium.time_slice{it}.boundary.outline.r = temp.r_lcfs.data(:,it);
  ids_equilibrium.time_slice{it}.boundary.outline.z = temp.z_lcfs.data(:,it);
  ids_equilibrium.time_slice{it}.boundary.lcfs.r = ids_equilibrium.time_slice{it}.boundary.outline.r;
  ids_equilibrium.time_slice{it}.boundary.lcfs.z = ids_equilibrium.time_slice{it}.boundary.outline.z;
  ids_equilibrium.time_slice{it}.boundary.geometric_axis.r = temp.rgeom.data(it);
  ids_equilibrium.time_slice{it}.boundary.geometric_axis.z = temp.zgeom.data(it);
  if temp.n_x_point.data(it) > 0
    ids_equilibrium.time_slice{it}.boundary.x_point(1:temp.n_x_point.data(it)) = ids_equilibrium.time_slice{it}.boundary.x_point(1);
    for i=1:length(temp.n_x_point.data(it))
      ids_equilibrium.time_slice{it}.boundary.x_point{i}.r = temp.r_x_point.data(i,it);
      ids_equilibrium.time_slice{it}.boundary.x_point{i}.z = temp.z_x_point.data(i,it);
    end
  end
end

%
%% profiles_1d (cannot use eqdsk since not same radial mesh)
%
% area = gdat(params_equilibrium.shot,'area');
% b_average = gdat(params_equilibrium.shot,'b_average');
% beta_pol = gdat(params_equilibrium.shot,'beta_pol');
% b_field_average = gdat(params_equilibrium.shot,'b_field_average');
% b_field_max = gdat(params_equilibrium.shot,'b_field_max');
% b_field_min = gdat(params_equilibrium.shot,'b_field_min');
% b_max = gdat(params_equilibrium.shot,'b_max');
% b_min = gdat(params_equilibrium.shot,'b_min');
% darea_dpsi = gdat(params_equilibrium.shot,'darea_dpsi');
% darea_drho_tor = gdat(params_equilibrium.shot,'darea_drho_tor');
profiles_1d.dpressure_dpsi = gdat(params_equilibrium.shot,'pprime');
% dpsi_drho_tor = gdat(params_equilibrium.shot,'dpsi_drho_tor');
% dvolume_dpsi = gdat(params_equilibrium.shot,'dvolume_dpsi');
% dvolume_drho_tor = gdat(params_equilibrium.shot,'dvolume_drho_tor');
% elongation = gdat(params_equilibrium.shot,'elongation');
profiles_1d.f_df_dpsi = gdat(params_equilibrium.shot,'ttprime');
profiles_1d.f = gdat(params_equilibrium.shot,'rbphi_rho');
% geometric_axis = gdat(params_equilibrium.shot,'geometric_axis');
% gm1 = gdat(params_equilibrium.shot,'gm1');
% gm2 = gdat(params_equilibrium.shot,'gm2');
% gm3 = gdat(params_equilibrium.shot,'gm3');
% gm4 = gdat(params_equilibrium.shot,'gm4');
% gm5 = gdat(params_equilibrium.shot,'gm5');
% gm6 = gdat(params_equilibrium.shot,'gm6');
% gm7 = gdat(params_equilibrium.shot,'gm7');
% gm8 = gdat(params_equilibrium.shot,'gm8');
% gm9 = gdat(params_equilibrium.shot,'gm9');
% j_parallel = gdat(params_equilibrium.shot,'j_parallel');
% j_tor = gdat(params_equilibrium.shot,'j_tor');
% magnetic_shear = gdat(params_equilibrium.shot,'magnetic_shear');
% mass_density = gdat(params_equilibrium.shot,'mass_density');
profiles_1d.phi = gdat(params_equilibrium.shot,'phi_tor');
profiles_1d.pressure = gdat(params_equilibrium.shot,'pressure');
% psi = gdat(params_equilibrium.shot,'psi_rho'); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
profiles_1d.q = gdat(params_equilibrium.shot,'q_rho');
profiles_1d.rho_tor = gdat(params_equilibrium.shot,'rhotor');
%rho_tor_norm = gdat(params_equilibrium.shot,'rhotor_norm'); % from rho_tor
profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,'rhovol');
% r_inboard = gdat(params_equilibrium.shot,'r_inboard');
% r_outboard = gdat(params_equilibrium.shot,'r_outboard');
% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner');
% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer');
% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner');
% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer');
% surface = gdat(params_equilibrium.shot,'surface');
% trapped_fraction = gdat(params_equilibrium.shot,'trapped_fraction');
% triangularity_lower = gdat(params_equilibrium.shot,'triangularity_lower');
% triangularity_upper = gdat(params_equilibrium.shot,'triangularity_upper');
profiles_1d.volume = gdat(params_equilibrium.shot,'volume_rho');

profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
for it=1:length(ids_equilibrium.time)
  for i=1:length(profiles_1d_fieldnames)
    if ~any(strcmp(profiles_1d_fieldnames{i},special_fields))
      if ~isstruct(ids_equilibrium.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}))
        if ~ischar(profiles_1d.(profiles_1d_fieldnames{i}).data) && ~isempty(profiles_1d.(profiles_1d_fieldnames{i}).data) ...
              && size(profiles_1d.(profiles_1d_fieldnames{i}).data,2)>=it
          ids_equilibrium.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}) = ...
              profiles_1d.(profiles_1d_fieldnames{i}).data(:,it);
        end
      else
        special_fields{end+1} = profiles_1d_fieldnames{i};
      end
    end
  end
end

% special cases
nrho = length(profiles_1d.rho_tor.x);
ntime = length(temp.psi_axis.data);
for it=1:length(ids_equilibrium.time)
  ids_equilibrium.time_slice{it}.profiles_1d.rho_tor_norm = ids_equilibrium.time_slice{it}.profiles_1d.rho_tor./ ...
      ids_equilibrium.time_slice{it}.profiles_1d.rho_tor(end);
  ids_equilibrium.time_slice{it}.profiles_1d.psi = (1-profiles_1d.rho_tor.x.^2).*temp.psi_axis.data(it) + ...
      global_quantities.psi_boundary.data(it);
end

%
%% profiles_2d{1} ala eqdsk, only this one thus grid_type=1
%
% b_field_r = gdat(params_equilibrium.shot,'b_field_r');
% b_field_tor = gdat(params_equilibrium.shot,'b_field_tor');
% b_field_z = gdat(params_equilibrium.shot,'b_field_z');
% b_r = gdat(params_equilibrium.shot,'b_r');
% b_tor = gdat(params_equilibrium.shot,'b_tor');
% b_z = gdat(params_equilibrium.shot,'b_z');
% grid = gdat(params_equilibrium.shot,'grid'); % special
profiles_2d.grid_type.name = 'rectangular';
profiles_2d.grid_type.index = 1;
profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk';
% j_parallel = gdat(params_equilibrium.shot,'j_parallel');
% j_tor = gdat(params_equilibrium.shot,'j_tor');
% phi = gdat(params_equilibrium.shot,'phi');
profiles_2d.psi = gdat(params_equilibrium.shot,'psi');
% r = gdat(params_equilibrium.shot,'r'); % not to be filled since in grid.dim1
% theta = gdat(params_equilibrium.shot,'theta');
% z = gdat(params_equilibrium.shot,'z'); % not to be filled since in grid.dim2

profiles_2d_fieldnames = fieldnames(profiles_2d);
special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments
for it=1:length(ids_equilibrium.time)
  for i=1:length(profiles_2d_fieldnames)
    if ~any(strcmp(profiles_2d_fieldnames{i},special_fields))
      if ~isstruct(ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}))
        if ~ischar(profiles_2d.(profiles_2d_fieldnames{i}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i}).data) ...
              && size(profiles_2d.(profiles_2d_fieldnames{i}).data,2)>=it
          ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}) = ...
              profiles_2d.(profiles_2d_fieldnames{i}).data(:,:,it);
        end
      else
        special_fields{end+1} = profiles_2d_fieldnames{i};
      end
    end
  end
end

% special cases
for it=1:length(ids_equilibrium.time)
  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name;
  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.index = profiles_2d.grid_type.index;
  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.description = profiles_2d.grid_type.description;
  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid.dim1 = profiles_2d.psi.dim{1};
  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid.dim2 = profiles_2d.psi.dim{2};
end
