geteqdskAUG.m

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'
% varargin{5:6}: 'fignb',ii (default is a new fig each time, if a number is given it over plots on this figure, needed if called in series)
%                        if <=0: no plots
%
% [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(upper(equil_source),'IDE');   equilpar_source = 'IDG'; end
if strcmp(upper(equil_source),'EQI');   equilpar_source = 'GQI'; end
if strcmp(upper(equil_source),'EQH');   equilpar_source = 'GQH'; end
if strcmp(upper(equil_source),'TRA');   equilpar_source = 'TRE'; 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 nargin >= 9 && ~isempty(varargin{5}) && strcmp(lower(varargin{5}),'fignb')
  if ~isempty(varargin{6})
    fignb = varargin{6};
  else
    disp(['Warning fignb in geteqdskAUG not defined']);
    return;
  end
else
  fignb = [];
end

if isnumeric(shot)
  equil=gdat_aug(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);
ixpoint = 2;
if abs(equil.Xpoints.psi(ixpoint,it)-eqdsk.psiedge) < 0.1.*abs(eqdsk.psiedge-eqdsk.psiaxis)
  eqdsk.psiedge = equil.Xpoints.psi(ixpoint,it);
end
% 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
% (psimesh from 0 to 1) 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.01,[1 2],[0 equil.pressure(end,it)]);
eqdsk.pprime = interpos(equil.rhopolnorm(:,it),equil.dpressuredpsi(:,it),eqdsk.rhopsi,-0.01,[1 2],[0 equil.dpressuredpsi(end,it)]);
eqdsk.FFprime = interpos(equil.rhopolnorm(:,it),equil.ffprime(:,it),eqdsk.rhopsi,-0.01,[1 2],[0 equil.ffprime(end,it)]);
if abs(equil.qvalue(end,it)) > 25
  eqdsk.q = interpos(21,equil.rhopolnorm(:,it),equil.qvalue(:,it),eqdsk.rhopsi);
else
  eqdsk.q = interpos(equil.rhopolnorm(:,it),equil.qvalue(:,it),eqdsk.rhopsi,-0.01,[1 2],[0 equil.qvalue(end,it)]);
end
if strcmp(upper(equil_source),'IDE');
  % q of IDE ids abs(q), add sign from EQI
  q95_eqi=gdat_aug(shot,{'GQI','q95'},'extra_arg_sf2sig',extra_arg_sf2sig);
  eqdsk.q = abs(eqdsk.q) * sign(q95_eqi.data(round(length(q95_eqi.data)/2)));
end
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.01);

b0=gdat_aug(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_aug(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_aug(shot,'zmag','source',equilpar_source,'extra_arg_sf2sig',extra_arg_sf2sig);
eqdsk.zaxis = zmag.data(itrmag) - eqdsk.zshift;

% get plasma boundary
doplot = 1;
if isnumeric(fignb) && fignb <= 0
  doplot = 0;
elseif isempty(fignb) || (isnumeric(fignb) && fignb < 1) || (~isnumeric(fignb) && ~ishandle(fignb)) 
  figure
  fignb_handle = gcf;
  set(gcf,'Name','from geteqdskAUG');
else
  figure(fignb)
  fignb_handle = gcf;
  clf
  set(gcf,'Name','from geteqdskAUG');
end

if doplot
  contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',100)
  hold on
end
psiedge_eff = 0.001*eqdsk.psiaxis + 0.999*eqdsk.psiedge; % note for IDE psi edge is already 99%
%ixpoint=iround_os(equil.Xpoints.psi(:,it),eqdsk.psiedge);
if abs(equil.Xpoints.psi(ixpoint,it)-eqdsk.psiedge) < 0.1.*abs(eqdsk.psiedge-eqdsk.psiaxis)
  psiedge_eff = 0.002*eqdsk.psiaxis + 0.998*equil.Xpoints.psi(ixpoint,it);
end
if doplot
  [hh1 hh2]=contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',[psiedge_eff psiedge_eff],'k');
  axis equal
else
  [hh1]=contourc(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',[psiedge_eff psiedge_eff]);
end
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}';
  if doplot
    plot(eqdsk.rplas,eqdsk.zplas,'k-')
  end
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);
if doplot
  plot(eqdsk.rlim,eqdsk.zlim,'k-')
  title(eqdsk.stitle)
end
eqdskAUG = eqdsk;