From bc16fc47c53fdb7326257bea508700f82b44485c Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Tue, 4 Jul 2017 13:05:44 +0000
Subject: [PATCH] added several new signals and tested geteqdskJET.m. This file
will be moved to CHEASEgui repository since it uses write_eqdsk and should be
in main paths
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7751 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/JET/gdat_jet.m | 39 +--
crpptbx/JET/get_grids_1d.m | 25 +-
crpptbx/JET/geteqdskJET.m | 537 +++++++++++++++++------------
crpptbx/JET/jet_requests_mapping.m | 28 +-
4 files changed, 381 insertions(+), 248 deletions(-)
diff --git a/crpptbx/JET/gdat_jet.m b/crpptbx/JET/gdat_jet.m
index 1bd04176..89c06ab8 100644
--- a/crpptbx/JET/gdat_jet.m
+++ b/crpptbx/JET/gdat_jet.m
@@ -244,7 +244,6 @@ gdat_data.gdat_params.help = jet_help_parameters(fieldnames(gdat_data.gdat_param
gdat_data.mapping_for.jet = mapping_for_jet;
gdat_params = gdat_data.gdat_params;
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1st treat the simplest method: "signal"
if strcmp(mapping_for_jet.method,'signal')
@@ -1193,19 +1192,20 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
phi = a;
phi_edge_2d = ones(length(x),1)*a(end,:);
gdat_data.rhotornorm = sqrt(phi./phi_edge_2d);
+ gdat_data.phi_tor = phi;
+ params_eff = gdat_data.gdat_params;
params_eff.data_request='b0';
b0=gdat_jet(shot,params_eff);
- if strcmp(data_request,'rhotor_edge')
- b0_phi=interpos(b0.t,b0.data,t,-1);
- gdat_data.data = sqrt(phi(end,:)./pi./reshape(b0_phi,1,length(b0_phi)));
+ gdat_data.b0 = interpos(b0.t,b0.data,t,-0.1);
+ if strcmp(data_request_eff,'rhotor_edge')
+ gdat_data.data = sqrt(phi(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)));
gdat_data.dim = {t};
gdat_data.t = gdat_data.dim{1};
gdat_data.data_fullpath='sqrt(Phi_edge/pi/B0) from {''ppf''},{''EFIT''},{''FTOR''}';
gdat_data.units = 'm';
gdat_data.dimunits{1} = 's';
- elseif strcmp(data_request,'rhotor')
- b0_phi=interpos(b0.t,b0.data,t,-1);
- gdat_data.data = sqrt(phi./(ones(length(x),1)*reshape(b0_phi,1,length(b0_phi)))./pi);
+ elseif strcmp(data_request_eff,'rhotor')
+ gdat_data.data = sqrt(phi./(ones(length(x),1)*reshape(gdat_data.b0,1,numel(gdat_data.b0)))./pi);
gdat_data.dim{1} = x;
gdat_data.dim{2} = t;
gdat_data.x = x;
@@ -1214,29 +1214,27 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.units = 'm';
gdat_data.dimunits{1} = 'rhopol\_norm';
gdat_data.dimunits{2} = 's';
- elseif strcmp(data_request,'rhotor_norm')
- b0_phi=interpos(b0.t,b0.data,t,-1);
+ elseif strcmp(data_request_eff,'rhotor_norm')
gdat_data.data = gdat_data.rhotornorm;
- gdat_data.rhotor_edge = sqrt(phi(end,:)./pi./reshape(b0_phi,1,length(b0_phi)));
+ gdat_data.rhotor_edge = sqrt(phi(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)));
gdat_data.dim{1} = x;
gdat_data.dim{2} = t;
gdat_data.x = x;
gdat_data.t = gdat_data.dim{2};
gdat_data.data_fullpath='sqrt(phitor/phitor_edge), rhotor_edge=sqrt(phitor/B0/pi) from {''ppf''},{''EFIT''},{''FTOR''}';
- gdat_data.data_fullpath='sqrt(phitor/pi/B0), rhotor_edge=sqrt(phitor/B0/pi)';
gdat_data.units = '';
gdat_data.dimunits{1} = 'rhopol\_norm';
gdat_data.dimunits{2} = 's';
end
case {'volume_rho', 'rhovol'}
- nodenameeff=[{'ppf'},{'EFIT'},{'VJAC'}];
+ nodenameeff=[{'ppf'},{'EFIT'},{'VJAC'}]; % dVdpsi?
ppftype = nodenameeff{1};
tracename = [nodenameeff{2} '/' nodenameeff{3}];
[a,x,t,d,e]=rda_jet(shot,ppftype,tracename);
nodenameeff2=[{'ppf'},{'EFIT'},{'VOLM'}];
- ppftype = nodenameeff{1};
- tracename = [nodenameeff{2} '/' nodenameeff{3}];
- [a2,x2,t2,d2,e2]=rda_jet(shot,ppftype,tracename);
+ ppftype2 = nodenameeff2{1};
+ tracename2 = [nodenameeff2{2} '/' nodenameeff2{3}];
+ [a2,x2,t2,d2,e2]=rda_jet(shot,ppftype2,tracename2);
if e==0 || e2==0
if gdat_params.nverbose>=3;
disp(['error after rda_jet in signal with nodenameeff= ' nodenameeff]);
@@ -1254,13 +1252,14 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.x = sqrt(x);
gdat_data.t = t;
gdat_data.dim = [{gdat_data.x}, {gdat_data.t}];
+ gdat_data.dimunits = [{'rhopol'}, {'s'}];
switch data_request_eff
case 'volume_rho'
gdat_data.data = volpsi;
- gdat_data.volume = volume;
+ gdat_data.volume_edge = volpsi(end,:);
case 'rhovol'
gdat_data.data = sqrt(volpsi./(ones(length(gdat_data.x),1)*reshape(volpsi(end,:),1,length(gdat_data.t))));
- gdat_data.volume = volpsi(end,:);
+ gdat_data.volume_edge = volpsi(end,:);
end
otherwise
disp(['data_request_eff = ' data_request_eff ' not defined in this section']);
@@ -1288,7 +1287,7 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
end
end
gdat_data.units = tracetdi.units;
- if strcmp(data_request,'volume')
+ if strcmp(data_request_eff,'volume')
gdat_data.data = tracetdi.data(end,:);
gdat_data.dim{1} = tracetdi.dim{2};
gdat_data.data_fullpath=['\results::psitbx:vol(end,:)'];
@@ -1298,10 +1297,10 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
gdat_data.data = tracetdi.data;
gdat_data.dim = tracetdi.dim;
gdat_data.dimunits = tracetdi.dimunits;
- if strcmp(data_request,'volume_rho')
+ if strcmp(data_request_eff,'volume_rho')
gdat_data.data_fullpath=['\results::psitbx:vol'];
gdat_data.request_description = 'volume(rho)';
- elseif strcmp(data_request,'rhovol')
+ elseif strcmp(data_request_eff,'rhovol')
gdat_data.volume_edge = gdat_data.data(end,:);
gdat_data.data = sqrt(gdat_data.data./repmat(reshape(gdat_data.volume_edge,1,size(gdat_data.data,2)),size(gdat_data.data,1),1));
gdat_data.data_fullpath='sqrt(\results::psitbx:vol/vol_edge)';
diff --git a/crpptbx/JET/get_grids_1d.m b/crpptbx/JET/get_grids_1d.m
index b5cb83b5..78343647 100644
--- a/crpptbx/JET/get_grids_1d.m
+++ b/crpptbx/JET/get_grids_1d.m
@@ -26,12 +26,14 @@ if (nopt == 0) || isempty(gdat_data.x) || isempty(gdat_data.t) || isempty(gdat_d
end
params_eff = gdat_data.gdat_params;params_eff.doplot=0;
-params_eff.data_request='rhotor';
-rhotor = gdat(gdat_data.shot,params_eff);
+params_eff.data_request='rhotor_norm';
+rhotor_norm = gdat(gdat_data.shot,params_eff)
params_eff.data_request='rhovol';
rhovol = gdat(gdat_data.shot,params_eff);
+params_eff.data_request='psi_axis';
+psi_axis = gdat(gdat_data.shot,params_eff);
+psi0 = interpos(psi_axis.t,psi_axis.data,gdat_data.t,-0.01);
-psi0 = interpos(rhotor.t',rhotor.psi_axis,gdat_data.t,-0.01);
if (nbdim_x == 1)
gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0,1,length(psi0));
elseif (nbdim_x == 2)
@@ -42,7 +44,8 @@ else
end
gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data));
gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data));
-it_rt = iround_os(rhotor.t,gdat_data.t);
+
+it_rt = iround_os(rhotor_norm.t,gdat_data.t);
it_vol = iround_os(rhovol.t,gdat_data.t);
for it=1:length(gdat_data.t)
% do an interpolation on closest point to avoid 2D interp
@@ -55,17 +58,21 @@ for it=1:length(gdat_data.t)
end
if (nbdim_x == 1)
if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
- gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor.x,rhotor.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
- gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm);
+ gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
+ if ~isempty(rhovol.x)
+ gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm);
+ end
end
else
if length(ii)==size(gdat_data.grids_1d.rhopolnorm,1)
- gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor.x,rhotor.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
- gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+ gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+ if ~isempty(rhovol.x)
+ gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+ end
end
end
end
-gdat_data.grids_1d.rhotor_edge=interpos(rhotor.t',rhotor.rhotor_edge,gdat_data.t',-0.01);
+gdat_data.grids_1d.rhotor_edge=interpos(rhotor_norm.t',rhotor_norm.rhotor_edge,gdat_data.t',-0.01);
gdat_data.grids_1d.volume_edge=interpos(rhovol.t',rhovol.volume_edge,gdat_data.t',-0.01);
diff --git a/crpptbx/JET/geteqdskJET.m b/crpptbx/JET/geteqdskJET.m
index fe024ccf..f4ed9d20 100644
--- a/crpptbx/JET/geteqdskJET.m
+++ b/crpptbx/JET/geteqdskJET.m
@@ -1,22 +1,22 @@
-function [efitdata]=geteqdskJET(shot,time,nrg,nzg,savedir,deltaz,efitlab,uid,seqd,varargin);
+function [efitdata,eqd]=geteqdskJET(shot,times_in,nrg,nzg,savedir,deltaz,efitlab,uid,seqd,varargin);
%
-% function [efitdata]=geteqdskJET(shot,time,nrg,nzg,savedir,deltaz,efitlab,uid,seqd,varargin);
+% function [efitdata]=geteqdskJET(shot,times_in,nrg,nzg,savedir,deltaz,efitlab,uid,seqd,varargin);
%
-% save eqdsk file to savedir/JET_EQDSK_shot_ttime
+% save eqdsk file to savedir/JET_EQDSK_shot_t(times_in(:))
% returns data used in structure efitdata (efitdata.shot, efitdata.Rbnd, efitdataZbnd....)
%
% examples:
-% efitdata=geteqdskJET(shot,time,33,65,'/tmp')
-% efitdata=geteqdskJET(shot,time,33,65,'/tmp',[],[],[],[],1,efitdata); % gives data so no need to load again
+% efitdata=geteqdskJET(shot,times_in,33,65,'/tmp')
+% efitdata=geteqdskJET(shot,times_in,33,65,'/tmp',[],[],[],[],1,efitdata); % gives data so no need to load again
%
% INPUTS :
% shot : shot number
-% time : single time slice of analysis
-% nrg, nzg: nb of R and Z points for R,Z grid
-% if nzg<0 force symmetric box aorund z=deltaz
+% times_in : single time slice of analysis
%
% OPTIONAL
-% savedir: directory to save eqdsk file (default: './')
+% nrg, nzg: nb of R and Z points for R,Z grid: default: 129, -129
+% if nzg<0 force symmetric box aorund z=deltaz
+% savedir: directory to save eqdsk file (default: '/tmp/username')
% deltaz: shift equilibrium vertically
% = 0: (default) no shift
% = -99: shift so that zmag=0 (thus deltaz=-zmag)
@@ -34,8 +34,17 @@ function [efitdata]=geteqdskJET(shot,time,nrg,nzg,savedir,deltaz,efitlab,uid,seq
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+if ~exist('nrg') || isempty(nrg)
+ nrg = 129;
+end
+if ~exist('nzg') || isempty(nzg)
+ nzg = -129;
+end
if nargin<=4 | isempty(savedir)
- savedir='./';
+ savedir=['/tmp/' getenv('USER')];
+ if ~exist(savedir,'dir')
+ unix(['mkdir ' savedir]);
+ end
end
if unix(['test -d ' savedir])
disp(['Problems in geteqdskJET, savedir=' savedir ' is not a directory'])
@@ -88,224 +97,316 @@ if nargin>=11 & ~isempty(varargin{2})
end
end
+do_plot=0;
+exp_machine='jet';
if iread==1
efitdata.shot=shot;
% get data needed
- efitdata.Rbnd=gdat(shot,['ppf/' efitlab '/RBND' s_extra{iexefit}],0,'JET');
+ efitdata.Rbnd=gdat(shot,[{'ppf'}, {efitlab}, {['RBND' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
efitdata.tefit=efitdata.Rbnd.t;
- efitdata.Zbnd=gdat(shot,['ppf/' efitlab '/ZBND' s_extra{iexefit}],0,'JET');
- efitdata.R0=gdat(shot,['ppf/' efitlab '/RGEO' s_extra{iexefit}],0,'JET');
- efitdata.a=gdat(shot,['ppf/' efitlab '/CR0' s_extra{iexefit}],0,'JET');
- % $$$ efitdata.Z0=gdat(shot,['ppf/' efitlab '/ZO' s_extra{iexefit}],0,'JET');
- efitdata.rmag=gdat(shot,['ppf/' efitlab '/rmag' s_extra{iexefit}],0,'JET');
- efitdata.zmag=gdat(shot,['ppf/' efitlab '/zmag' s_extra{iexefit}],0,'JET');
- efitdata.faxs=gdat(shot,['ppf/' efitlab '/faxs' s_extra{iexefit}],0,'JET');
- efitdata.fbnd=gdat(shot,['ppf/' efitlab '/fbnd' s_extra{iexefit}],0,'JET');
- efitdata.bvac=gdat(shot,['ppf/' efitlab '/bvac' s_extra{iexefit}],0,'JET');
- efitdata.ip=gdat(shot,['ppf/' efitlab '/xip' s_extra{iexefit}],0,'JET');
- efitdata.F=gdat(shot,['ppf/' efitlab '/F' s_extra{iexefit}],0,'JET');
+ efitdata.Zbnd=gdat(shot,[{'ppf'}, {efitlab}, {['ZBND' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.R0=gdat(shot,[{'ppf'}, {efitlab}, {['RGEO' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.a=gdat(shot,[{'ppf'}, {efitlab}, {['CR0' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.Z0=gdat(shot,[{'ppf'}, {efitlab}, {['ZGEO' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.rmag=gdat(shot,[{'ppf'}, {efitlab}, {['rmag' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.zmag=gdat(shot,[{'ppf'}, {efitlab}, {['zmag' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.faxs=gdat(shot,[{'ppf'}, {efitlab}, {['faxs' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.fbnd=gdat(shot,[{'ppf'}, {efitlab}, {['fbnd' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.bvac=gdat(shot,[{'ppf'}, {efitlab}, {['bvac' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.ip=gdat(shot,[{'ppf'}, {efitlab}, {['xip' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.F=gdat(shot,[{'ppf'}, {efitlab}, {['F' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
efitdata.psin=efitdata.F.x;
- efitdata.P=gdat(shot,['ppf/' efitlab '/P' s_extra{iexefit}],0,'JET');
- efitdata.Q=gdat(shot,['ppf/' efitlab '/Q' s_extra{iexefit}],0,'JET');
- efitdata.kappa=gdat(shot,['ppf/' efitlab '/ELON' s_extra{iexefit}],0,'JET');
- efitdata.q95=gdat(shot,['ppf/' efitlab '/q95' s_extra{iexefit}],0,'JET');
- efitdata.btpd=gdat(shot,['ppf/' efitlab '/btpd' s_extra{iexefit}],0,'JET');
- efitdata.bttd=gdat(shot,['ppf/' efitlab '/bttd' s_extra{iexefit}],0,'JET');
- efitdata.btnd=gdat(shot,['ppf/' efitlab '/btnd' s_extra{iexefit}],0,'JET');
- efitdata.btpm=gdat(shot,['ppf/' efitlab '/btpm' s_extra{iexefit}],0,'JET');
- efitdata.bttm=gdat(shot,['ppf/' efitlab '/bttm' s_extra{iexefit}],0,'JET');
- efitdata.btnm=gdat(shot,['ppf/' efitlab '/btnm' s_extra{iexefit}],0,'JET');
- efitdata.xli=gdat(shot,['ppf/' efitlab '/xli' s_extra{iexefit}],0,'JET');
- efitdata.sspr=gdat(shot,['ppf/' efitlab '/sspr' s_extra{iexefit}],0,'JET');
- efitdata.sspi=gdat(shot,['ppf/' efitlab '/sspi' s_extra{iexefit}],0,'JET');
- % add for profiles
- efitdata.ti=gdat(shot,['ppf/TION/TI' s_extra{iexchain2}],0,'JET');
- efitdata.p_tion=gdat(shot,['ppf/TION/p' s_extra{iexchain2}],0,'JET');
- efitdata.pi=gdat(shot,['ppf/NION/DD' s_extra{iexchain2}],0,'JET');
- efitdata.zef=gdat(shot,['ppf/NION/ZEF' s_extra{iexchain2}],0,'JET');
- % add for calculating NTM parameters
- efitdata.bpol=gdat(shot,['ppf/equi/bpol' s_extra{iexequi}],0,'JET');
- efitdata.bpo2=gdat(shot,['ppf/equi/bpo2' s_extra{iexequi}],0,'JET');
- efitdata.qmag=gdat(shot,['ppf/' efitlab '/qmag' s_extra{iexefit}],0,'JET');
- efitdata.lidrpe=gdat(shot,['ppf/lidr/pe'],0,'JET');
- efitdata.nexav=gdat(shot,['ppf/nex/av'],0,'JET');
- efitdata.nbi=gdat(shot,['ppf/nbi/ptot'],0,'JET');
- efitdata.icrh=gdat(shot,['ppf/icrh/ptot'],0,'JET');
- efitdata.ptot=gdat(shot,['ppf/mg3/yto'],0,'JET');
- efitdata.halpha=gdat(shot,['jpf/dd/s3-ad35'],0,'JET');
- efitdata.n1=gdat(shot,['jpf/da/c1-g101'],0,'JET');
- efitdata.n2=gdat(shot,['jpf/da/c1-g102'],0,'JET');
+ efitdata.P=gdat(shot,[{'ppf'}, {efitlab}, {['P' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.Q=gdat(shot,[{'ppf'}, {efitlab}, {['Q' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.kappa=gdat(shot,[{'ppf'}, {efitlab}, {['ELON' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.q95=gdat(shot,[{'ppf'}, {efitlab}, {['q95' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.btpd=gdat(shot,[{'ppf'}, {efitlab}, {['btpd' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.bttd=gdat(shot,[{'ppf'}, {efitlab}, {['bttd' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.btnd=gdat(shot,[{'ppf'}, {efitlab}, {['btnd' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.btpm=gdat(shot,[{'ppf'}, {efitlab}, {['btpm' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.bttm=gdat(shot,[{'ppf'}, {efitlab}, {['bttm' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.btnm=gdat(shot,[{'ppf'}, {efitlab}, {['btnm' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.xli=gdat(shot,[{'ppf'}, {efitlab}, {['li3m' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.sspr=gdat(shot,[{'ppf'}, {efitlab}, {['sspr' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+ efitdata.sspi=gdat(shot,[{'ppf'}, {efitlab}, {['sspi' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ % add for profiles
+% $$$ efitdata.ti=gdat(shot,[{'ppf'},{'TION'},{['TI' s_extra{iexchain2}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.p_tion=gdat(shot,[{'ppf'},{'TION'},{['p' s_extra{iexchain2}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.pi=gdat(shot,[{'ppf'},{'NION'},{['DD' s_extra{iexchain2}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.zef=gdat(shot,[{'ppf'},{'NION'},{['ZEF' s_extra{iexchain2}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ % add for calculating NTM parameters
+% $$$ efitdata.bpol=gdat(shot,[{'ppf'},{'equi'},{['bpol' s_extra{iexequi}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.bpo2=gdat(shot,[{'ppf'},{'equi'},{['bpo2' s_extra{iexequi}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.qmag=gdat(shot,[{'ppf'},{'' efitlab ''},{['qmag' s_extra{iexefit}]}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.lidrpe=gdat(shot,[{'ppf'},{'lidr'},{['pe']}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.nexav=gdat(shot,[{'ppf'},{'nex'},{['av']}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.nbi=gdat(shot,[{'ppf'},{'nbi'},{['ptot']}],'machine',exp_machine,'doplot',do_plot);
+% $$$ efitdata.icrh=gdat(shot,[{'ppf'},{'icrh'},{['ptot']}],'machine',exp_machine,'doplot',do_plot);
+% $$$ % old, not working anymore as of June 2017?
+% $$$ % $$$ efitdata.ptot=gdat(shot,[{'ppf'},{'mg3'},{['yto']}],'machine',exp_machine,'doplot',do_plot);
+% $$$ % $$$ efitdata.halpha=gdat(shot,[{'jpf'},{'dd'},{['s3-ad35']}],'machine',exp_machine,'doplot',do_plot);
+% $$$ % $$$ efitdata.n1=gdat(shot,[{'jpf'},{'da'},{['c1-g101']}],'machine',exp_machine,'doplot',do_plot);
+% $$$ % $$$ efitdata.n2=gdat(shot,[{'jpf'},{'da'},{['c1-g102']}],'machine',exp_machine,'doplot',do_plot);
end
-
tefit=efitdata.tefit;
-[zz index_efit]=min(abs(tefit-time));
-time_efit=tefit(index_efit);
-disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
-disp(['efit at t=' num2str(time_efit)])
-disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
-
-Rbnd=efitdata.Rbnd.data(:,index_efit);
-Zbnd=efitdata.Zbnd.data(:,index_efit);
-R0=efitdata.R0.data(index_efit);
-% $$$ Z0=efitdata.Z0.data(index_efit);
-rmag=efitdata.rmag.data(index_efit);
-zmag=efitdata.zmag.data(index_efit);
-if deltaz==-99
- deltaz=-zmag;
-end
-faxs=efitdata.faxs.data(index_efit);
-fbnd=efitdata.fbnd.data(index_efit);
-bvac=-efitdata.bvac.data(index_efit);
-ip=-efitdata.ip.data(index_efit);
-psi_efit=efitdata.psin;
-F=efitdata.F.data(:,index_efit);
-P=efitdata.P.data(:,index_efit);
-Q=efitdata.Q.data(:,index_efit);
-kappa=efitdata.kappa.data(index_efit);
-q95=efitdata.q95.data(index_efit);
-if index_efit<=length(efitdata.btpd.data); btpd=efitdata.btpd.data(index_efit); else; btpd=0; end
-if index_efit<=length(efitdata.bttd.data); bttd=efitdata.bttd.data(index_efit); else; bttd=0; end
-if index_efit<=length(efitdata.btnd.data); btnd=efitdata.btnd.data(index_efit); else; btnd=0; end
-if index_efit<=length(efitdata.btpm.data); btpm=efitdata.btpm.data(index_efit); else; btpm=0; end
-if index_efit<=length(efitdata.bttm.data); bttm=efitdata.bttm.data(index_efit); else; bttm=0; end
-if index_efit<=length(efitdata.btnm.data); btnm=efitdata.btnm.data(index_efit); else; btnm=0; end
-xli=efitdata.xli.data(index_efit);
-
-[R,Z,psinrz]=psirz(shot,time,nrg,nzg,efitlab,uid{iexefit},seqd(iexefit),ncont,efitdata.sspr,efitdata.sspi,deltaz);
-nzg=abs(nzg);
-
-% define file name
-s = sprintf('%.6f',time);
-fname=sprintf('%s/EQDSK.%st%s',savedir,num2str(shot),s);
-fid=fopen(fname,'w');
-
-% 1st eqdsk line: 48 characters for file description and then, 3, nr, nz
-tdate=date;
-ss=['JET #' num2str(shot) ' t= ' num2str(time_efit) ' ' efitlab '?uid=' uid{iexefit} '+seq=' num2str(seqd(iexefit)) ', ' tdate];
-if length(ss)<48
- ss(end:48)=' ';
-else
- ss=ss(1:48);
-end
-
-fprintf(fid,'%s%4d%4d%4d\n',ss,3,nrg,nzg);
-
-% 2nd line: rboxlen, zboxlen, r0, rboxlft, zboxmid
-rmin=min(R)
-rmax=max(R)
-zmin=min(Z)
-zmax=max(Z)
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',rmax-rmin,zmax-zmin,R0,rmin,0.5*(zmin+zmax)+deltaz);
+index_efit = iround_os(tefit,times_in);
+for it=1:numel(index_efit)
+ index_efit_eff=index_efit(it);
+ time_efit=tefit(index_efit_eff);
-% 3rd line: rmag, zmag, psimag, psiedge, B0
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',rmag,zmag+deltaz,faxs,fbnd,bvac);
-
-% 4th line: Ip, psiax1, psiax2, raxis1, raxis2
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',ip,faxs,0.,rmag,0.);
-
-% 5th line: zaxis1, zaxis2, psi_sep, R_xpoint, Z_xpoint
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',zmag+deltaz,0.0,fbnd,-1,-1);
-
-% 6th entry: F(psi) on nr equidistant psi mesh
-psieq=[0:1/(nrg-1):1];
-psi_efit_eff=faxs+psi_efit.*(fbnd-faxs);
-psieq_eff=faxs+psieq.*(fbnd-faxs);
-[G Gprime]=interpos(13,psi_efit_eff,-F,psieq_eff);
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',G(1:end-1));
-fprintf(fid,'%16.9e\n',G(end));
-
-% 7th entry: p(psi) on nr equidistant psi mesh
-[press pressprime]=interpos(13,psi_efit_eff,P,psieq_eff);
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',press(1:end-1));
-fprintf(fid,'%16.9e\n',press(end));
-
-% 8th entry: FF'(psi) on nr equidistant psi mesh
-y=G.*Gprime;
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',y(1:end-1));
-fprintf(fid,'%16.9e\n',y(end));
-
-% 9th entry: p'(psi) on nr equidistant psi mesh (in MKSA)
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',pressprime(1:end-1));
-fprintf(fid,'%16.9e\n',pressprime(end));
-
-% 10th entry: psi(i,j)
-psirz=faxs+psinrz.*(fbnd-faxs);
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',psirz);
-if mod(nrg*nzg,5)~=0
- fprintf(fid,'\n');
-end
-
-% 11th entry: q profile on nr equidistant psi mesh
-y=interpos(13,psi_efit,Q,psieq,1e-6);
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',y(1:end-1));
-fprintf(fid,'%16.9e\n',y(end));
-
-% 12th entry: (R,Z) plasma boundary and wall position
-npts=length(Rbnd);
-fprintf(fid,'%5d%5d\n',npts,5);
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',[Rbnd Zbnd+deltaz]');
-if mod(2*npts,5) ~= 0
- fprintf(fid,'\n');
-end
-rdel=+0.8*(R(2)-R(1));
-zdel=+0.8*(Z(2)-Z(1));
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',min(Rbnd)-rdel,min(Zbnd)+deltaz-zdel, ...
- max(Rbnd)+rdel,min(Zbnd)+deltaz-zdel,max(Rbnd)+rdel,max(Zbnd)+deltaz+zdel, ...
- min(Rbnd)-rdel,max(Zbnd)+deltaz+zdel,min(Rbnd)-rdel,min(Zbnd)+deltaz-zdel);
+ Rbnd=efitdata.Rbnd.data(:,index_efit_eff);
+ Zbnd=efitdata.Zbnd.data(:,index_efit_eff);
+ R0=efitdata.R0.data(index_efit_eff);
+ Z0=efitdata.Z0.data(index_efit_eff);
+ rmag=efitdata.rmag.data(index_efit_eff);
+ zmag=efitdata.zmag.data(index_efit_eff);
+ if deltaz==-99
+ deltaz=-zmag;
+ end
+ faxs=efitdata.faxs.data(index_efit_eff);
+ fbnd=efitdata.fbnd.data(index_efit_eff);
+ bvac=efitdata.bvac.data(index_efit_eff);
+ ip=efitdata.ip.data(index_efit_eff);
+ psi_efit=efitdata.psin;
+ F=efitdata.F.data(:,index_efit_eff);
+ P=efitdata.P.data(:,index_efit_eff);
+ Q=efitdata.Q.data(:,index_efit_eff);
+ kappa=efitdata.kappa.data(index_efit_eff);
+ q95=efitdata.q95.data(index_efit_eff);
+ if index_efit_eff<=length(efitdata.btpd.data); btpd=efitdata.btpd.data(index_efit_eff); else; btpd=0; end
+ if index_efit_eff<=length(efitdata.bttd.data); bttd=efitdata.bttd.data(index_efit_eff); else; bttd=0; end
+ if index_efit_eff<=length(efitdata.btnd.data); btnd=efitdata.btnd.data(index_efit_eff); else; btnd=0; end
+ if index_efit_eff<=length(efitdata.btpm.data); btpm=efitdata.btpm.data(index_efit_eff); else; btpm=0; end
+ if index_efit_eff<=length(efitdata.bttm.data); bttm=efitdata.bttm.data(index_efit_eff); else; bttm=0; end
+ if index_efit_eff<=length(efitdata.btnm.data); btnm=efitdata.btnm.data(index_efit_eff); else; btnm=0; end
+ xli=efitdata.xli.data(index_efit_eff);
+ [R,Z,psinrz]=psinrzjet(shot,time_efit,nrg,nzg,efitlab,uid{iexefit},seqd(iexefit),ncont,efitdata.sspr,efitdata.sspi,deltaz);
+ nzg=abs(nzg);
-% Some useful data to compare with recomputed equilibria
-fprintf(fid,'%18.8e psiedge-psiax\n',fbnd-faxs);
-fprintf(fid,'%18.8e r-magaxe\n',rmag);
-fprintf(fid,'%18.8e z-magaxe\n',zmag+deltaz);
-fprintf(fid,'%18.8e z0 (zaver)\n',0.5.*(min(Zbnd)+max(Zbnd))+deltaz);
-fprintf(fid,'%18.8e r-major\n',R0);
-fprintf(fid,'%18.8e B0\n',bvac);
-fprintf(fid,'%18.8e CURRT -> I-p [A]: %18.8e \n',ip/4e-7/pi,ip);
-fprintf(fid,'%18.8e kappa\n',kappa);
-fprintf(fid,'%18.8e q_0\n',Q(1));
-y=interpos(13,psi_efit,Q,[0:0.01:1],1e-6);
-fprintf(fid,'%18.8e q_edge, from int: %18.8e\n',Q(end),y(end));
-fprintf(fid,'%18.8e q_95, from int: %18.8e\n',q95,y(96));
-[qmin ind]=min(y);
-fprintf(fid,'%18.8e q_min, at psinorm= %18.8e\n',qmin,y(ind));
-fprintf(fid,'%18.8e beta_pol(wdia)\n',btpd);
-fprintf(fid,'%18.8e beta_pol(efit)\n',btpm);
-fprintf(fid,'%18.8e beta_tor(wdia)\n',bttd);
-fprintf(fid,'%18.8e beta_tor(efit)\n',bttm);
-fprintf(fid,'%18.8e beta_N(wdia)\n',btnd);
-fprintf(fid,'%18.8e beta_N(efit)\n',btnm);
-fprintf(fid,'%18.8e li\n',xli);
+ % define file name
+ s = sprintf('%.6f',time_efit);
+ fname=['EQDSK.' num2str(shot) 't' s];
+ fnamefull=fullfile(savedir,fname);
-fprintf(fid,'\n%18.8e time\n',time);
-fprintf(fid,' %d shot number\n',shot);
-fprintf(fid,' efit version : %s\n',efitlab);
-fprintf(fid,' uid : %s, %s, %s\n',uid{1},uid{2},uid{3});
-fprintf(fid,' seq : %d, %d, %d\n',seqd(1:3));
+ % fill in eqdsk structures
+ eqd{it}.fnamefull = fnamefull;
+ eqd{it}.fname = fname;
+ eqd{it}.cocos = 7;
+ % 1st eqdsk line: 48 characters for file description and then, 3, nr, nz
+ tdate=date;
+ ss=['JET #' num2str(shot) ' t= ' num2str(time_efit) ' ' efitlab '?uid=' uid{iexefit} '+seq=' num2str(seqd(iexefit)) ', ' tdate];
+ if length(ss)<48
+ ss(end:48)=' ';
+ else
+ ss=ss(1:48);
+ end
+ eqd{it}.stitle = ss;
+ eqd{it}.nr = nrg;
+ eqd{it}.nz = nzg;
+ rmin=min(R);
+ rmax=max(R);
+ zmin=min(Z);
+ zmax=max(Z);
+ eqd{it}.rboxlen = rmax-rmin;
+ eqd{it}.zboxlen = zmax-zmin;
+ eqd{it}.r0 = R0;
+ eqd{it}.rboxleft = rmin;
+ eqd{it}.zmid = 0.5*(zmin+zmax)+deltaz;
+ eqd{it}.rmesh=linspace(eqd{it}.rboxleft,eqd{it}.rboxleft+eqd{it}.rboxlen,eqd{it}.nr)';
+ eqd{it}.zmesh=linspace(eqd{it}.zmid-eqd{it}.zboxlen/2,eqd{it}.zmid+eqd{it}.zboxlen/2,eqd{it}.nz)';
+ eqd{it}.raxis = rmag;
+ eqd{it}.zaxis = zmag+deltaz;
+ eqd{it}.psiaxis = faxs;
+ eqd{it}.psiedge = fbnd;
+ eqd{it}.b0 = bvac;
+ eqd{it}.ip = ip;
+ rdel=+0.8*(R(2)-R(1));
+ zdel=+0.8*(Z(2)-Z(1));
+ eqd{it}.rlim=[eqd{it}.rboxleft-rdel ; eqd{it}.rboxleft+eqd{it}.rboxlen+rdel ; ...
+ eqd{it}.rboxleft+eqd{it}.rboxlen+rdel ; eqd{it}.rboxleft-rdel ; ...
+ eqd{it}.rboxleft-rdel];
+ eqd{it}.zlim=[eqd{it}.zmid-eqd{it}.zboxlen/2-zdel ; eqd{it}.zmid-eqd{it}.zboxlen/2-zdel ; ...
+ eqd{it}.zmid+eqd{it}.zboxlen/2+zdel ; eqd{it}.zmid+eqd{it}.zboxlen/2+zdel ; ...
+ eqd{it}.zmid-eqd{it}.zboxlen/2-zdel];
+ eqd{it}.nblim = 5;
+ psieq=[0:1/(nrg-1):1];
+ eqd{it}.psimesh = psieq;
+ eqd{it}.rhopsi = sqrt(eqd{it}.psimesh);
+ psi_efit_eff=faxs+psi_efit.*(fbnd-faxs);
+ psieq_eff=faxs+psieq.*(fbnd-faxs);
+ [G Gprime]=interpos(13,psi_efit_eff,F,psieq_eff);
+ eqd{it}.F = G;
+ eqd{it}.FFprime = G.*Gprime;
+ [press pressprime]=interpos(13,psi_efit_eff,P,psieq_eff);
+ eqd{it}.p = press;
+ eqd{it}.pprime = pressprime;
+ psirz=faxs+psinrz.*(fbnd-faxs);
+ eqd{it}.psi = psirz;
+ eqd{it}.psirz = reshape(eqd{it}.psi,prod(size(eqd{it}.psi)),1);
+ y=interpos(13,psi_efit,Q,psieq,-0.01);
+ eqd{it}.q = y;
+ npts=length(Rbnd);
+ eqd{it}.nbbound = npts;
+ eqd{it}.rplas = Rbnd;
+ eqd{it}.zplas = Zbnd+deltaz;
+ % Some useful data to compare with recomputed equilibria
+ eqd{it}.extralines = [];
+ eqd{it}.extralines{end+1} = [' ' num2str(fbnd-faxs) ' psiedge-psiax '];
+ eqd{it}.extralines{end+1} = [' ' num2str(rmag) ' r-magaxe '];
+ eqd{it}.extralines{end+1} = [' ' num2str(zmag+deltaz) ' z-magaxe '];
+ eqd{it}.extralines{end+1} = [' ' num2str(0.5.*(min(Zbnd)+max(Zbnd))+deltaz) ' z0 (zaver) '];
+ eqd{it}.extralines{end+1} = [' ' num2str(R0) ' r-major '];
+ eqd{it}.extralines{end+1} = [' ' num2str(bvac) ' B0 '];
+ eqd{it}.extralines{end+1} = [' ' num2str(ip/4e-7/pi) ' CURRT -> I-p [A]: ' num2str(ip/4e-7/pi,ip)];
+ eqd{it}.extralines{end+1} = [' ' num2str(kappa) ' kappa '];
+ eqd{it}.extralines{end+1} = [' ' num2str(Q(1)) ' q_0 '];
+ y=interpos(13,psi_efit,Q,[0:0.01:1],1e-6);
+ eqd{it}.extralines{end+1} = [' ' num2str(Q(end)) ' q_edge, from int: ' num2str(y(end))];
+ eqd{it}.extralines{end+1} = [' ' num2str(q95) ' q_95, from int: ' num2str(y(96))];
+ [qmin ind]=min(y);
+ eqd{it}.extralines{end+1} = [' ' num2str(qmin) ' q_min, at psinorm= ' num2str(y(ind))];
+ eqd{it}.extralines{end+1} = [' ' num2str(btpd) ' beta_pol(wdia) '];
+ eqd{it}.extralines{end+1} = [' ' num2str(btpm) ' beta_pol(efit) '];
+ eqd{it}.extralines{end+1} = [' ' num2str(bttd) ' beta_tor(wdia) '];
+ eqd{it}.extralines{end+1} = [' ' num2str(bttm) ' beta_tor(efit) '];
+ eqd{it}.extralines{end+1} = [' ' num2str(btnd) ' beta_N(wdia) '];
+ eqd{it}.extralines{end+1} = [' ' num2str(btnm) ' beta_N(efit) '];
+ eqd{it}.extralines{end+1} = [' ' num2str(xli) ' li '];
+
+ eqd{it}.extralines{end+1} = [' ' num2str(time_efit) 'time '];
+ eqd{it}.extralines{end+1} = [' ' num2str(shot) 'shot number '];
+ eqd{it}.extralines{end+1} = ['efit version : %s ' num2str(efitlab)];
+ eqd{it}.extralines{end+1} = ['uid : ' num2str(uid{iexefit})];
+ eqd{it}.extralines{end+1} = ['seq : ' num2str(seqd(iexefit))];
+ % write to file (it write the positive Ip, B0 as well automatically
+ write_eqdsk(eqd{it}.fnamefull,eqd{it}); % cocos contained in eqdsk
-fclose(fid);
-disp(['wrote ',fname]);
+% $$$ fid=fopen(fnamefull,'w');
+% $$$
+% $$$ % 1st eqdsk line: 48 characters for file description and then, 3, nr, nz
+% $$$ tdate=date;
+% $$$ ss=['JET #' num2str(shot) ' t= ' num2str(time_efit) ' ' efitlab '?uid=' uid{iexefit} '+seq=' num2str(seqd(iexefit)) ', ' tdate];
+% $$$ if length(ss)<48
+% $$$ ss(end:48)=' ';
+% $$$ else
+% $$$ ss=ss(1:48);
+% $$$ end
+% $$$
+% $$$ fprintf(fid,'%s%4d%4d%4d\n',ss,3,nrg,nzg);
+% $$$
+% $$$ % 2nd line: rboxlen, zboxlen, r0, rboxlft, zboxmid
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',rmax-rmin,zmax-zmin,R0,rmin,0.5*(zmin+zmax)+deltaz);
+% $$$
+% $$$ % 3rd line: rmag, zmag, psimag, psiedge, B0
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',rmag,zmag+deltaz,faxs,fbnd,bvac);
+% $$$
+% $$$ % 4th line: Ip, psiax1, psiax2, raxis1, raxis2
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',ip,faxs,0.,rmag,0.);
+% $$$
+% $$$ % 5th line: zaxis1, zaxis2, psi_sep, R_xpoint, Z_xpoint
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',zmag+deltaz,0.0,fbnd,-1,-1);
+% $$$
+% $$$ % 6th entry: F(psi) on nr equidistant psi mesh
+% $$$ psieq=[0:1/(nrg-1):1];
+% $$$ psi_efit_eff=faxs+psi_efit.*(fbnd-faxs);
+% $$$ psieq_eff=faxs+psieq.*(fbnd-faxs);
+% $$$ [G Gprime]=interpos(13,psi_efit_eff,F,psieq_eff);
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',G(1:end-1));
+% $$$ fprintf(fid,'%16.9e\n',G(end));
+% $$$
+% $$$ % 7th entry: p(psi) on nr equidistant psi mesh
+% $$$ [press pressprime]=interpos(13,psi_efit_eff,P,psieq_eff);
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',press(1:end-1));
+% $$$ fprintf(fid,'%16.9e\n',press(end));
+% $$$
+% $$$ % 8th entry: FF'(psi) on nr equidistant psi mesh
+% $$$ y=G.*Gprime;
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',y(1:end-1));
+% $$$ fprintf(fid,'%16.9e\n',y(end));
+% $$$
+% $$$ % 9th entry: p'(psi) on nr equidistant psi mesh (in MKSA)
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',pressprime(1:end-1));
+% $$$ fprintf(fid,'%16.9e\n',pressprime(end));
+% $$$
+% $$$ % 10th entry: psi(i,j)
+% $$$ psirz=faxs+psinrz.*(fbnd-faxs);
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',psirz);
+% $$$ if mod(nrg*nzg,5)~=0
+% $$$ fprintf(fid,'\n');
+% $$$ end
+% $$$
+% $$$ % 11th entry: q profile on nr equidistant psi mesh
+% $$$ y=interpos(13,psi_efit,Q,psieq,1e-6);
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',y(1:end-1));
+% $$$ fprintf(fid,'%16.9e\n',y(end));
+% $$$
+% $$$ % 12th entry: (R,Z) plasma boundary and wall position
+% $$$ npts=length(Rbnd);
+% $$$ fprintf(fid,'%5d%5d\n',npts,5);
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',[Rbnd Zbnd+deltaz]');
+% $$$ if mod(2*npts,5) ~= 0
+% $$$ fprintf(fid,'\n');
+% $$$ end
+% $$$ rdel=+0.8*(R(2)-R(1));
+% $$$ zdel=+0.8*(Z(2)-Z(1));
+% $$$ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',min(Rbnd)-rdel,min(Zbnd)+deltaz-zdel, ...
+% $$$ max(Rbnd)+rdel,min(Zbnd)+deltaz-zdel,max(Rbnd)+rdel,max(Zbnd)+deltaz+zdel, ...
+% $$$ min(Rbnd)-rdel,max(Zbnd)+deltaz+zdel,min(Rbnd)-rdel,min(Zbnd)+deltaz-zdel);
+% $$$
+% $$$ % Some useful data to compare with recomputed equilibria
+% $$$ fprintf(fid,'%18.8e psiedge-psiax\n',fbnd-faxs);
+% $$$ fprintf(fid,'%18.8e r-magaxe\n',rmag);
+% $$$ fprintf(fid,'%18.8e z-magaxe\n',zmag+deltaz);
+% $$$ fprintf(fid,'%18.8e z0 (zaver)\n',0.5.*(min(Zbnd)+max(Zbnd))+deltaz);
+% $$$ fprintf(fid,'%18.8e r-major\n',R0);
+% $$$ fprintf(fid,'%18.8e B0\n',bvac);
+% $$$ fprintf(fid,'%18.8e CURRT -> I-p [A]: %18.8e \n',ip/4e-7/pi,ip);
+% $$$ fprintf(fid,'%18.8e kappa\n',kappa);
+% $$$ fprintf(fid,'%18.8e q_0\n',Q(1));
+% $$$ y=interpos(13,psi_efit,Q,[0:0.01:1],1e-6);
+% $$$ fprintf(fid,'%18.8e q_edge, from int: %18.8e\n',Q(end),y(end));
+% $$$ fprintf(fid,'%18.8e q_95, from int: %18.8e\n',q95,y(96));
+% $$$ [qmin ind]=min(y);
+% $$$ fprintf(fid,'%18.8e q_min, at psinorm= %18.8e\n',qmin,y(ind));
+% $$$ fprintf(fid,'%18.8e beta_pol(wdia)\n',btpd);
+% $$$ fprintf(fid,'%18.8e beta_pol(efit)\n',btpm);
+% $$$ fprintf(fid,'%18.8e beta_tor(wdia)\n',bttd);
+% $$$ fprintf(fid,'%18.8e beta_tor(efit)\n',bttm);
+% $$$ fprintf(fid,'%18.8e beta_N(wdia)\n',btnd);
+% $$$ fprintf(fid,'%18.8e beta_N(efit)\n',btnm);
+% $$$ fprintf(fid,'%18.8e li\n',xli);
+% $$$
+% $$$ fprintf(fid,'\n%18.8e time\n',time_efit);
+% $$$ fprintf(fid,' %d shot number\n',shot);
+% $$$ fprintf(fid,' efit version : %s\n',efitlab);
+% $$$ fprintf(fid,' uid : %s, %s, %s\n',uid{1},uid{2},uid{3});
+% $$$ fprintf(fid,' seq : %d, %d, %d\n',seqd(1:3));
+% $$$
+% $$$ fclose(fid);
+% $$$ disp(['wrote ',fnamefull]);
-if ncont>0
- figure;
- pos=get(gcf,'position');
- set(gcf,'position',[pos(1)+0.5*pos(3) 0.8*abs(pos(2)-pos(4)) pos(3) 2*pos(4)])
- subplot(3,1,1)
- plot(psi_efit,P/P(1),'-')
- hold on
- plot(psi_efit,F/F(end),'r-')
- ss=sprintf('%.4f',time);
- title(['JET #' num2str(shot) ' t= ' ss])
- legend('P/P(0)','F/F(edge)',3)
- subplot(3,1,2)
- plot(psieq,pressprime/abs(pressprime(1)),'-')
- hold on
- plot(psieq,G.*Gprime/abs(G(1).*Gprime(1)),'r-')
- legend('Pprime/|Pprime(0)|','F*Fprime/|F*Fprime(0)|',2)
- subplot(3,1,3)
- plot(psi_efit,Q,'-')
- hold on
- aa=axis;
- plot([0.95 0.95],[aa(3) aa(4)],'k--')
- grid on
- xlabel('\psi/\psi_{edge}')
+ if ncont>0
+ figure;
+ pos=get(gcf,'position');
+ set(gcf,'position',[pos(1)+0.5*pos(3) 0.8*abs(pos(2)-pos(4)) pos(3) 2*pos(4)])
+ subplot(3,1,1)
+ plot(psi_efit,P/P(1),'-')
+ hold on
+ plot(psi_efit,F/F(end),'r-')
+ ss=sprintf('%.4f',time_efit);
+ title(['JET #' num2str(shot) ' t= ' ss])
+ legend('P/P(0)','F/F(edge)',3)
+ subplot(3,1,2)
+ plot(psieq,pressprime/abs(pressprime(1)),'-')
+ hold on
+ plot(psieq,G.*Gprime/abs(G(1).*Gprime(1)),'r-')
+ legend('Pprime/|Pprime(0)|','F*Fprime/|F*Fprime(0)|',2)
+ subplot(3,1,3)
+ plot(psi_efit,Q,'-')
+ hold on
+ aa=axis;
+ plot([0.95 0.95],[aa(3) aa(4)],'k--')
+ grid on
+ xlabel('\psi/\psi_{edge}')
+ end
end
diff --git a/crpptbx/JET/jet_requests_mapping.m b/crpptbx/JET/jet_requests_mapping.m
index 9a4f0a62..f843d20f 100644
--- a/crpptbx/JET/jet_requests_mapping.m
+++ b/crpptbx/JET/jet_requests_mapping.m
@@ -74,7 +74,8 @@ switch lower(data_request)
mapping.method = 'expression';
mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''delta_bottom''; ' ...
'gdat_tmp=gdat_jet(shot,params_eff);params_eff.data_request=''delta_top'';' ...
- 'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.data = 0.5.*(gdat_tmp.data+gdat_tmp2.data);'];
+ 'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.data = 0.5.*(gdat_tmp.data+gdat_tmp2.data);' ...
+ 'gdat_tmp.label=gdat_data.label; gdat_tmp.gdat_request=gdat_data.gdat_request;'];
case 'delta_top'
mapping.label = 'delta\_top';
mapping.timedim = 1;
@@ -120,6 +121,11 @@ switch lower(data_request)
mapping.method = 'signal';
mapping.expression = [{'ppf'},{'EHTR'},{'ELON'}];
mapping.expression = [{'ppf'},{'EFIT'},{'ELON'}];
+ case {'li', 'l_i'}
+ mapping.timedim = 1;
+ mapping.label = 'l_i_3';
+ mapping.method = 'signal';
+ mapping.expression = [{'ppf'},{'EFIT'},{'LI3M'}];
case 'mhd'
mapping.timedim = 1;
mapping.label = 'n=1 and n=2 amplitude';
@@ -192,6 +198,21 @@ switch lower(data_request)
'gdat_tmp.data = gdat_tmp.data./(tmp_data2+1e-5);'];
case 'ni'
mapping.method = 'switchcase'; % especially since might have option fit, etc
+ case {'phi_tor', 'phitor', 'toroidal_flux'}
+ mapping.label = 'toroidal_flux';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'PPF'},{'EQUI'},{'PHIT'}]; % note this is only chain2, so should check with efit...
+ case 'p_ohmic'
+ mapping.label = 'p\_ohmic';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'PPF'},{'EFIT'},{'POHM'}];
+ case {'p_rad', 'prad'}
+ mapping.label = 'p\_rad';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'PPF'},{'PRAD'},{'PRAD'}];
case 'powers'
mapping.timedim = 1;
mapping.label = 'various powers';
@@ -292,6 +313,11 @@ switch lower(data_request)
mapping.timedim = 2;
mapping.label = 'volume\_norm';
mapping.method = 'switchcase';
+ case {'wmhd', 'w_mhd'}
+ mapping.label = 'W_mhd';
+ mapping.timedim = 1;
+ mapping.method = 'signal';
+ mapping.expression = [{'PPF'},{'EFIT'},{'WP'}];
case 'zeff'
mapping.label = 'zeff from KS3';
mapping.timedim = 1;
--
GitLab