三维射线追踪

matlab环境下，三维射线追踪，包含梯度法射线路径求取

function [PtraceT] = alxray（ xxyyzzshujumnl）
%UNtitleD3 Summary of this function goes here
%   Detailed explanation goes here
%生成走时模型


A=shuju;
% A=zeros（mnl）;%设置模型大小，充零
% for k=1:1:l
%     for i=1:1:m
%         for j=1:1:n
%             A（ijk）=sqrt（（i-1）^2+（j-1）^2+（k-1）^2）/2000;%计算解析走时
%         end
%     end
% end
% ptrace=zeros（l*53）;%最临近追踪点
% Ptrace=zeros（l*53）;%实际追踪点
Ptrace（1:）=[xx yy zz];%初始化要追踪的检波点的坐标
%%%%%%%%网格点走时贡献项初始化%%%%%%%%

T=zeros（1l*5）;
T0=zeros（l*564）;Tx0=zeros（l*564）;Ty0=zeros（l*564）;Tz0=zeros（l*564）;
Tx=zeros（1l*5）;Ty=zeros（1l*5）;Tz=zeros（1l*5）;G=zeros（l*53）;G0=zeros（l*53）;
buchang=zeros（1l*5）;%开辟存放步长的空间
WG0={l*564};WG={l*564};
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%射线追踪%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ptrace（1:）=[370 250 120];%初始化要追踪的检波点的坐标
ptrace（1:）=[ceil（Ptrace（11））-1 ceil（Ptrace（12））-1 ceil（Ptrace（13））-1];%同步初始化实际追踪点坐标
% for t=1:1:l*5 %循环次数取大一点，保证追踪到震源附近，大多数时候选的比网格点多一点就行了
for t=1:1:l*5
%     for i=1:1:64
        for ii=1:1:64
c{1}=[-1 -1 -1];c{2}=[0 -1 -1];c{3}=[1 -1 -1];c{4}=[2 -1 -1];
c{5}=[-1 0 -1];c{6}=[0 0 -1];c{7}=[1 0 -1];c{8}=[2 0 -1];
c{9}=[-1 1 -1];c{10}=[0 1 -1];c{11}=[1 1 -1];c{12}=[2 1 -1];
c{13}=[-1 2 -1];c{14}=[0 2 -1];c{15}=[1 2 -1];c{16}=[2 2 -1];
c{17}=[-1 -1 0];c{18}=[0 -1 0];c{19}=[1 -1 0];c{20}=[2 -1 0];
c{21}=[-1 0 0];c{22}=[0 0 0];c{23}=[1 0 0];c{24}=[2 0 0];
c{25}=[-1 1 0];c{26}=[0 1 0];c{27}=[1 1 0];c{28}=[2 1 0];
c{29}=[-1 2 0];c{30}=[0 2 0];c{31}=[1 2 0];c{32}=[2 2 0];
c{33}=[-1 -1 1];c{34}=[0 -1 1];c{35}=[1 -1 1];c{36}=[2 -1 1];
c{37}=[-1 0 1];c{38}=[0 0 1];c{39}=[1 0 1];c{40}=[2 0 1];
c{41}=[-1 1 1];c{42}=[0 1 1];c{43}=[1 1 1];c{44}=[2 1 1];
c{45}=[-1 2 1];c{46}=[0 2 1];c{47}=[1 2 1];c{48}=[2 2 1];
c{49}=[-1 -1 2];c{50}=[0 -1 2];c{51}=[1 -1 2];c{52}=[2 -1 2];
c{53}=[-1 0 2];c{54}=[0 0 2];c{55}=[1 0 2];c{56}=[2 0 2];
c{57}=[-1 1 2];c{58}=[0 1 2];c{59}=[1 1 2];c{60}=[2 1 2];
c{61}=[-1 2 2];c{62}=[0 2 2];c{63}=[1 2 2];c{64}=[2 2 2];
WG0{tii}=c{ii};%生成网格坐标基本cell
WG{tii}=WG0{tii}+ptrace（t:）;%生成网格坐标cell，即输出。
        end
    %%%%%%%%第一层%%%%%%%%%%
    T0（t1）=（1-（Ptrace（t1）-WG{t22}（11）））^3*（1-（Ptrace（t2）-WG{t22}（12）））^3*（1-（Ptrace（t3）-WG{t22}（13）））^3*A（WG{t1}（11）WG{t1}（12）WG{t1}（13））/（6*6*6）;
    T0（t2）=（3*（Ptrace（t1）-WG{t22}（11））^3-6*（Ptrace（t1）-WG{t22}（11））^2+4）*（1-（Ptrace（t2）-WG{t22}（12）））^3*...
        （1-（Ptrace（t3）-WG{t22}（13）））^3*A（WG{t2}（11）WG{t2}（12）WG{t2}（13））/（6*6*6）;
    T0（t3）=（-3*（Ptrace（t1）-WG{t22}（11））^3+3*（Ptrace（t1）-WG{t22}（11））^2+3*（Ptrace（t1）-WG{t22}（11））+1）*（1-（Ptrace（t2）-WG{t22}（12）））^3*...
        （1-（Ptrace（t3）-WG{t22}（13）））^3*A（WG{t3}（11）WG{t3}（12）WG{t3}（13））/（6*6*6）;
    T0（t4）=（Ptrace（t1）-WG{t22}（11））^3*（1-（Ptrace（t2）-WG{t22}（12）））^3*...
        （1-（Ptrace（t3）-WG{t22}（13）））^3*A（WG{t4}（11）WG{t4}（12）WG{t4}（13））/（6*6*6）;
    T0（t5）=（1-（Ptrace（t1）-WG{t22}（11）））^3*（3*（Ptrace（t2）-WG{t22}（12））^3-6*（Ptra

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件       68152  2017-08-07 10:23  alxray.m
     文件        1740  2017-08-06 17:02  huatu.m
     文件      354028  2017-01-09 11:31  alx.mat