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

faster mapflux and add diag option is psirz 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1855 d63d8f72-b253-0410-a779-e742ad2e26cf
parent c9b1bad8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
JET/mapflux_fast.m 0 → 100644
+ 98
− 0
View file @ 049fdbf0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function flux=mapflux_fast(cij,ri,zi,r,z,varargin)
% mapflux : carte de flux normalise (r,z) sur jet
%
% syntaxe : flux=mapflux(cij,ri,zi,r,z)
%
% entrees :
% cij : coefficients bspline (efit/sspr(nr+nz+1:nr+nz+nc))
% ri : points de controle pour r (efit/sspr(1:nr))
% zi : points de controle pour z (efit/sspr(nr+1:nr+nz))
% r : grille en r pour la carte de flux (vecteur)
% z : grille en z pour la carte de flux (vecteur)
%
% varargin{1}=1: output as diagonal of flux to get psin on (R,Z) set of points
% 0: output full matrix (default)
% varargin{1}=1: limits spline calculation to relevant ri,zi points
% (useful when computing on a few (r,z) points)
% (default)
% =0: compute all splines (slightly faster when whole grid needed)
%
% sortie :
% flux : matrice donnant le flux normalise sur la grille (r,z)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% auteur : patrick maget
% date : 13/10/2000
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
nr=length(ri);
nz=length(zi);
nrout=length(r);
nzout=length(z);
absplr=zeros(nrout,nr-4);
if nargin>6 & varargin{2}==0
i1=1;
i2=nr-4;
else
ii=find(ri>=min(r));i1=ii(1);
ii=find(ri>max(r));i2=ii(1);
end
for i=max(1,i1-4):min(nr-4,i2-1)
absplr(:,i)=bsplinev(ri,r,3,i)';
end
absplz=zeros(nz-4,nzout);
if nargin>6 & varargin{2}==0
i1=1;
i2=nz-4;
else
ii=find(zi>=min(z));i1=ii(1);
ii=find(zi>max(z));i2=ii(1);
end
for j=max(1,i1-4):min(nz-4,i2-1)
absplz(j,:)=bsplinev(zi,z,3,j);
end
cij2D=reshape(cij,nz-4,nr-4)';
flux=absplr*(cij2D*absplz);
if ~isempty(varargin) & varargin{1}
flux=diag(flux);
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% fonction b-spline
% bki=bspline(ti,t,k,i)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% auteur : patrick maget
% date : 14/12/99
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% retourne pour la valeur de i,
% le vecteur de poids correspondant au vecteur de demande t
% seuls les elements de t situes entre ti(i) et ti(i+k+1) ont un poids non nul
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bki=bsplinev(ti,t,k,i)
if (i>=length(ti))
bki=0.0;
elseif (k==0)
bki= ((t>=ti(i))&(t<ti(i+1)));
% if ((t<ti(i)) | (t>=ti(i+1)))
% bki=0.0;
% else
% bki=1.0;
% end
else
if (ti(i+k)==ti(i))
fac1=0.0;
else
fac1=(t-ti(i))/(ti(i+k)-ti(i));
end
if (ti(i+k+1)==ti(i+1))
fac2=0.0;
else
fac2=(ti(i+k+1)-t)/(ti(i+k+1)-ti(i+1));
end
bki=fac1.*bsplinev(ti,t,k-1,i)+fac2.*bsplinev(ti,t,k-1,i+1);
end
JET/psinrzjet.m 0 → 100644
+ 128
− 0
View file @ 049fdbf0
function [r,z,psinrz,sspr,sspi]=psinrzjet(shot,time,nrg_rg,nzg_zg,efitlab,uid,seqd,varargin);
%
% psirz : reconstruction des surfaces de flux
%
% ce programme utilise les donnees de efit ou eftm
%
% examples:
% [r,z,psinrz,sspr,sspi]=psinrzjet(shot,time,nrg_rg,nzg_zg,[efitlab,uid,seq,ncont]);
% [r,z,psinrz,sspr,sspi]=psinrzjet(50814,60,65,65,[],[],[],60,sspr,sspi); % to get plot and give sspr,sspi
% [r,z,psinrz]=psinrzjet(50814,60.4,[3 3.2],[0 0.1],[],[],[],0,sspr,sspi,[],1); %
%
% entrees :
% shot : numero du choc jet
% time : time de l'analyse
% nrg_rg, nzg_zg: nb de points de la grille en r (resp. en z) sur laquelle on fait la
% reconstruction des surfaces de flux.
% if nzg_zg is negative, make symmetric box around zero
% if array, assumes rout and zout given firectly
% varargin{1}: plot option: 0: do not plot contours, >0 plot contour with varargin{1} nb of contours (60 is good)
% varargin{2}: sspr from ppf/efit/sspr, structure containing sspr.data and sspr.t
% varargin{3}: sspi from ppf/efit/sspi, structure containing sspi.data and sspi.t
% varargin{4}: deltaz
% varargin{5}: idiag: =1 : gives diag(psinrz) as output to give psi at (r,z) points
% idiag=0 (default) full matrix as output
%
% facultatifs :
% efitlab : efit ou eftm (efitlab=efit par defaut)
% uid : eventuellement, donnees ppf privees (uid='jetppf' par defaut).
% seq : "sequence number" de l'uid (0 par defaut -> version la plus recente)
%
% sorties :
% r, z : vecteurs de la grille (r,z)
% psinrz : matrice psin(r,z) (flux normalise)
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if nargin<=4 | isempty(efitlab)
efitlab='efit';
end
if nargin<=5 | isempty(uid)
uid='jetppf';
end
if nargin<=6 | isempty(seqd)
seqd='0';
end
ncont=0;
if nargin>=8 & ~isempty(varargin{1})
ncont=varargin{1};
end
% equilibre magnetique
%---------------------
if nargin>=9 & ~isempty(varargin{2})
sspr=varargin{2};
else
sspr=gdat(shot,'ppf',[efitlab '/sspr?uid=' uid '+seq=' seqd]);
end
ssprs=sspr.data;
tpefit=sspr.t;
if nargin>=10 & ~isempty(varargin{3})
sspi=varargin{3};
else
sspi=gdat(shot,'ppf',[efitlab '/sspi?uid=' uid '+seq=' seqd]);
end
sspis=sspi.data;
tpefit=sspi.t;
if nargin>=11 & ~isempty(varargin{4})
deltaz=varargin{4};
else
deltaz=0;
end
if nargin>=12 & ~isempty(varargin{5})
idiag=varargin{5};
else
idiag=0;
end
[x,ind]=min(abs(time-tpefit));
sspr_t=ssprs(:,ind);
sspi_t=sspis(:,ind);
nr=sspi_t(1);
nz=sspi_t(2);
nc=sspi_t(3);
ri=sspr_t(1:nr);
zi=sspr_t((nr+1):(nr+nz));
cij=sspr_t((nr+nz+1):(nr+nz+nc));
rmin=0.01*sspr_t(1);
rmax=0.01*sspr_t(nr);
zmin=0.01*sspr_t(nr+1);
zmax=0.01*sspr_t(nr+nz);
rmin=rmin+1e-4;
rmax=rmax-1e-4;
zmin=zmin+1e-4;
zmax=zmax-1e-4;
if length(nrg_rg)==1 & length(nzg_zg)==1
r=linspace(rmin,rmax,nrg_rg);
z=linspace(zmin,zmax,abs(nzg_zg));
if nzg_zg<0
zlim=max(abs(zmin-deltaz),abs(zmax+deltaz));
z=linspace(-zlim-deltaz,zlim-deltaz,abs(nzg_zg));
end
else
r=nrg_rg;
z=nzg_zg;
end
% mapflux contruit la carte de flux psin sur (r,z)
psinrz=mapflux_fast(cij,0.01*ri,0.01*zi,r,z,idiag);
if ncont>0
figure
contour(r,z,psinrz',ncont);
hold on
[h1, h2]=contour(r,z,psinrz',[1 1],'k');
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment