%--CTSimpson3D積分 ,對三維火焰溫度場進行投影的數值模擬 --2005/9/16
%--p1投影值（S×K×S矩陣），K--投影方向數，S--每方向投影數（射線數），此處沿Z軸所取層數也取值S

function [p1,K,S]=CTSimpson3D(K,S)

%--da間隔角度，xra xrb  x軸起止點，yra yrb  y軸起止點，y軸間隔d，zra zrb  z軸起止點，z軸間隔dz,h積分步長，p投影向量
da=pi/K;xra=-sqrt(2)/2;xrb=sqrt(2)/2;yra=-0.5;yrb=0.5;
d=(yrb-yra)/(S-1);L=51;h=(xrb-xra)/(L-1);p=[];
zra=0;zrb=1;dz=(zrb-zra)/(S-1);

for zs=0:S-1
    zr=dz*zs;
    for k=0:K-1
        for s=0:S-1
            a=k*da;%確定射線方向
            yr=yra+s*d;%確定射線沿著yr軸的位置
            v1=[cos(a) sin(a);-sin(a) cos(a)];   %旋轉向量使積分線以da大小旋轉
            v=inv(v1);
            f1=0;f2=0;
            for l=0:L-1
                xr=xra+(2*l+1)*h;
                xy=v*[xr;yr];
                f1=f1+ThrDMod(xy(1),xy(2),zr);
            end
            for l=1:L-1
                xr=xra+2*l*h;
                xy=v*[xr;yr];
                f2=f2+ThrDMod(xy(1),xy(2),zr);
            end
            xy=v*[xra;yr];
            fa=ThrDMod(xy(1),xy(2),zr);
            xy=v*[xrb;yr];
            fb=ThrDMod(xy(1),xy(2),zr);
            p1(s+1,k+1,zs+1)=h/3*(fa+4*f1+2*f2+fb);
        end
    end
end