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

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);
r=reshape(r,1,length(r));
z=reshape(z,1,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));
  if ~isempty(ii)
    i2=ii(1);
  else
    i2=nr;
  end
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));
  if ~isempty(ii)
    i2=ii(1);
  else
    i2=nz;
  end
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