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

rm geteqdskJET.m and add to CHEASEgui subfolder for consistency 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7752 d63d8f72-b253-0410-a779-e742ad2e26cf
parent bc16fc47
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
crpptbx/JET/geteqdskJET.m deleted 100644 → 0
+ 0
− 412
View file @ bc16fc47
function [efitdata,eqd]=geteqdskJET(shot,times_in,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_t(times_in(:))
% returns data used in structure efitdata (efitdata.shot, efitdata.Rbnd, efitdataZbnd....)
%
% examples:
% 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
% times_in : single time slice of analysis
%
% OPTIONAL
% 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)
% otherwise: shift Z position by this amount
% efitlab : 'efit' or 'eftm' (default: 'efit')
% uid{i}: user id if different from main database (default: 'jetppf')
% seq(i): sequence number if not last one required (default: 0)
% i: 1 (efit), 2 (chain2: lid2), 3 (equi)
%
% varargin{1}: plot option: 0: do not plot contours, >0 plot contour with varargin{1} nb of contours (60 is good)
% varargin{2}: efitdata in input (no need to load again if shot=efitdata.shot)
%
% OUTPUTS
% efitdata: structure containing all data vs time so that one can call geteqdsk for another time quicker for same shot
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if ~exist('nrg') || isempty(nrg)
nrg = 129;
end
if ~exist('nzg') || isempty(nzg)
nzg = -129;
end
if nargin<=4 | isempty(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'])
return
end
if nargin<=5 | isempty(deltaz)
deltaz=0;
end
if nargin<=6 | isempty(efitlab)
efitlab='efit';
end
iexefit=1;
iexchain2=2;
iexequi=3;
uid0{iexefit}='jetppf';
uid0{iexchain2}='jetppf';
uid0{iexequi}='jetppf';
seqd0(1:iexequi)=0;
if nargin>=8 & ~isempty(uid)
for i=length(uid)+1:length(uid0)
uid{i}=uid0{i}
end
else
uid=uid0;
end
if nargin>=9 & ~isempty(seqd)
for i=length(seqd)+1:length(seqd0)
seqd(i)=seqd0(i);
end
else
seqd=seqd0;
end
for i=1:length(uid)
s_extra{i}=['?uid=' uid{i} '+seq=' num2str(seqd(i))];
end
ncont=0;
if nargin>=10 & ~isempty(varargin{1})
ncont=varargin{1};
end
iread=1;
if nargin>=11 & ~isempty(varargin{2})
efitdata=varargin{2};
if shot==efitdata.shot
iread=0;
else
disp(['shot=' num2str(shot) ' different from efitdata.shot=' num2str(efitdata.shot)])
warning('reload data with new shot')
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}]}],'machine',exp_machine,'doplot',do_plot);
efitdata.tefit=efitdata.Rbnd.t;
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}]}],'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;
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);
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);
% define file name
s = sprintf('%.6f',time_efit);
fname=['EQDSK.' num2str(shot) 't' s];
fnamefull=fullfile(savedir,fname);
% 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
% $$$ 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_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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment