From 02fa09416b98a2998152800e74922a6290d0e099 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Mon, 12 Feb 2018 16:49:30 +0000
Subject: [PATCH] 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
---
crpptbx/AUG/aug_requests_mapping.m | 1 +
crpptbx/AUG/geteqdskAUG.m | 20 ++++++++++++++------
2 files changed, 15 insertions(+), 6 deletions(-)
diff --git a/crpptbx/AUG/aug_requests_mapping.m b/crpptbx/AUG/aug_requests_mapping.m
index 92a991d3..c234dae4 100644
--- a/crpptbx/AUG/aug_requests_mapping.m
+++ b/crpptbx/AUG/aug_requests_mapping.m
@@ -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';
diff --git a/crpptbx/AUG/geteqdskAUG.m b/crpptbx/AUG/geteqdskAUG.m
index 928de02a..3ca27016 100644
--- a/crpptbx/AUG/geteqdskAUG.m
+++ b/crpptbx/AUG/geteqdskAUG.m
@@ -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)
--
GitLab