function [eqdsk_out,varargout] = ids2eqdsk(ids_equilibrium,ids_wall,time_ref,eqdsk_filename_suffix,write_opt)
%
% [eqdsk_out, varout] = ids2eqdsk(ids_equilibrium, ids_wall, time_ref, eqdsk_filename_suffix);
%
% ids_equilibrium: ids equilibrium structure to extract eqdsk related quantities
% ids_wall:        ids wall structure to extract limiter and vessel points if provided
% time_ref:        time to extract data
% eqdsk_filename_suffix: string to be used for end of eqdsk filename, e.g. shot_run_text
% write_opt:       1 (default) do write eqdsk to /tmp/username, 0: no writing
%
% OUTPUT
% eqdsk_out:       eqdsk structure
%
if ~exist('time_ref','var') || isempty(time_ref)
  % take 1st time point
  if ids_equilibrium.ids_properties.homogeneous_time == 1
    time_ref = ids_equilibrium.time(1);
  else
    time_ref = ids_equilibrium.time_slice{1}.time;
  end
end
if ~exist('eqdsk_filename','var') || isempty(eqdsk_filename)
  eqdsk_filename = ['t' num2str(time_ref)];
end
if ~exist('write_opt','var') || isempty(write_opt)
  write_opt = 1;
end

% find time_slice for equilibrium ids
if ids_equilibrium.ids_properties.homogeneous_time == 1
  i_equil = iround_os(ids_equilibrium.time,time_ref);
else
  for j=1:numel(ids_equilibrium.time_slice{j})
    time_array(j) = ids_equilibrium.time_slice{j}.time;
  end
  i_equil = iround_os(time_array,time_ref);
end
% find time_slice for wall ids
if ids_wall.ids_properties.homogeneous_time == 1
  i_wall = iround_os(ids_wall.time,time_ref);
else
  for j=1:numel(ids_wall.time_slice{j})
    time_array(j) = ids_wall.time_slice{j}.time;
  end
  i_wall = iround_os(time_array,time_ref);
end

eqdsk_out.cocos = 11; % assume 11 within IMAS at this stage, should come from ids_equilibrium metadata

% check if needed fields are present and non empty
if ids_equilibrium.time_slice{i_equil}.global_quantities.ip ~= -9.0000e+40
  eqdsk_out.ip = ids_equilibrium.time_slice{i_equil}.global_quantities.ip;
else
  error(['ids_equilibrium.time_slice{i_equil}.global_quantities.ip not provided']);
end
if ids_equilibrium.vacuum_toroidal_field.b0 ~= -9.0000e+40
  eqdsk_out.b0 = ids_equilibrium.vacuum_toroidal_field.b0;
else
  error(['ids_equilibrium.vacuum_toroidal_field.b0 not provided']);
end
if ids_equilibrium.vacuum_toroidal_field.r0 ~= -9.0000e+40
  eqdsk_out.r0 = ids_equilibrium.vacuum_toroidal_field.r0;
else
  error(['ids_equilibrium.vacuum_toroidal_field.r0 not provided']);
end
if ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.r ~= -9.0000e+40
  eqdsk_out.raxis = ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.r;
else
  error(['ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.r not provided']);
end
if ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.z ~= -9.0000e+40
  eqdsk_out.zaxis = ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.z;
else
  error(['ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.z not provided']);
end
eqdsk_out.zshift = 0.; % can add an option to shift eqdsk but can use eqdsk_transform

if ids_equilibrium.time_slice{i_equil}.profiles_1d.dpressure_dpsi ~= -9.0000e+40
  dpressure_dpsi = ids_equilibrium.time_slice{i_equil}.profiles_1d.dpressure_dpsi;
else
  error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.dpressure_dpsi not provided']);
end
if ids_equilibrium.time_slice{i_equil}.profiles_1d.pressure ~= -9.0000e+40
  pressure = ids_equilibrium.time_slice{i_equil}.profiles_1d.pressure;
else
  error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.pressure not provided']);
end
if ids_equilibrium.time_slice{i_equil}.profiles_1d.f ~= -9.0000e+40
  f_rbphi = ids_equilibrium.time_slice{i_equil}.profiles_1d.f;
else
  error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.f not provided']);
end
if ids_equilibrium.time_slice{i_equil}.profiles_1d.f_df_dpsi ~= -9.0000e+40
  f_df_dpsi = ids_equilibrium.time_slice{i_equil}.profiles_1d.f_df_dpsi;
else
  error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.f_df_dpsi not provided']);
end
if ids_equilibrium.time_slice{i_equil}.global_quantities.psi_axis ~= -9.0000e+40
  psi_axis = ids_equilibrium.time_slice{i_equil}.global_quantities.psi_axis;
else
  error(['ids_equilibrium.time_slice{i_equil}.global_quantities.psi_axis not provided']);
end
if ids_equilibrium.time_slice{i_equil}.global_quantities.psi_boundary ~= -9.0000e+40
  psi_boundary = ids_equilibrium.time_slice{i_equil}.global_quantities.psi_boundary;
else
  error(['ids_equilibrium.time_slice{i_equil}.global_quantities.psi_boundary not provided']);
end
if ids_equilibrium.time_slice{i_equil}.profiles_1d.psi ~= -9.0000e+40
  psi = ids_equilibrium.time_slice{i_equil}.profiles_1d.psi;
  rhopol_norm = sqrt((psi-psi(1))./(psi(end)-psi(1)));
  rhopol_norm(1) = 0.; rhopol_norm(end) = 1.;
else
  error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.psi not provided, could use rho_tor_norm, not implemented yet']);
end
if ids_equilibrium.time_slice{i_equil}.profiles_1d.q ~= -9.0000e+40
  q_profile = ids_equilibrium.time_slice{i_equil}.profiles_1d.q;
else
  error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.q not provided']);
end
if numel(ids_equilibrium.time_slice{i_equil}.profiles_2d) >= 1
  i_2d = 0;
  % find eqdsk type grid
  for j=1:numel(ids_equilibrium.time_slice{i_equil}.profiles_2d)
    if ids_equilibrium.time_slice{i_equil}.profiles_2d{j}.grid_type.index == 1
      i_2d = j;
    end
  end
  if i_2d > 0
    % grid ala eqdsk so can use psi, dim1, dim2
    rmesh = ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.grid.dim1;
    zmesh = ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.grid.dim2;
    if ~isempty(ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.psi)
      psi_2d = ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.psi;
    else
      error('ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.psi not provided');
    end
  else
    error([' cannot find ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.grid_type.index == 1']);
  end
end
if ~isempty(ids_equilibrium.time_slice{i_equil}.boundary.outline.r)
  r_lcfs = ids_equilibrium.time_slice{i_equil}.boundary.outline.r;
else
  error(['ids_equilibrium.time_slice{i_equil}.boundary.outline.r not provided']);
end
if ~isempty(ids_equilibrium.time_slice{i_equil}.boundary.outline.z)
  z_lcfs = ids_equilibrium.time_slice{i_equil}.boundary.outline.z;
else
  error(['ids_equilibrium.time_slice{i_equil}.boundary.outline.z not provided']);
end

eqdsk_out.cocos = 11;
eqdsk_out.nr = numel(rmesh);
eqdsk_out.nz = numel(zmesh);
eqdsk_out.rmesh = rmesh;
eqdsk_out.zmesh = zmesh;
eqdsk_out.zshift = 0.;
eqdsk_out.psi = psi_2d;
eqdsk_out.psirz = psi_2d;
eqdsk_out.rboxlen = eqdsk_out.rmesh(end) - eqdsk_out.rmesh(1);
eqdsk_out.rboxleft=eqdsk_out.rmesh(1);
eqdsk_out.zboxlen = eqdsk_out.zmesh(end) - eqdsk_out.zmesh(1);
eqdsk_out.zmid = 0.5*(eqdsk_out.zmesh(end) + eqdsk_out.zmesh(1)) ;
eqdsk_out.psiaxis = psi_axis
eqdsk_out.psiedge = psi_boundary;
if ids_equilibrium.time_slice{i_equil}.boundary_separatrix.type == 1
  % diverted
  % could check psi value with psi boundary
end

% need psimesh ot be equidistant and same nb of points as nr
eqdsk_out.psimesh=linspace(0,1,eqdsk_out.nr)';
eqdsk_out.rhopsi = sqrt(eqdsk_out.psimesh);
% psimesh assumed from 0 on axis (so psi_edge is delta_psi) to have correct d/dpsi but to have sign(dpsi) easily
% (psimesh from 0 to 1) eqdsk_out.psimesh = (eqdsk_out.psiedge-eqdsk_out.psiaxis) .* eqdsk_out.psimesh;
nrho = size(pressure,1);
eqdsk_out.p = interpos(rhopol_norm,pressure,eqdsk_out.rhopsi,-0.01,[1 2],[0 pressure(end)]);
eqdsk_out.pprime = interpos(rhopol_norm,dpressure_dpsi,eqdsk_out.rhopsi,-0.01,[1 2],[0 dpressure_dpsi(end)]);
eqdsk_out.F = interpos(rhopol_norm,f_rbphi,eqdsk_out.rhopsi,-0.01,[1 2],[0 f_rbphi(end)]);
eqdsk_out.FFprime = interpos(rhopol_norm,f_df_dpsi,eqdsk_out.rhopsi,-0.01,[1 2],[0 f_df_dpsi(end)]);
eqdsk_out.q = interpos(rhopol_norm,q_profile,eqdsk_out.rhopsi,-0.01,[1 2],[0 q_profile(end)]);

eqdsk_out.stitle = [eqdsk_filename_suffix ' tref=' num2str(time_ref)];

eqdsk_out.nbbound = numel(ids_equilibrium.time_slice{i_equil}.boundary.outline.r);
if ~isempty(ids_equilibrium.time_slice{i_equil}.boundary.outline.r) && eqdsk_out.nbbound > 1 && ...
      eqdsk_out.nbbound == numel(ids_equilibrium.time_slice{i_equil}.boundary.outline.z)
  eqdsk_out.rplas = ids_equilibrium.time_slice{i_equil}.boundary.outline.r;
  eqdsk_out.zplas = ids_equilibrium.time_slice{i_equil}.boundary.outline.z;
else
  error('Plasma boundary not provided in time_slice{i_equil}.boundary.outline.r');
end

% extract limiter + double ids_wall if available
if numel(ids_wall.description_2d)>0 && numel(ids_wall.description_2d{1}.limiter.unit)>0 && ...
      ~isempty(ids_wall.description_2d{1}.limiter.unit{1}.outline.r) && ...
      numel(ids_wall.description_2d{1}.vessel.unit)>=2 && ...
      ~isempty(ids_wall.description_2d{1}.vessel.unit{1}.annular.centreline.r) && ...
      ~isempty(ids_wall.description_2d{1}.vessel.unit{2}.annular.centreline.r)
  Rlimiter = ids_wall.description_2d{1}.limiter.unit{1}.outline.r;
  Zlimiter = ids_wall.description_2d{1}.limiter.unit{1}.outline.z;
  if numel(ids_wall.description_2d{1}.limiter.unit)>=2 && ~isempty(ids_wall.description_2d{1}.limiter.unit{2}.outline.r)
    nbpts2 = numel(ids_wall.description_2d{1}.limiter.unit{2}.outline.r);
    Rlimiter(end+1:end+nbpts2) = ids_wall.description_2d{1}.limiter.unit{2}.outline.r(nbpts2:-1:1);
    Zlimiter(end+1:end+nbpts2) = ids_wall.description_2d{1}.limiter.unit{2}.outline.z(nbpts2:-1:1);
  end
  Rlimiter(end+1) = Rlimiter(1);
  Zlimiter(end+1) = Zlimiter(1);
  
  % vessel
  Rvessel_1 = ids_wall.description_2d{1}.vessel.unit{1}.annular.centreline.r;
  Zvessel_1 = ids_wall.description_2d{1}.vessel.unit{1}.annular.centreline.z;
  Rvessel_2 = ids_wall.description_2d{1}.vessel.unit{2}.annular.centreline.r;
  Zvessel_2 = ids_wall.description_2d{1}.vessel.unit{2}.annular.centreline.z;
  
  aa = [Rlimiter Zlimiter];
  nlim = size(aa,1);
  nvessel_1 = numel(Rvessel_1);
  nvessel_2 = numel(Rvessel_2);

  d_lim_vessel_1 = sqrt((Rvessel_1-Rlimiter(end)).^2 + (Zvessel_1-Zlimiter(end)).^2);
  [zzz,idmin] = min(d_lim_vessel_1);
  aa(nlim+1:nlim+(nvessel_1-idmin+1),1) = Rvessel_1(idmin:end);
  aa(nlim+(nvessel_1-idmin+1)+1:nlim+nvessel_1+1,1) = Rvessel_1(1:idmin); % end return point
  aa(nlim+1:nlim+(nvessel_1-idmin+1),2) = Zvessel_1(idmin:end);
  aa(nlim+(nvessel_1-idmin+1)+1:nlim+nvessel_1+1,2) = Zvessel_1(1:idmin);
  
  d_lim_vessel_2 = sqrt((Rvessel_2-Rlimiter(end)).^2 + (Zvessel_2-Zlimiter(end)).^2);
  [zzz,idmin] = min(d_lim_vessel_2);
  aa(nlim+nvessel_1+1+1:nlim+nvessel_1+1+(nvessel_2-idmin+1),1) = Rvessel_2(idmin:end);
  aa(nlim+nvessel_1+1+(nvessel_2-idmin+1)+1:nlim+nvessel_1+1+nvessel_2+1,1) = Rvessel_2(1:idmin);
  aa(nlim+nvessel_1+1+1:nlim+nvessel_1+1+(nvessel_2-idmin+1),2) = Zvessel_2(idmin:end);
  aa(nlim+nvessel_1+1+(nvessel_2-idmin+1)+1:nlim+nvessel_1+1+nvessel_2+1,2) = Zvessel_2(1:idmin);
  
else
  % read from previously defined boundary
  try
    aa = load('ITERlimiter.data');
  catch
    aa = [];
  end
end

if ~isempty(aa)
  eqdsk_out.nblim = size(aa,1);
  eqdsk_out.rlim = aa(:,1);
  eqdsk_out.zlim = aa(:,2);
else
  eqdsk_out.nblim = 0;
  eqdsk_out.rlim = [];
  eqdsk_out.zlim = [];
end

if write_opt == 1
  write_eqdsk(['EQDSK_COCOS' num2str(eqdsk_out.cocos) '_' eqdsk_filename_suffix],eqdsk_out,[eqdsk_out.cocos eqdsk_out.cocos]);
end
