From dbee37dbfffe70ab7826a66134158dd171e578d2 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 10 Oct 2001 15:11:26 +0000
Subject: [PATCH] geteqdsk works and minor bug in loadJETdata fixed
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1844 d63d8f72-b253-0410-a779-e742ad2e26cf
---
JET/geteqdskJET.m | 44 +++++++++++++++++++++++++++++---------------
JET/loadJETdata.m | 9 +++++++--
gdat.m | 2 +-
3 files changed, 37 insertions(+), 18 deletions(-)
diff --git a/JET/geteqdskJET.m b/JET/geteqdskJET.m
index 567f0c26..6c5d2692 100644
--- a/JET/geteqdskJET.m
+++ b/JET/geteqdskJET.m
@@ -28,7 +28,7 @@ function geteqdskJET(shot,time,nrg,nzg,savedir,deltaz,efitlab,uid,seqd,varargin)
if nargin<=4 | isempty(savedir)
savedir='./';
end
-if ~unix(['test -d ' savedir])
+if unix(['test -d ' savedir])
disp(['Problems in geteqdskJET, savedir=' savedir ' is not a directory'])
return
end
@@ -51,6 +51,8 @@ if nargin>=10 & ~isempty(varargin{1})
ncont=varargin{1};
end
+iread=1;
+if iread==1
% get data needed
Rbnd=gdat(shot,['ppf/' efitlab '/RBND' s_extra],0,'JET');
@@ -60,12 +62,14 @@ time_efit=tefit(index_efit);
disp(['efit at t=' num2str(time_efit)])
Rbnd=Rbnd.data(:,index_efit);
-Zbnd=gdat(shot,['ppf/' efitlab '/RBND' s_extra],0,'JET');
+Zbnd=gdat(shot,['ppf/' efitlab '/ZBND' s_extra],0,'JET');
Zbnd=Zbnd.data(:,index_efit);
R0=gdat(shot,['ppf/' efitlab '/RGEO' s_extra],0,'JET');
R0=R0.data(index_efit);
-Z0=gdat(shot,['ppf/' efitlab '/ZO' s_extra],0,'JET');
-Z0=Z0.data(index_efit);
+% $$$ Z0=gdat(shot,['ppf/' efitlab '/ZO' s_extra],0,'JET');
+% $$$ Z0=Z0.data(index_efit);
+rmag=gdat(shot,['ppf/' efitlab '/rmag' s_extra],0,'JET');
+rmag=rmag.data(index_efit);
zmag=gdat(shot,['ppf/' efitlab '/zmag' s_extra],0,'JET');
zmag=zmag.data(index_efit);
if deltaz==-99
@@ -76,9 +80,9 @@ faxs=faxs.data(index_efit);
fbnd=gdat(shot,['ppf/' efitlab '/fbnd' s_extra],0,'JET');
fbnd=fbnd.data(index_efit);
bvac=gdat(shot,['ppf/' efitlab '/bvac' s_extra],0,'JET');
-bvac=bvac.data(index_efit);
+bvac=-bvac.data(index_efit);
xip=gdat(shot,['ppf/' efitlab '/xip' s_extra],0,'JET');
-xip=xip.data(index_efit);
+xip=-xip.data(index_efit);
F=gdat(shot,['ppf/' efitlab '/F' s_extra],0,'JET');
psi_efit=F.x;
F=F.data(:,index_efit);
@@ -86,7 +90,7 @@ P=gdat(shot,['ppf/' efitlab '/P' s_extra],0,'JET');
P=P.data(:,index_efit);
Q=gdat(shot,['ppf/' efitlab '/Q' s_extra],0,'JET');
Q=Q.data(:,index_efit);
-kappa=gdat(shot,['ppf/' efitlab '/kappa' s_extra],0,'JET');
+kappa=gdat(shot,['ppf/' efitlab '/ELON' s_extra],0,'JET');
kappa=kappa.data(index_efit);
q95=gdat(shot,['ppf/' efitlab '/q95' s_extra],0,'JET');
q95=q95.data(index_efit);
@@ -104,6 +108,11 @@ sspi=gdat(shot,['ppf/' efitlab '/sspi' s_extra],0,'JET');
[R,Z,psinrz]=psirz(shot,time,nrg,nzg,efitlab,uid,seqd,0,sspr,sspi);
+% $$$ save geteqdskmat
+else
+% $$$ load geteqdskmat
+end
+
% define file name
s = sprintf('%.6f',time);
fname=sprintf('%s/JET_EQDSK_%st%s',savedir,num2str(shot),s);
@@ -111,8 +120,13 @@ 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 '+seq=' seqd ', ' tdate];
-ss(end:48)=' ';
+ss=['JET #' num2str(shot) ' t= ' num2str(time_efit) ' ' efitlab '?uid=' uid '+seq=' num2str(seqd) ', ' 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
@@ -123,17 +137,17 @@ 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);
% 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);
+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',-xip,faxs,0.,rmag,0.);
+fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',xip,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];
-[G Gprime]=interpos(13,psi_efit,F,psieq);
+[G Gprime]=interpos(13,psi_efit,-F,psieq);
fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',G(1:end-1));
fprintf(fid,'%16.9e\n',G(end));
@@ -152,8 +166,8 @@ 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;
-fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',psirz(:,:));
+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
@@ -177,7 +191,7 @@ fprintf(fid,'%16.9e%16.9e%16.9e%16.9e%16.9e\n',min(Rbnd)-rdel,min(Zbnd)+deltaz-z
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 psimax\n',fbnd);
+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);
diff --git a/JET/loadJETdata.m b/JET/loadJETdata.m
index caed7615..a1e31541 100644
--- a/JET/loadJETdata.m
+++ b/JET/loadJETdata.m
@@ -231,7 +231,12 @@ trace.name=[];
if size(data_type_eff,1)==2
% in case node name was given in 2 parts directly (as old way)
ii1=strmatch(data_type_eff(1),JETsiglocation(1,:),'exact');
- index=find(strmatch(data_type_eff_noext,JETsiglocation(2,ii1),'exact'));
+ iiindex=strmatch(data_type_eff_noext,JETsiglocation(2,ii1),'exact');
+ if ~isempty(iiindex)
+ index=ii1(iiindex);
+ else
+ index=[];
+ end
if isempty(index)
disp('********************')
disp('trace not yet registered.')
@@ -475,7 +480,7 @@ switch JETkeywrdcase{index}
for i=starti:endi
trace.R(i,:)=interp1(radius.t,radius.data(i,:),trace.t);
end
-keyboard
+
otherwise
disp('case not yet defined')
diff --git a/gdat.m b/gdat.m
index 230e0084..f59ada50 100644
--- a/gdat.m
+++ b/gdat.m
@@ -117,7 +117,7 @@ else
end
% PLOT DATA (if required)
-if doplot==1
+if doplot==1 & length(trace.data)>1
figure;zoom on
if length(size(trace.data))<=2
plot(trace.t,trace.data);
--
GitLab