function [r,z,psinrz,sspr,sspi]=psirz(shot,time,nrg,nzg,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]=psirz(shot,time,nrg,nzg,[efitlab,uid,seq,ncont]);
%  [r,z,psinrz,sspr,sspi]=psirz(50814,60,65,65,[],[],[],60,sspr,sspi); % to get plot and give sspr,sspi
%
% entrees :
% shot : numero du choc jet
% time   : time de l'analyse
% nrg, nzg: nb de points de la grille en r (resp. en z) sur laquelle on fait la 
%           reconstruction des surfaces de flux.
%           if nzg is negative, make symmetric box around zero
% 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
%
% 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
  
[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;

r=linspace(rmin,rmax,nrg);
z=linspace(zmin,zmax,abs(nzg));
if nzg<0
  zlim=max(abs(zmin-deltaz),abs(zmax+deltaz));
  z=linspace(-zlim-deltaz,zlim-deltaz,abs(nzg));
end

% mapflux contruit la carte de flux psin sur (r,z)
psinrz=mapflux(cij,0.01*ri,0.01*zi,r,z);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function flux=mapflux(cij,ri,zi,r,z)

% 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)
%
% 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);
fpq=zeros(length(r),length(z));
for i=1:(nr-4)
   for j=1:(nz-4)
      fpq=fpq+cij(j+(i-1)*(nz-4))*bsplinev(ri,r,3,i)'*bsplinev(zi,z,3,j);
   end
end
flux=fpq;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 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