Skip to content
Snippets Groups Projects
Commit 02fa0941 authored by Olivier Sauter's avatar Olivier Sauter
Browse files

use X-point psi to find plasma boundary for eqdsk, might need to check when...

use X-point psi to find plasma boundary for eqdsk, might need to check when limited. Use FPC for deafault Ip, closer to EQI/IDE

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@9334 d63d8f72-b253-0410-a779-e742ad2e26cf
parent e7a528d8
No related branches found
No related tags found
No related merge requests found
......@@ -149,6 +149,7 @@ switch lower(data_request)
mapping.label = 'Plasma current';
mapping.method = 'signal';
mapping.expression = [{'MAG'},{'Ipa'}];
mapping.expression = [{'FPC'},{'IpiFP'}];
case 'kappa'
mapping.timedim = 1;
mapping.label = '\kappa';
......
......@@ -82,19 +82,23 @@ 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
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.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.03,[1 2],[0 equil.qvalue(end,it)]);
eqdsk.q = interpos(equil.rhopolnorm(:,it),equil.qvalue(:,it),eqdsk.rhopsi,-0.01,[1 2],[0 equil.qvalue(end,it)]);
end
eqdsk.ip = equil.Ip(it);
......@@ -102,7 +106,7 @@ eqdsk.stitle=['#' num2str(shot) ' t=' num2str(time_eff) ' from ' equil.gdat_para
eqdsk.ind1=1;
psisign = sign(eqdsk.psimesh(end)-eqdsk.psimesh(1));
[dum1,dum2,dum3,F2_05]=interpos(psisign.*eqdsk.psimesh,eqdsk.FFprime,-0.1);
[dum1,dum2,dum3,F2_05]=interpos(psisign.*eqdsk.psimesh,eqdsk.FFprime,-0.01);
b0=gdat(shot,'b0');
[zz itb0]=min(abs(b0.t-time_eff));
......@@ -122,7 +126,11 @@ eqdsk.zaxis = zmag.data(itrmag) - eqdsk.zshift;
figure
contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',100)
hold on
psiedge_eff = 0.001*eqdsk.psiaxis + 0.999*eqdsk.psiedge;
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
[hh1 hh2]=contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',[psiedge_eff psiedge_eff],'k');
axis equal
if ~isempty(hh1)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment