function [ids_equilibrium,ids_equilibrium_description,varargout] = aug_get_ids_equilibrium(shot,ids_equil_empty,varargin);
%
%  [ids_equilibrium,ids_equilibrium_description,varargout] = aug_get_ids_equilibrium(shot,ids_equil_empty,varargin);
% 
%

machine = 'aug';

[ids_equilibrium, params_equilibrium] = aug_ids_headpart(shot,ids_equil_empty,'equilibrium',varargin{:});

% 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 and q_rho time base
%
% time may be "too long" in shotfile compared to effective equilibrium reconstructed times

% use "equil" for AUG as reference at this stage, it contains most info, so no need for individual calls
temp_2d.equil = gdat(params_equilibrium.shot,'equil','machine',machine);
R0exp = 1.65; % to decide B0 vacuum from Jpol

% temp.q_rho = gdat(params_equilibrium.shot,'q_rho','machine',machine);
ids_equilibrium.time = temp_2d.equil.t;
ids_equilibrium_description.time = '.t subfield from: [''equil'']';
%
tens_time = -1;

% $$$ vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','machine',machine);
% $$$ vacuum_toroidal_field_desc.b0 = ['''b0'' with interpos using tens_time = ' num2str(tens_time)];
% $$$ vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'']';
% $$$ ids_equilibrium.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
% $$$ ids_equilibrium.vacuum_toroidal_field.b0 = interpos(vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data, ids_equilibrium.time,tens_time);
% $$$ ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
vacuum_toroidal_field.b0.data = temp_2d.equil.jpol(end,:).*2e-7 / R0exp;
vacuum_toroidal_field.b0.t = temp_2d.equil.t;
vacuum_toroidal_field_desc.b0 = ['''equil.jpol'' *2e-7 divided by R0exp=' num2str(R0exp)];
vacuum_toroidal_field_desc.r0 = 'set by hand in aug_get_ids_equilibrium';
ids_equilibrium.vacuum_toroidal_field.r0 = R0exp;
ids_equilibrium.vacuum_toroidal_field.b0 = vacuum_toroidal_field.b0.data;
ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;

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,'machine',machine);
% $$$ end
% $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_equilibrium.time(it)';

global_quantities.area.data = temp_2d.equil.area(end,:);
global_quantities.area.t = temp_2d.equil.t;
global_quantities_desc.area = 'area_edge';
global_quantities.beta_normal = gdat(params_equilibrium.shot,'betan','machine',machine);
global_quantities_desc.beta_normal = 'betan';
global_quantities.beta_pol = gdat(params_equilibrium.shot,'betap','machine',machine);
global_quantities_desc.beta_pol = 'betap';
global_quantities.beta_tor = gdat(params_equilibrium.shot,'beta','machine',machine);
global_quantities_desc.beta_tor = 'beta';
global_quantities.energy_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',machine);
global_quantities_desc.energy_mhd = 'w_mhd';
global_quantities.ip = gdat(params_equilibrium.shot,'ip','machine',machine);
global_quantities_desc.ip = 'ip';
% length_pol = gdat(params_equilibrium.shot,'length_pol','machine',machine); % to be added
global_quantities.li_3 = gdat(params_equilibrium.shot,'li','machine',machine);
global_quantities_desc.li_3 = 'li';
temp.r_magnetic_axis = gdat(params_equilibrium.shot,'r_axis','machine',machine);
temp_desc.r_magnetic_axis = 'r_axis';
temp.z_magnetic_axis = gdat(params_equilibrium.shot,'z_axis','machine',machine);
temp_desc.z_magnetic_axis = 'z_axis';
global_quantities.psi_axis = gdat(params_equilibrium.shot,'psi_axis','machine',machine);
global_quantities_desc.psi_axis = 'psi_axis';
global_quantities.psi_boundary = gdat(params_equilibrium.shot,'psi_edge','machine',machine);
global_quantities_desc.psi_boundary = 'psi_edge';
global_quantities.q_95 = gdat(params_equilibrium.shot,'q95','machine',machine);
global_quantities_desc.q_95 = 'q95';
global_quantities.q_axis = gdat(params_equilibrium.shot,'q0','machine',machine); % will be checked with q_rho?
global_quantities_desc.q_axis = 'q0';
% surface = gdat(params_equilibrium.shot,'surface','machine',machine); % to be added
global_quantities.volume = gdat(params_equilibrium.shot,'volume','machine',machine);
global_quantities_desc.volume = 'volume';
global_quantities.w_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',machine);
global_quantities_desc.w_mhd = 'w_mhd';

global_quantities_fieldnames = fieldnames(global_quantities);
special_fields = {'magnetic_axis', 'q_min'}; % fields needing non-automatic treatments
for i=1:length(global_quantities_fieldnames)
  try
    if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
      if ~isstruct(ids_equilibrium.time_slice{1}.global_quantities.(global_quantities_fieldnames{i}))
	if length(global_quantities.(global_quantities_fieldnames{i}).data) > length(ids_equilibrium_description.time)
	  % assume need interpolation with constant extrapolation
	  aa = interpos(63,global_quantities.(global_quantities_fieldnames{i}).t, ...
	  global_quantities.(global_quantities_fieldnames{i}).data,ids_equilibrium.time,tens_time);
	  for it=1:length(ids_equilibrium.time)
	    ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = aa(it);
	  end
	elseif length(global_quantities.(global_quantities_fieldnames{i}).data) == length(ids_equilibrium_description.time)
	  % assume same time
	  if sum(abs(global_quantities.(global_quantities_fieldnames{i}).t-ids_equilibrium_description.time))>0
	    warning(['not same time for ' global_quantities_fieldnames{i}])
	  else
	    for it=1:length(ids_equilibrium.time)
	      ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = aa(it);
	    end
	  end
	elseif isempty(global_quantities.(global_quantities_fieldnames{i}).data)
	  % do not do anything, no data
	else
	  warning(['less points than expected, cannot treat for ' global_quantities_fieldnames{i}])
	end
      else
	warning([global_quantities_fieldnames{i} ' is a structure should be in special cases?']);
      end
    else
      special_fields{end+1} = global_quantities_fieldnames{i};
    end
  catch ME_global
    rethrow(ME_global)
    keyboard
  end
end
% special case
aar = interpos(63,temp.r_magnetic_axis.t,temp.r_magnetic_axis.data,ids_equilibrium.time,tens_time);
aaz = interpos(63,temp.z_magnetic_axis.t,temp.z_magnetic_axis.data,ids_equilibrium.time,tens_time);
for it=1:length(ids_equilibrium.time)
  ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r = aar(it);
  ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.z = aaz(it);
  [zz,izz]=min(temp_2d.equil.qvalue(:,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_2d.equil.rhotornorm(izz,it);
end

% for boundary in addition to lcfs
% active_limiter_point = gdat(params_equilibrium.shot,'active_limiter_point','machine',machine);
boundary.elongation = gdat(params_equilibrium.shot,'kappa','machine',machine);
boundary_desc.elongation = 'kappa';
% elongation_lower = gdat(params_equilibrium.shot,'elongation_lower','machine',machine);
% elongation_upper = gdat(params_equilibrium.shot,'elongation_upper','machine',machine);
boundary.minor_radius = gdat(params_equilibrium.shot,'a_minor','machine',machine);
boundary_desc.minor_radius = 'a_minor';
% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',machine);
% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer','machine',machine);
% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner','machine',machine);
% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer','machine',machine);
% strike_point = gdat(params_equilibrium.shot,'strike_point','machine',machine);
boundary.triangularity = gdat(params_equilibrium.shot,'delta','machine',machine);
boundary_desc.triangularity = 'delta';
boundary.triangularity_lower = gdat(params_equilibrium.shot,'delta_bottom','machine',machine);
boundary_desc.triangularity_lower = 'delta_bottom';
boundary.triangularity_upper = gdat(params_equilibrium.shot,'delta_top','machine',machine);
boundary_desc.triangularity_upper = 'delta_top';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TCV LIUQE specific, to deal later...
% $$$ temp.n_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''n_xpts'''',''''liuqe.m'''')','machine',machine);
% $$$ 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'''')','machine',machine);
% $$$ 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'''')','machine',machine);
% $$$ temp_desc.z_x_point = '''tcv_eq(''''z_xpts'''',''''liuqe.m'''')''';
temp.rgeom = gdat(params_equilibrium.shot,'rgeom','machine',machine);
temp_desc.rgeom = 'rgeom';
temp.zgeom = gdat(params_equilibrium.shot,'zgeom','machine',machine);
temp_desc.zgeom = 'zgeom';
% $$$ temp.r_lcfs = gdat(params_equilibrium.shot,'r_contour_edge','machine',machine);
% $$$ temp_desc.r_lcfs = 'r_contour_edge';
% $$$ temp.z_lcfs = gdat(params_equilibrium.shot,'z_contour_edge','machine',machine);
% $$$ temp_desc.z_lcfs = 'z_contour_edge';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% at this stage use geteqdsk to get R,Z boundary since already included
for it=1:length(temp_2d.equil.t)
  [eqdskAUG{it}, dummy, equil_t_index]=geteqdskAUG(temp_2d.equil,temp_2d.equil.t(it),[],[],[],[],[],'fignb',999);
end

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

% special cases
for it=1:length(ids_equilibrium.time)
  if ~isempty(eqdskAUG{it}.rplas)
    ids_equilibrium.time_slice{it}.boundary.outline.r = eqdskAUG{it}.rplas;
    ids_equilibrium.time_slice{it}.boundary.outline.z = eqdskAUG{it}.zplas;
    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;
  end
end
if ~isempty(temp.rgeom.data)
  aa = interpos(63,temp.rgeom.t,temp.rgeom.data,ids_equilibrium.time,tens_time);
  for it=1:length(ids_equilibrium.time)
    ids_equilibrium.time_slice{it}.boundary.geometric_axis.r = aa(it);
  end
end
if ~isempty(temp.zgeom.data)
  aa = interpos(63,temp.zgeom.t,temp.zgeom.data,ids_equilibrium.time,tens_time);
  for it=1:length(ids_equilibrium.time)
    ids_equilibrium.time_slice{it}.boundary.geometric_axis.z = aa(it);
  end
end
% could use temp_2d.equil.Xpoints but not sure which points to take
% $$$   if ~isempty(temp.n_x_point.data) && 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

%
%% profiles_1d (cannot use eqdsk since not same radial mesh), but can interpolate and use base from equil
%
% area = gdat(params_equilibrium.shot,'area','machine',machine);
area.data = temp_2d.equil.area;
% b_average = gdat(params_equilibrium.shot,'b_average','machine',machine);
% obsolescent: b_average.data = temp_2d.equil.bave;
% (profile) beta_pol = gdat(params_equilibrium.shot,'','machine',machine);
b_field_average = temp_2d.equil.bave;
% b_field_max = gdat(params_equilibrium.shot,'b_field_max','machine',machine);
% b_field_min = gdat(params_equilibrium.shot,'b_field_min','machine',machine);
% b_max = gdat(params_equilibrium.shot,'b_max','machine',machine);
% b_min = gdat(params_equilibrium.shot,'b_min','machine',machine);
% darea_dpsi = gdat(params_equilibrium.shot,'darea_dpsi','machine',machine);
% darea_drho_tor = gdat(params_equilibrium.shot,'darea_drho_tor','machine',machine);
profiles_1d.dpressure_dpsi.data = temp_2d.equil.dpressuredpsi;
% dpsi_drho_tor = gdat(params_equilibrium.shot,'dpsi_drho_tor','machine',machine);
dvolume_dpsi = temp_2d.equil.dvoldpsi;
% dvolume_drho_tor = gdat(params_equilibrium.shot,'dvolume_drho_tor','machine',machine);
% elongation = gdat(params_equilibrium.shot,'elongation','machine',machine);
profiles_1d.f_df_dpsi.data = temp_2d.equil.ffprime;
profiles_1d.f.data = temp_2d.equil.jpol*2e-7;
% geometric_axis = gdat(params_equilibrium.shot,'geometric_axis','machine',machine);
% gm1 = gdat(params_equilibrium.shot,'gm1','machine',machine);
% gm2 = gdat(params_equilibrium.shot,'gm2','machine',machine);
% gm3 = gdat(params_equilibrium.shot,'gm3','machine',machine);
% gm4 = gdat(params_equilibrium.shot,'gm4','machine',machine);
% gm5 = gdat(params_equilibrium.shot,'gm5','machine',machine);
% gm6 = gdat(params_equilibrium.shot,'gm6','machine',machine);
% gm7 = gdat(params_equilibrium.shot,'gm7','machine',machine);
% gm8 = gdat(params_equilibrium.shot,'gm8','machine',machine);
% gm9 = gdat(params_equilibrium.shot,'gm9','machine',machine);
% j_parallel = gdat(params_equilibrium.shot,'j_parallel','machine',machine);
% j_tor = gdat(params_equilibrium.shot,'j_tor','machine',machine);
% magnetic_shear = gdat(params_equilibrium.shot,'magnetic_shear','machine',machine);
% mass_density = gdat(params_equilibrium.shot,'mass_density','machine',machine);
profiles_1d.phi.data = temp_2d.equil.pressure;
profiles_1d.pressure.data = temp_2d.equil.pressure;
% psi = gdat(params_equilibrium.shot,'psi_rho','machine',machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
profiles_1d.psi.data = temp_2d.equil.psi;
profiles_1d.q.data = temp_2d.equil.qvalue;
for it=1:length(temp_2d.equil.t)
  profiles_1d.rho_tor.data(:,it) = sqrt(temp_2d.equil.phi(:,it)./vacuum_toroidal_field.b0.data(it)/pi);
end
%rho_tor_norm = gdat(params_equilibrium.shot,'rhotor_norm','machine',machine); % from rho_tor
profiles_1d.rho_volume_norm.data = temp_2d.equil.rhovolnorm;
% r_inboard = gdat(params_equilibrium.shot,'r_inboard','machine',machine);
% r_outboard = gdat(params_equilibrium.shot,'r_outboard','machine',machine);
% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',machine);
% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer','machine',machine);
% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner','machine',machine);
% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer','machine',machine);
% surface = gdat(params_equilibrium.shot,'surface','machine',machine);
% trapped_fraction = gdat(params_equilibrium.shot,'trapped_fraction','machine',machine);
% triangularity_lower = gdat(params_equilibrium.shot,'triangularity_lower','machine',machine);
% triangularity_upper = gdat(params_equilibrium.shot,'triangularity_upper','machine',machine);
profiles_1d.volume.data = temp_2d.equil.vol;

profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
special_fields = {'geometric_axis', 'rho_tor_norm'}; % 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
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 = temp.psi_axis.data(it) + ...
  %     profiles_1d.rho_tor.x(:,it).^2.*(global_quantities.psi_boundary.data(it)-temp.psi_axis.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','machine',machine);
% b_field_tor = gdat(params_equilibrium.shot,'b_field_tor','machine',machine);
% b_field_z = gdat(params_equilibrium.shot,'b_field_z','machine',machine);
% b_r = gdat(params_equilibrium.shot,'b_r','machine',machine);
% b_tor = gdat(params_equilibrium.shot,'b_tor','machine',machine);
% b_z = gdat(params_equilibrium.shot,'b_z','machine',machine);
% grid = gdat(params_equilibrium.shot,'grid','machine',machine); % 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','machine',machine);
% j_tor = gdat(params_equilibrium.shot,'j_tor','machine',machine);
% phi = gdat(params_equilibrium.shot,'phi','machine',machine);
% profiles_2d.psi = gdat(params_equilibrium.shot,'psi','machine',machine); % add psi_bound in a second step in special cases
% r = gdat(params_equilibrium.shot,'r','machine',machine); % not to be filled since in grid.dim1
% theta = gdat(params_equilibrium.shot,'theta','machine',machine);
% z = gdat(params_equilibrium.shot,'z','machine',machine); % not to be filled since in grid.dim2

profiles_2d_fieldnames = fieldnames(profiles_2d);
special_fields = {'grid', 'grid_type', 'psi'}; % 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,3)>=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 = temp_2d.equil.Rmesh(:,it);
  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid.dim2 = temp_2d.equil.Zmesh(:,it);
  ids_equilibrium.time_slice{it}.profiles_2d{1}.psi(:,:) =  temp_2d.equil.psi2D(:,:,it);
end
