From 702dae64594b905799cc318131efab027ce8102e Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Thu, 11 Oct 2001 09:10:24 +0000
Subject: [PATCH] fix minor bug with geteqdsk
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1847 d63d8f72-b253-0410-a779-e742ad2e26cf
---
JET/geteqdskJET.m | 167 +++++++++++++++++++++++++---------------------
1 file changed, 90 insertions(+), 77 deletions(-)
diff --git a/JET/geteqdskJET.m b/JET/geteqdskJET.m
index 74902f6a..5585659f 100644
--- a/JET/geteqdskJET.m
+++ b/JET/geteqdskJET.m
@@ -1,9 +1,11 @@
function [efitdata]=geteqdskJET(shot,time,nrg,nzg,savedir,deltaz,efitlab,uid,seqd,varargin);
%
% save eqdsk file to savedir/JET_EQDSK_shot_ttime
+% returns data used in structure efitdata (efitdata.shot, efitdata.Rbnd, efitdataZbnd....)
%
% examples:
-% geteqdsk(shot,time,33,65,'/tmp')
+% efitdata=geteqdskJET(shot,time,33,65,'/tmp')
+% efitdata=geteqdskJET(shot,time,33,65,'/tmp',[],[],[],[],1,efitdata); % gives data so no need to load again
%
% INPUTS :
% shot : shot number
@@ -56,17 +58,19 @@ end
iread=1;
if nargin>=11 & ~isempty(varargin{2})
efitdata=varargin{2};
- if shot~=efitdata.shot
+ if shot==efitdata.shot
+ iread=0;
+ else
disp(['shot=' num2str(shot) ' different from efitdata.shot=' num2str(efitdata.shot)])
- warning('in geteqdskJET')
- return
+ warning('reload data with new shot')
end
- iread=0;
end
if iread==1
+ efitdata.shot=shot;
% get data needed
efitdata.Rbnd=gdat(shot,['ppf/' efitlab '/RBND' s_extra],0,'JET');
+ efitdata.tefit=efitdata.Rbnd.t;
efitdata.Zbnd=gdat(shot,['ppf/' efitlab '/ZBND' s_extra],0,'JET');
efitdata.R0=gdat(shot,['ppf/' efitlab '/RGEO' s_extra],0,'JET');
% $$$ efitdata.Z0=gdat(shot,['ppf/' efitlab '/ZO' s_extra],0,'JET');
@@ -77,57 +81,57 @@ if iread==1
efitdata.bvac=gdat(shot,['ppf/' efitlab '/bvac' s_extra],0,'JET');
efitdata.ip=gdat(shot,['ppf/' efitlab '/xip' s_extra],0,'JET');
efitdata.F=gdat(shot,['ppf/' efitlab '/F' s_extra],0,'JET');
+ efitdata.psin=efitdata.F.x;
efitdata.P=gdat(shot,['ppf/' efitlab '/P' s_extra],0,'JET');
efitdata.Q=gdat(shot,['ppf/' efitlab '/Q' s_extra],0,'JET');
efitdata.kappa=gdat(shot,['ppf/' efitlab '/ELON' s_extra],0,'JET');
efitdata.q95=gdat(shot,['ppf/' efitlab '/q95' s_extra],0,'JET');
efitdata.btpd=gdat(shot,['ppf/' efitlab '/btpd' s_extra],0,'JET');
- efitdata.bttd=gdat(shot,['ppf/' efitlab '/efitdata.bttd' s_extra],0,'JET');
+ efitdata.bttd=gdat(shot,['ppf/' efitlab '/bttd' s_extra],0,'JET');
efitdata.btnd=gdat(shot,['ppf/' efitlab '/btnd' s_extra],0,'JET');
efitdata.btpm=gdat(shot,['ppf/' efitlab '/btpm' s_extra],0,'JET');
- efitdat.bttm=gdat(shot,['ppf/' efitlab '/bttm' s_extra],0,'JET');
+ efitdata.bttm=gdat(shot,['ppf/' efitlab '/bttm' s_extra],0,'JET');
efitdata.btnm=gdat(shot,['ppf/' efitlab '/btnm' s_extra],0,'JET');
efitdata.xli=gdat(shot,['ppf/' efitlab '/xli' s_extra],0,'JET');
efitdata.sspr=gdat(shot,['ppf/' efitlab '/sspr' s_extra],0,'JET');
efitdata.sspi=gdat(shot,['ppf/' efitlab '/sspi' s_extra],0,'JET');
end
-tefit=efitdata.Rbnd.t;
+tefit=efitdata.tefit;
[zz index_efit]=min(abs(tefit-time));
time_efit=tefit(index_efit);
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
disp(['efit at t=' num2str(time_efit)])
disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
-efitdata.Rbnd=efitdata.Rbnd.data(:,index_efit);
-efitdata.Zbnd=efitdata.Zbnd.data(:,index_efit);
-efitdata.R0=efitdata.R0.data(index_efit);
-% $$$ efitdata.Z0=efitdata.Z0.data(index_efit);
-efitdata.rmag=efitdata.rmag.data(index_efit);
-efitdata.zmag=efitdata.zmag.data(index_efit);
+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=-efitdata.zmag;
end
-efitdata.faxs=efitdata.faxs.data(index_efit);
-efitdata.fbnd=efitdata.fbnd.data(index_efit);
-efitdata.bvac=-efitdata.bvac.data(index_efit);
-efitdata.ip=-efitdata.ip.data(index_efit);
-psi_efit=efitdata.F.x;
-efitdata.F=efitdata.F.data(:,index_efit);
-efitdata.P=efitdata.P.data(:,index_efit);
-efitdata.Q=efitdata.Q.data(:,index_efit);
-efitdata.kappa=efitdata.kappa.data(index_efit);
-efitdata.q95=efitdata.q95.data(index_efit);
-efitdata.btpd=efitdata.btpd.data(index_efit);
-efitdata.bttd=efitdata.bttd.data(index_efit);
-efitdata.btnd=efitdata.btnd.data(index_efit);
-efitdata.btpm=efitdata.btpm.data(index_efit);
-efitdat.bttm=efitdat.bttm.data(index_efit);
-efitdata.btnm=efitdata.btnm.data(index_efit);
-efitdata.xli=efitdata.xli.data(index_efit);
+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);
+btpd=efitdata.btpd.data(index_efit);
+bttd=efitdata.bttd.data(index_efit);
+btnd=efitdata.btnd.data(index_efit);
+btpm=efitdata.btpm.data(index_efit);
+bttm=efitdata.bttm.data(index_efit);
+btnm=efitdata.btnm.data(index_efit);
+xli=efitdata.xli.data(index_efit);
-
-[R,Z,psinrz]=psirz(shot,time,nrg,nzg,efitlab,uid,seqd,0,efitdata.sspr,efitdata.sspi);
+[R,Z,psinrz]=psirz(shot,time,nrg,nzg,efitlab,uid,seqd,ncont,efitdata.sspr,efitdata.sspi);
% define file name
s = sprintf('%.6f',time);
@@ -150,27 +154,27 @@ 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,efitdata.R0,rmin,0.5*(zmin+zmax)+deltaz);
+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',efitdata.rmag,efitdata.zmag+deltaz,efitdata.faxs,efitdata.fbnd,efitdata.bvac);
+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',efitdata.ip,efitdata.faxs,0.,efitdata.rmag,0.);
+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',efitdata.zmag+deltaz,0.0,efitdata.fbnd,-1,-1);
+fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',zmag+deltaz,0.0,fbnd,-1,-1);
-% 6th entry: efitdata.F(psi) on nr equidistant psi mesh
+% 6th entry: F(psi) on nr equidistant psi mesh
psieq=[0:1/(nrg-1):1];
-psi_efit_eff=efitdata.faxs+psi_efit.*(efitdata.fbnd-efitdata.faxs);
-psieq_eff=efitdata.faxs+psieq.*(efitdata.fbnd-efitdata.faxs);
-[G Gprime]=interpos(13,psi_efit_eff,-efitdata.F,psieq_eff);
+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,efitdata.P,psieq_eff);
+[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));
@@ -184,52 +188,52 @@ 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=efitdata.faxs+psinrz.*(efitdata.fbnd-efitdata.faxs);
+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,efitdata.Q,psieq,1e-6);
+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(efitdata.Rbnd);
+npts=length(Rbnd);
fprintf(fid,'%5d%5d\n',npts,5);
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',[efitdata.Rbnd efitdata.Zbnd+deltaz]');
+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=2*(R(2)-R(1));
zdel=2*(Z(2)-Z(1));
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',min(efitdata.Rbnd)-rdel,min(efitdata.Zbnd)+deltaz-zdel, ...
- max(efitdata.Rbnd)+rdel,min(efitdata.Zbnd)+deltaz-zdel,max(efitdata.Rbnd)+rdel,max(efitdata.Zbnd)+deltaz+zdel, ...
- min(efitdata.Rbnd)-rdel,max(efitdata.Zbnd)+deltaz+zdel,min(efitdata.Rbnd)-rdel,min(efitdata.Zbnd)+deltaz-zdel);
+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',efitdata.fbnd-efitdata.faxs);
-fprintf(fid,'%18.8e r-magaxe\n',efitdata.rmag);
-fprintf(fid,'%18.8e z-magaxe\n',efitdata.zmag+deltaz);
-fprintf(fid,'%18.8e z0 (zaver)\n',0.5.*(min(efitdata.Zbnd)+max(efitdata.Zbnd))+deltaz);
-fprintf(fid,'%18.8e r-major\n',efitdata.R0);
-fprintf(fid,'%18.8e B0\n',efitdata.bvac);
-fprintf(fid,'%18.8e CURRT -> I-p [A]: %18.8e \n',efitdata.ip/4e-7/pi,efitdata.ip);
-fprintf(fid,'%18.8e kappa\n',efitdata.kappa);
-fprintf(fid,'%18.8e q_0\n',efitdata.Q(1));
-y=interpos(13,psi_efit,efitdata.Q,[0:0.01:1],1e-6);
-fprintf(fid,'%18.8e q_edge, from int: %18.8e\n',efitdata.Q(end),y(end));
-fprintf(fid,'%18.8e q_95, from int: %18.8e\n',efitdata.q95,y(96));
+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',efitdata.btpd);
-fprintf(fid,'%18.8e beta_pol(efit)\n',efitdata.btpm);
+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',efitdat.bttm);
-fprintf(fid,'%18.8e beta_N(wdia)\n',efitdata.btnd);
-fprintf(fid,'%18.8e beta_N(efit)\n',efitdata.btnm);
-fprintf(fid,'%18.8e li\n',efitdata.xli);
+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);
fprintf(fid,' %d shot number\n',shot);
@@ -241,17 +245,26 @@ fclose(fid);
disp(['wrote ',fname]);
if ncont>0
- figure
- contour(r,z,psinrz',ncont);
+ 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
- [h1, h2]=contour(r,z,psinrz',[1 1]);
- for i=1:length(h2)
- set(h2(i),'LineWidth',2);
- end
+ plot(psi_efit,F/F(end),'r-')
ss=sprintf('%.4f',time);
title(['JET #' num2str(shot) ' t= ' ss])
- xlabel('R [m]')
- ylabel('Z [m]')
- axis equal
- axis([rmin rmax zmin zmax])
+ 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
--
GitLab