Diagnostics/DiagUtility/calc_Bolus_CS.m

function [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS);
% [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS);
%-- GM-Bolus transp. : 
%   from Kwx & Kwy (Skew-Flx => /2), 
%    compute Volume Stream function psiX,psiY above uVel.vVel 
%     (at interface between 2 levels), units=m^3/s :
%     psiX=(rAc*kwx)_i / dXc ; psiY=(rAc*kwy)_j / dYc ;
%   and then the bolus velocity (m/s):
%     ub = d_k(psiX)/dYg/drF ; vb = d_k(psiY)/dXg/drF ;
%---------------------------------------------------------------------
% set_axis
[ncx,nc]=size(rAc); nr=length(delR);
nPg=6*nc*nc; ncp=nc+1;

rAu=dXc.*dYg; rAu=reshape(rAu,nPg,1);
rAv=dYc.*dXg; rAv=reshape(rAv,nPg,1);

%-- K*rAc + add 1 overlap :
tmp=reshape(rAc,nPg,1)*ones(1,nr);
kwx=reshape(kwx,nPg,nr);
kwy=reshape(kwy,nPg,nr);
tmp=0.5*tmp;
kwx=tmp.*kwx;
kwy=tmp.*kwy;
kwx=reshape(kwx,ncx,nc,nr);
kwy=reshape(kwy,ncx,nc,nr);
v6X=split_C_cub(kwx,1);
v6Y=split_C_cub(kwy,1);
k6x=v6X(:,[2:ncp],:,:);
k6y=v6Y([2:ncp],:,:,:);

%--- recip_hFac & mask :
hFacW=reshape(hFacW,nPg,nr);
hFacS=reshape(hFacS,nPg,nr);
rhFcW=1./hFacW; rhFcW(find(hFacW==0))=0; 
rhFcS=1./hFacS; rhFcS(find(hFacS==0))=0; 
hFacW=ceil(hFacW); hFacW=min(1,hFacW);
hFacS=ceil(hFacS); hFacS=min(1,hFacS);

%----------------- 

 v6X=zeros(nc,nc,nr,6);
 v6X([1:nc],:,:,:)=k6x([2:ncp],:,:,:)+k6x([1:nc],:,:,:);
 v6X=v6X/2;
 psiX=zeros(ncx,nc,nr+1); 
for n=1:6
 is=1+nc*(n-1);ie=nc*n;
 psiX([is:ie],[1:nc],[1:nr])=v6X([1:nc],[1:nc],[1:nr],n);
end
 psiX=reshape(psiX,nPg,nr+1); 
 psiX(:,[1:nr])=hFacW.*psiX(:,[1:nr]);
 ub=zeros(nPg,nr);
 ub=psiX(:,[2:nr+1])-psiX(:,[1:nr]);
 delR = reshape(delR,[1,length(delR)]);
 tmp=rAu*delR;
 ub=ub./tmp; ub=rhFcW.*ub;
 ub=reshape(ub,ncx,nc,nr);
 
%----------------- 

 v6Y=zeros(nc,nc,nr,6);
 v6Y(:,[1:nc],:,:)=k6y(:,[2:ncp],:,:)+k6y(:,[1:nc],:,:);
 v6Y=v6Y/2;
 psiY=zeros(ncx,nc,nr+1);
for n=1:6
 is=1+nc*(n-1);ie=nc*n;
 psiY([is:ie],[1:nc],[1:nr])=v6Y([1:nc],[1:nc],[1:nr],n);
end
 psiY=reshape(psiY,nPg,nr+1); 
 psiY(:,[1:nr])=hFacS.*psiY(:,[1:nr]);
 vb=zeros(nPg,nr);
 vb=psiY(:,[2:nr+1])-psiY(:,[1:nr]);
 tmp=rAv*delR;
 vb=vb./tmp; vb=rhFcS.*vb;
 vb=reshape(vb,ncx,nc,nr);
 
%----------------- 
return
1	function [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS);
2	% [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS);
3	%-- GM-Bolus transp. :
4	% from Kwx & Kwy (Skew-Flx => /2),
5	% compute Volume Stream function psiX,psiY above uVel.vVel
6	% (at interface between 2 levels), units=m^3/s :
7	% psiX=(rAckwx)_i / dXc ; psiY=(rAckwy)_j / dYc ;
8	% and then the bolus velocity (m/s):
9	% ub = d_k(psiX)/dYg/drF ; vb = d_k(psiY)/dXg/drF ;
10	%---------------------------------------------------------------------
11	% set_axis
12	[ncx,nc]=size(rAc); nr=length(delR);
13	nPg=6ncnc; ncp=nc+1;
14
15	rAu=dXc.*dYg; rAu=reshape(rAu,nPg,1);
16	rAv=dYc.*dXg; rAv=reshape(rAv,nPg,1);
17
18	%-- K*rAc + add 1 overlap :
19	tmp=reshape(rAc,nPg,1)*ones(1,nr);
20	kwx=reshape(kwx,nPg,nr);
21	kwy=reshape(kwy,nPg,nr);
22	tmp=0.5*tmp;
23	kwx=tmp.*kwx;
24	kwy=tmp.*kwy;
25	kwx=reshape(kwx,ncx,nc,nr);
26	kwy=reshape(kwy,ncx,nc,nr);
27	v6X=split_C_cub(kwx,1);
28	v6Y=split_C_cub(kwy,1);
29	k6x=v6X(:,[2:ncp],:,:);
30	k6y=v6Y([2:ncp],:,:,:);
31
32	%--- recip_hFac & mask :
33	hFacW=reshape(hFacW,nPg,nr);
34	hFacS=reshape(hFacS,nPg,nr);
35	rhFcW=1./hFacW; rhFcW(find(hFacW==0))=0;
36	rhFcS=1./hFacS; rhFcS(find(hFacS==0))=0;
37	hFacW=ceil(hFacW); hFacW=min(1,hFacW);
38	hFacS=ceil(hFacS); hFacS=min(1,hFacS);
39
40	%-----------------
41
42	v6X=zeros(nc,nc,nr,6);
43	v6X([1:nc],:,:,:)=k6x([2:ncp],:,:,:)+k6x([1:nc],:,:,:);
44	v6X=v6X/2;
45	psiX=zeros(ncx,nc,nr+1);
46	for n=1:6
47	is=1+nc(n-1);ie=ncn;
48	psiX([is:ie],[1:nc],[1:nr])=v6X([1:nc],[1:nc],[1:nr],n);
49	end
50	psiX=reshape(psiX,nPg,nr+1);
51	psiX(:,[1:nr])=hFacW.*psiX(:,[1:nr]);
52	ub=zeros(nPg,nr);
53	ub=psiX(:,[2:nr+1])-psiX(:,[1:nr]);
54	delR = reshape(delR,[1,length(delR)]);
55	tmp=rAu*delR;
56	ub=ub./tmp; ub=rhFcW.*ub;
57	ub=reshape(ub,ncx,nc,nr);
58
59	%-----------------
60
61	v6Y=zeros(nc,nc,nr,6);
62	v6Y(:,[1:nc],:,:)=k6y(:,[2:ncp],:,:)+k6y(:,[1:nc],:,:);
63	v6Y=v6Y/2;
64	psiY=zeros(ncx,nc,nr+1);
65	for n=1:6
66	is=1+nc(n-1);ie=ncn;
67	psiY([is:ie],[1:nc],[1:nr])=v6Y([1:nc],[1:nc],[1:nr],n);
68	end
69	psiY=reshape(psiY,nPg,nr+1);
70	psiY(:,[1:nr])=hFacS.*psiY(:,[1:nr]);
71	vb=zeros(nPg,nr);
72	vb=psiY(:,[2:nr+1])-psiY(:,[1:nr]);
73	tmp=rAv*delR;
74	vb=vb./tmp; vb=rhFcS.*vb;
75	vb=reshape(vb,ncx,nc,nr);
76
77	%-----------------
78	return