-
Olivier Sauter authored
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@9099 d63d8f72-b253-0410-a779-e742ad2e26cf
Olivier Sauter authoredgit-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@9099 d63d8f72-b253-0410-a779-e742ad2e26cf
geteqdskAUG.m 5.17 KiB
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;