diff --git a/JET/mapflux_fast.m b/JET/mapflux_fast.m
new file mode 100644
index 0000000000000000000000000000000000000000..4cb5fee37af39cd8f7b0d8d3795e4c929ff8a41b
--- /dev/null
+++ b/JET/mapflux_fast.m
@@ -0,0 +1,98 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+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
diff --git a/JET/psinrzjet.m b/JET/psinrzjet.m
new file mode 100644
index 0000000000000000000000000000000000000000..09c82d122b0163b3dc511893a065479526752f4e
--- /dev/null
+++ b/JET/psinrzjet.m
@@ -0,0 +1,128 @@
+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