to create JET eqdsk

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1843 d63d8f72-b253-0410-a779-e742ad2e26cf

JET/geteqdskJET.m

+function geteqdskJET(shot,time,nrg,nzg,savedir,deltaz,efitlab,uid,seqd,varargin);
+%
+% save eqdsk file to savedir/JET_EQDSK_shot_ttime
+% 
+% examples:
+%          geteqdsk(shot,time,33,65,'/tmp')
+%
+% INPUTS :
+% shot : shot number
+% time   : single time slice of analysis
+% nrg, nzg: nb of R and Z points for R,Z grid
+% OPTIONAL
+% savedir: directory to save eqdsk file (default: './')
+% 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: user id if different from main database (default: 'jetppf')
+% seq: sequence number if not last one required (default: 0)
+%
+% varargin{1}: plot option: 0: do not plot contours, >0 plot contour with varargin{1} nb of contours (60 is good)
+%
+% OUTPUTS
+%
+%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+if nargin<=4 | isempty(savedir)
+  savedir='./';
+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
+if nargin<=7 | isempty(uid)
+  uid='jetppf';
+end
+if nargin<=8 | isempty(seqd)
+  seqd=0;
+end
+s_extra=['?uid=' uid '+seq=' num2str(seqd)];
+
+ncont=0;
+if nargin>=10 & ~isempty(varargin{1})
+  ncont=varargin{1};
+end
+
+% get data needed
+Rbnd=gdat(shot,['ppf/' efitlab '/RBND' s_extra],0,'JET');
+
+tefit=Rbnd.t;
+[zz index_efit]=min(abs(tefit-time));
+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=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);
+zmag=gdat(shot,['ppf/' efitlab '/zmag' s_extra],0,'JET');
+zmag=zmag.data(index_efit);
+if deltaz==-99
+  deltaz=-zmag;
+end
+faxs=gdat(shot,['ppf/' efitlab '/faxs' s_extra],0,'JET');
+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);
+xip=gdat(shot,['ppf/' efitlab '/xip' s_extra],0,'JET');
+xip=xip.data(index_efit);
+F=gdat(shot,['ppf/' efitlab '/F' s_extra],0,'JET');
+psi_efit=F.x;
+F=F.data(:,index_efit);
+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=kappa.data(index_efit);
+q95=gdat(shot,['ppf/' efitlab '/q95' s_extra],0,'JET');
+q95=q95.data(index_efit);
+btpd=gdat(shot,['ppf/' efitlab '/btpd' s_extra],0,'JET');
+btpd=btpd.data(index_efit);
+bttd=gdat(shot,['ppf/' efitlab '/bttd' s_extra],0,'JET');
+bttd=bttd.data(index_efit);
+btnd=gdat(shot,['ppf/' efitlab '/btnd' s_extra],0,'JET');
+btnd=btnd.data(index_efit);
+xli=gdat(shot,['ppf/' efitlab '/xli' s_extra],0,'JET');
+xli=xli.data(index_efit);
+
+sspr=gdat(shot,['ppf/' efitlab '/sspr' s_extra],0,'JET');
+sspi=gdat(shot,['ppf/' efitlab '/sspi' s_extra],0,'JET');
+
+[R,Z,psinrz]=psirz(shot,time,nrg,nzg,efitlab,uid,seqd,0,sspr,sspi);
+
+% define file name
+s = sprintf('%.6f',time);
+fname=sprintf('%s/JET_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 '+seq=' seqd ', ' tdate];
+ss(end:48)=' ';
+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);
+
+% 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',-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);
+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,P,psieq);
+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;
+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=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(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   psimax\n',fbnd);
+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',xip/4e-7/pi,xip);
+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\n',btpd);
+fprintf(fid,'%18.8e   beta_tor\n',bttd);
+fprintf(fid,'%18.8e   beta_N\n',btnd);
+fprintf(fid,'%18.8e   li\n',xli);
+
+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\n',uid);
+fprintf(fid,'  seq : %d\n',seqd);
+
+fclose(fid);
+disp(['wrote ',fname]);
+
+if ncont>0
+  figure
+  contour(r,z,psinrz',ncont);
+  hold on
+  [h1, h2]=contour(r,z,psinrz',[1 1]);
+  for i=1:length(h2)
+    set(h2(i),'LineWidth',2);
+  end
+  ss=sprintf('%.4f',time);
+  title(['JET #' num2str(shot) ' t= ' ss])
+  xlabel('R [m]')
+  ylabel('Z [m]')
+  axis equal
+  axis([rmin rmax zmin zmax])
+end