第13章 编制和调试MATLAB程序

fdsadfdsid

%目标面高度10米，长40米（X轴步长200），宽10米（Y轴步长10）
clear;
clc;
format long;            % 长精度显示15位双精度，7位单精度
data_phy=zeros(500,100);      % 创建一个500行100列的0矩阵
E_dx=0.5;               % X轴方向b（目标面Y轴的的长）*k(步长)
E_dx_dy=0.0005;           % Y轴方向k（步长）*k(步长)100/4/(500*100)
Dphy=100;               %光通量
I0=Dphy/pi;               %中心光强：光源在指定方向的单位立体角内发出的光通量
for i=1:1            
v_N_x(1)=0;
v_N_z(1)=(v_Out_z(1)-1.4935*v_In_z(1))/sqrt(1+1.4935^2-2*1.4935*v_Out_z(1)*v_In_z(1));
x(1)=0;
z(1)=10;
x(2)=z(1)/cot(theta(1));
z(2)=10;
for i=1:99
    v_Out_x(i+1)=(200*i-x(i+1))/sqrt((200*i-x(i+1))^2+(10000-z(i+1))^2);
    v_Out_z(i+1)=(10000-z(i+1))/sqrt((200*i-x(i+1))^2+(10000-z(i+1))^2);
    v_In_x(i+1)=x(i+1)/sqrt(x(i+1)^2+z(i+1)^2);
    v_In_z(i+1)=z(i+1)/sqrt(x(i+1)^2+z(i+1)^2);
    v_N_x(i+1)=(v_Out_x(i+1)-1.4935*v_In_x(i+1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_x(i+1)*v_In_x(i+1)+v_Out_z(i+1)*v_In_z(i+1)));
    v_N_z(i+1)=(v_Out_z(i+1)-1.4935*v_In_z(i+1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_x(i+1)*v_In_x(i+1)+v_Out_z(i+1)*v_In_z(i+1)));
    x(i+2)=(v_N_x(i+1)*x(i+1)+v_N_z(i+1)*z(i+1))/(v_N_x(i+1)+v_N_z(i+1)*cot(theta(i+1)));%利用x(2)=z(1)/cot(theta(1))和向量公式
    z(i+2)=cot(theta(i+1))*x(i+2);
end
data_xyz_theta=zeros(101,3);
data_xyz_theta(:,1)=x';
data_xyz_theta(:,3)=z';
data_xyz_theta
save data_xyz_theta.txt data_xyz_theta -ascii -double;
plot3(data_xyz_theta(:,1),data_xyz_theta(:,2),data_xyz_theta(:,3),'ro');
xlabel('x');
ylabel('y');
zlabel('z');
title('四分之一自由曲面面型数据点');
hold on;
v_Out_y(1)=0;
v_Out_z(1)=1;
v_In_y(1)=0;
v_In_z(1)=1;
v_N_y(1)=0;
v_N_z(1)=(v_Out_z(1)-1.4935*v_In_z(1))/sqrt(1+1.4935^2-2*1.4935*v_Out_z(1)*v_In_z(1));
y(1)=0;
z(1)=10;
y(2)=z(1)/tan(phy(1));
z(2)=10;
for i=1:499
    v_Out_y(i+1)=(15*i-y(i+1))/sqrt((15*i-y(i+1))^2+(10000-z(i+1))^2);
    v_Out_z(i+1)=(10000-z(i+1))/sqrt((15*i-y(i+1))^2+(10000-z(i+1))^2);
    v_In_y(i+1)=y(i+1)/sqrt(y(i+1)^2+z(i+1)^2);
    v_In_z(i+1)=z(i+1)/sqrt(y(i+1)^2+z(i+1)^2);
    v_N_y(i+1)=(v_Out_y(i+1)-1.4935*v_In_y(i+1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_y(i+1)*v_In_y(i+1)+v_Out_z(i+1)*v_In_z(i+1)));
    v_N_z(i+1)=(v_Out_z(i+1)-1.4935*v_In_z(i+1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_y(i+1)*v_In_y(i+1)+v_Out_z(i+1)*v_In_z(i+1)));
    y(i+2)=(v_N_y(i+1)*y(i+1)+v_N_z(i+1)*z(i+1))/(v_N_y(i+1)+v_N_z(i+1)*tan(phy(i+1)));
    z(i+2)=tan(phy(i+1))*y(i+2);
end
data_xyz_phy=zeros(501,3);
data_xyz_phy(:,2)=y';
data_xyz_phy(:,3)=z';
data_xyz_phy
save data_xyz_phy.txt data_xyz_phy -ascii -double;
plot3(data_xyz_phy(:,1),data_xyz_phy(:,2),data_xyz_phy(:,3),'g+');
data_xyz_t_p=zeros(3,501);
for j=1:100
    for i=1:499
    v_Out_x(1)=(200*j-data_xyz_theta(j+1,1))/sqrt((200*j-data_xyz_theta(j+1,1))^2+(10000-data_xyz_theta(j+1,3))^2);
    v_Out_y(1)=0;
    v_Out_z(1)=(10000-data_xyz_theta(j+1,3))/sqrt((200*j-data_xyz_theta(j+1,1))^2+(10000-data_xyz_theta(j+1,3))^2);
    v_In_x(1)=data_xyz_theta(j+1,1)/sqrt(data_xyz_theta(j+1,1)^2+data_xyz_theta(j+1,3)^2);
    v_In_y(1)=0;
    v_In_z(1)=data_xyz_theta(j+1,3)/sqrt(data_xyz_theta(j+1,1)^2+data_xyz_theta(j+1,3)^2);
    v_N_x(1)=(v_Out_x(1)-1.4935*v_In_x(1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_x(1)*v_In_x(1)+v_Out_z(1)*v_In_z(1)));
    v_N_y(1)=0;
    v_N_z(1)=(v_Out_z(1)-1.4935*v_In_z(1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_x(1)*v_In_x(1)+v_Out_z(1)*v_In_z(1)));
    x(1)=data_xyz_theta(j+1,1);
    y(1)=0;
    z(1)=data_xyz_theta(j+1,3);
    x(2)=(v_N_x(1)*x(1)+v_N_z(1)*z(1))/(v_N_x(1)*sin(theta(j))*sin(data_phy(1,1))+v_N_z(1)*cos(theta(j))*sin(data_phy(1,1)))*sin(theta(j))*sin(data_phy(1,1));
    y(2)=(v_N_x(1)*x(1)+v_N_z(1)*z(1))/(v_N_x(1)*sin(theta(j))*sin(data_phy(1,1))+v_N_z(1)*cos(theta(j))*sin(data_phy(1,1)))*cos(data_phy(1,1));
    z(2)=(v_N_x(1)*x(1)+v_N_z(1)*z(1))/(v_N_x(1)*sin(theta(j))*sin(data_phy(1,1))+v_N_z(1)*cos(theta(j))*sin(data_phy(1,1)))*cos(theta(j))*sin(data_phy(1,1));
    v_Out_x(i+1)=(200*j-x(i+1))/sqrt((200*j-x(i+1))^2+(15*i-y(i+1))^2+(10000-z(i+1))^2);
    v_Out_y(i+1)=(15*i-y(i+1))/sqrt((200*j-x(i+1))^2+(15*i-y(i+1))^2+(10000-z(i+1))^2);
    v_Out_z(i+1)=(10000-z(i+1))/sqrt((200*j-x(i+1))^2+(15*i-y(i+1))^2+(10000-z(i+1))^2);
    v_In_x(i+1)=x(i+1)/sqrt(x(i+1)^2+y(i+1)^2+z(i+1)^2);
    v_In_y(i+1)=y(i+1)/sqrt(x(i+1)^2+y(i+1)^2+z(i+1)^2);
    v_In_z(i+1)=z(i+1)/sqrt(x(i+1)^2+y(i+1)^2+z(i+1)^2);
    v_N_x(i+1)=(v_Out_x(i+1)-1.4935*v_In_x(i+1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_x(i+1)*v_In_x(i+1)+v_Out_y(i+1)*v_In_y(i+1)+v_Out_z(i+1)*v_In_z(i+1)));
    v_N_y(i+1)=(v_Out_y(i+1)-1.4935*v_In_y(i+1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_x(i+1)*v_In_x(i+1)+v_Out_y(i+1)*v_In_y(i+1)+v_Out_z(i+1)*v_In_z(i+1)));
    v_N_z(i+1)=(v_Out_z(i+1)-1.4935*v_In_z(i+1))/sqrt(1+1.4935^2-2*1.4935*(v_Out_x(i+1)*v_In_x(i+1)+v_Out_y(i+1)*v_In_y(i+1)+v_Out_z(i+1)*v_In_z(i+1)));
    x(i+2)=(v_N_x(i+1)*x(i+1)+v_N_y(i+1)*y(i+1)+v_N_z(i+1)*z(i+1))/(v_N_x(i+1)*sin(theta(j))*sin(data_phy(i+1,1))+v_N_y(i+1)*cos(data_phy(i+1,1))+v_N_z(i+1)*cos(theta(j))*sin(data_phy(i+1,1)))*sin(theta(j))*sin(data_phy(i+1,1));
    y(i+2)=(v_N_x(i+1)*x(i+1)+v_N_y(i+1)*y(i+1)+v_N_z(i+1)*z(i+1))/(v_N_x(i+1)*sin(theta(j))*sin(data_phy(i+1,1))+v_N_y(i+1)*cos(data_phy(i+1,1))+v_N_z(i+1)*cos(theta(j))*sin(data_phy(i+1,1)))*cos(data_phy(i+1,1));
    z(i+2)=(v_N_x(i+1)*x(i+1)+v_N_y(i+1)*y(i+1)+v_N_z(i+1)*z(i+1))/(v_N_x(i+1)*sin(theta(j))*sin(data_phy(i+1,1))+v_N_y(i+1)*cos(data_phy(i+1,1))+v_N_z(i+1)*cos(theta(j))*sin(data_phy(i+1,1)))*cos(theta(j))*sin(data_phy(i+1,1));
    end
    data_xyz_t_p(1,=x;
    data_xyz_t_p(2,=y;
    data_xyz_t_p(3,=z;
    plot3(data_xyz_t_p(1,,data_xyz_t_p(2,,data_xyz_t_p(3,,'b.');
    temp=data_xyz_t_p'
    temp=data_xyz_t_p;
    str=['data_xyz_t_p(' num2str(j) ').txt'];
    fid=fopen(str,'w');
    fprintf(fid,'%f %f %f\n',temp);
    fclose(fid);
end