Newer
Older
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
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