You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 

72 lines
2.7 KiB

function [ch_da]=interp_ch_da(a0,b0,c0)
load Cda_components_largerange.mat aa bb cc...
c11 c12 c13 c21 c22 c23 c31 c32 c33 c44 c55 c66
ha = aa(2)-aa(1);
hb = bb(2)-bb(1);
hc = cc(2)-cc(1);
ja=min(floor((a0-aa(1))/ha)+1,size(aa,2)-1);
jb=min(floor((b0-bb(1))/hb)+1,size(bb,2)-1);
jc=min(floor((c0-cc(1))/hc)+1,size(cc,2)-1);
ja=max(ja,1);jb=max(jb,1);jc=max(jc,1);
N=zeros(1,8);
ch_da=zeros([12,size(a0)]);
nz=size(a0,1);
ny=size(a0,2);
nx=size(a0,3);
for iz=1:nz
for iy=1:ny
for ix=1:nx
a=a0(iz,iy,ix);
b=b0(iz,iy,ix);
c=c0(iz,iy,ix);
a00=aa(ja(iz,iy,ix));
b00=bb(jb(iz,iy,ix));
c00=cc(jc(iz,iy,ix));
N(1)=(ha+a00-a)*(hb+b00-b)*(hc+c00-c)/(ha*hb*hc);
N(2)=(a-a00)*(hb+b00-b)*(hc+c00-c)/(ha*hb*hc);
N(3)=(ha+a00-a)*(b-b00)*(hc+c00-c)/(ha*hb*hc);
N(4)=(a-a00)*(b-b00)*(hc+c00-c)/(ha*hb*hc);
N(5)=(ha+a00-a)*(hb+b00-b)*(c-c00)/(ha*hb*hc);
N(6)=(a-a00)*(hb+b00-b)*(c-c00)/(ha*hb*hc);
N(7)=(ha+a00-a)*(b-b00)*(c-c00)/(ha*hb*hc);
N(8)=(a-a00)*(b-b00)*(c-c00)/(ha*hb*hc);
R=compR(c11,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(1,iz,iy,ix)=N*R;
R=compR(c12,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(2,iz,iy,ix)=N*R;
R=compR(c13,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(3,iz,iy,ix)=N*R;
R=compR(c21,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(4,iz,iy,ix)=N*R;
R=compR(c22,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(5,iz,iy,ix)=N*R;
R=compR(c23,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(6,iz,iy,ix)=N*R;
R=compR(c31,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(7,iz,iy,ix)=N*R;
R=compR(c32,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(8,iz,iy,ix)=N*R;
R=compR(c33,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(9,iz,iy,ix)=N*R;
R=compR(c44,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(10,iz,iy,ix)=N*R;
R=compR(c55,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(11,iz,iy,ix)=N*R;
R=compR(c66,ja(iz,iy,ix),jb(iz,iy,ix),jc(iz,iy,ix));
ch_da(12,iz,iy,ix)=N*R;
end
end
end
end
function R=compR(c,aj,bj,cj)
R=zeros(8,1);
R(1)=c(aj,bj,cj);
R(2)=c(aj+1,bj,cj);
R(3)=c(aj,bj+1,cj);
R(4)=c(aj+1,bj+1,cj);
R(5)=c(aj,bj,cj+1);
R(6)=c(aj+1,bj,cj+1);
R(7)=c(aj,bj+1,cj+1);
R(8)=c(aj+1,bj+1,cj+1);
end