function [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time,zshift,varargin);
%
% [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time,zshift,varargin);
%
% if shot is a structure assumes obtained from gdat(shot,'equil',...);
%
% varargin{1:2}: 'source','EQI' (default) or 'source','IDE'  (used in gdat(...,'equil') )
% varargin{3:4}: 'extra_arg_sf2sig','[]' (default) or 'extra_arg_sf2sig','''-ed'',2'
%
% [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time);
% [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(shot,time,[],'source','IDE');
%
% you can then do:
%     write_eqdsk('fname',eqdskAUG,17); % EQI is COCOS=17 by default
%     write_eqdsk('fname',eqdskAUG,[17 11]); % if you want an ITER version with COCOS=11
%


if ~exist('shot') || isempty(shot); 
  disp('requires a shot or equil structure from gdat(shot,''equil'',...) in geteqdskAUG')
  return
end

if ~exist('time'); 
  time_eff = 2.0;
else
  time_eff = time;
end

if ~exist('zshift') || isempty(zshift) || ~isnumeric(zshift)
  zshift_eff = 0.; % no shift
else
  zshift_eff = zshift;
end

if nargin >= 5 && ~isempty(varargin{1}) && strcmp(lower(varargin{1}),'source')
  if ~isempty(varargin{2})
    equil_source = varargin{2}
  else
    disp(['Warning source for equil in geteqdskAUG not defined']);
    return;
  end
else
  equil_source = 'EQI';
end
equilpar_source = 'FPG';
if strcmp(lower(equil_source),'ide');   equilpar_source = 'IDG'; end

if nargin >= 7 && ~isempty(varargin{3}) && strcmp(lower(varargin{3}),'extra_arg_sf2sig')
  if ~isempty(varargin{4})
    extra_arg_sf2sig = varargin{4}
  else
    disp(['Warning extra_arg_sf2sig in geteqdskAUG not defined']);
    return;
  end
else
  extra_arg_sf2sig = '[]';
end

if isnumeric(shot)
  equil=gdat(shot,'equil','equil',equil_source,'extra_arg_sf2sig',extra_arg_sf2sig);
else
  equil = shot;
  shot = equil.shot;
end
[zz it]=min(abs(equil.t-time_eff));

equil_all_t = equil;
equil_t_index = it;

eqdsk.cocos=17;
eqdsk.nr = size(equil.Rmesh,1);
eqdsk.nz = size(equil.Zmesh,1);
eqdsk.rmesh = equil.Rmesh(:,it);
eqdsk.zmesh = equil.Zmesh(:,it) - zshift_eff;
eqdsk.zshift = zshift_eff;
eqdsk.psi = equil.psi2D(:,:,it);
eqdsk.psirz = equil.psi2D(:,:,it);
eqdsk.rboxlen = equil.Rmesh(end,it) - equil.Rmesh(1,it) ;
eqdsk.rboxleft=eqdsk.rmesh(1);
eqdsk.zboxlen = equil.Zmesh(end,it) - equil.Zmesh(1,it) ;
eqdsk.zmid = 0.5*(equil.Zmesh(end,it) + equil.Zmesh(1,it)) ;
eqdsk.psiaxis = equil.psi_axis(it);
eqdsk.psiedge = equil.psi_lcfs(it);
% need psimesh ot be equidistant and same nb of points as nr
eqdsk.psimesh=linspace(0,1,eqdsk.nr)';
eqdsk.rhopsi = sqrt(eqdsk.psimesh);
% psimesh assumed from 0 on axis (so psi_edge is delta_psi) to have correct d/dpsi but to have sign(dpsi) easily
eqdsk.psimesh = (eqdsk.psiedge-eqdsk.psiaxis) .* eqdsk.psimesh;
nrho = size(equil.pressure,1);
eqdsk.p = interpos(equil.rhopolnorm(:,it),equil.pressure(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.pressure(end,it)]);
eqdsk.pprime = interpos(equil.rhopolnorm(:,it),equil.dpressuredpsi(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.dpressuredpsi(end,it)]);
eqdsk.FFprime = interpos(equil.rhopolnorm(:,it),equil.ffprime(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.ffprime(end,it)]);
eqdsk.q = interpos(equil.rhopolnorm(:,it),equil.qvalue(:,it),eqdsk.rhopsi,-0.03,[1 2],[0 equil.qvalue(end,it)]);

eqdsk.ip = equil.Ip(it);
eqdsk.stitle=['#' num2str(shot) ' t=' num2str(time_eff) ' from ' equil.gdat_params.exp_name '/' equil.gdat_params.equil];
eqdsk.ind1=1;

psisign = sign(eqdsk.psimesh(end)-eqdsk.psimesh(1));
[dum1,dum2,dum3,F2_05]=interpos(psisign.*eqdsk.psimesh,eqdsk.FFprime,-0.1);

b0=gdat(shot,'b0');
[zz itb0]=min(abs(b0.t-time_eff));
eqdsk.b0 = b0.data(itb0);
eqdsk.r0 = 1.65;
fedge=eqdsk.r0.*eqdsk.b0;
F2 = psisign.*2.*F2_05 + fedge.^2;
eqdsk.F = sqrt(F2)*sign(eqdsk.b0);

rmag=gdat(shot,'rmag','source',equilpar_source,'extra_arg_sf2sig',extra_arg_sf2sig);
[zz itrmag]=min(abs(rmag.t-time_eff));
eqdsk.raxis = rmag.data(itrmag);
zmag=gdat(shot,'zmag','source',equilpar_source,'extra_arg_sf2sig',extra_arg_sf2sig);
eqdsk.zaxis = zmag.data(itrmag) - eqdsk.zshift;

% get plasma boundary
figure
contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',100)
hold on
psiedge_eff = 0.001*eqdsk.psiaxis + 0.999*eqdsk.psiedge;
[hh1 hh2]=contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',[psiedge_eff psiedge_eff],'k');
axis equal
if ~isempty(hh1)
  ij=1;
  ij_prev = 0;
  nbhh1(ij)=hh1(2,ij_prev+1);
  Rbnd{ij}=hh1(1,2:1+nbhh1(ij));
  Zbnd{ij}=hh1(2,2:1+nbhh1(ij));
  ij_prev = ij_prev+1+nbhh1(ij);
  while (ij_prev+1<size(hh1,2))
    % next
    ij=ij + 1;
    nbhh1(ij)=hh1(2,ij_prev+1);
    Rbnd{ij}=hh1(1,ij_prev+2:ij_prev+1+nbhh1(ij));
    Zbnd{ij}=hh1(2,ij_prev+2:ij_prev+1+nbhh1(ij));
    ij_prev = ij_prev+1+nbhh1(ij);
  end
  % assume LCFS with most points
  [zzz irz]=max(nbhh1);
  eqdsk.nbbound = nbhh1(irz);
  eqdsk.rplas = Rbnd{irz}';
  eqdsk.zplas = Zbnd{irz}';
  plot(eqdsk.rplas,eqdsk.zplas,'k-')
else
  eqdsk.nbbound = [];
  eqdsk.rplas = [];
  eqdsk.zplas = [];
end

[aget]=which('geteqdskAUG');
[path1,name2,ext3]=fileparts(aget);
eval(['load ' fullfile(path1,'AUG_innerouterwall.data')])
eqdsk.nblim=size(AUG_innerouterwall,1);
eqdsk.rlim=AUG_innerouterwall(:,1);
eqdsk.zlim=AUG_innerouterwall(:,2);
plot(eqdsk.rlim,eqdsk.zlim,'k-')

title(eqdsk.stitle)

eqdskAUG = eqdsk;
