clear;
clc;
format long;
%新建目标区域
data_phy=zeros(500,100);
E_dx=0.5;
E_dx_dy=0.0005;
Dphy=100;
I0=Dphy/pi;
%下述法新建的theta角为经度角,经度角可以进行平均划分 0.01 - pi/2度
for i=1:100
theta(i)=asin(2*i*E_dx/Dphy);
end
theta=theta'
%下述的phi为纬度角,此处由于光源强度分布未知,故未知其准确性,但其依然使用迭代方法计算phi0~phiN,其中N=499,表示这个方向分成500个区域进行计算
for j=1:100
for i=1:499
phy(1)=pi/2-j*E_dx_dy/(I0*sin(theta(j)));
phy(i+1)=phy(i)-j*E_dx_dy/(I0*sin(theta(j))*sin(phy(i))^2);
end
%纬度角phi序列重复(W/N)次并转置形成500x(W/N)矩阵
data_phy(:,j)=phy';
end
%初始化光线及自由曲面的法线单位向量
v_Out_x(1)=0;
v_Out_z(1)=1;
v_In_x(1)=0;
v_In_z(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)=x(1)+z(1)/cot(theta(1));
z(2)=10;
%初始化照明区域
H=7000;
L=3500;
W=3500;
M=length(theta);
N=length(phy);
for i=1:99
%设置目标照明位置距离光源z=10000,范围为(L/M)00,共分成(W/N)个小区域(theta有(W/N)份)进行能量对应,故每个区域(L/M),分别计算单位向量Vout,x
%Vout,z Vin,x Vin,z 及自由曲面的法向量
v_Out_x(i+1)=((L/M)*i-x(i+1))/sqrt(((L/M)*i-x(i+1))^2+(10000-z(i+1))^2);
v_Out_z(i+1)=(10000-z(i+1))/sqrt(((L/M)*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)));
z(i+2)=cot(theta(i+1))*x(i+2);
%最后迭代出theta所在面的所有线点
end
%在XZ面绘制种子线
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;
%计算yz面种子线
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
%设置目标照明位置距离光源z=10000,范围为5000,共分成500个小区域(phy有500份)进行能量对应,故每个区域10,分别计算单位向量Vout,y
%Vout,z Vin,y Vin,z 及自由曲面的法向量
v_Out_y(i+1)=(10*i-y(i+1))/sqrt((10*i-y(i+1))^2+(10000-z(i+1))^2);
v_Out_z(i+1)=(10000-z(i+1))/sqrt((10*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
%重复xz面种子线求解过程
v_Out_x(1)=((L/M)*j-data_xyz_theta(j+1,1))/sqrt(((L/M)*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(((L/M)*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);
%YZ面种子线作为基准,以XZ面种子线各点作为步长,求解种子线
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)=((L/M)*j-x(i+1))/sqrt(((L/M)*j-x(i+1))^2+(10*i-y(i+1))^2+(10000-z(i+1))^2);
v_Out_y(i+1)=(10*i-y(i+1))/sqrt(((L/M)*j-x(i+1))^2+(10*i-y(i+1))^2+(10000-z(i+1))^2);
v_Out_z(i+1)=(10000-z(i+1))/sqrt(((L/M)*j-x(i+1))^2+(10*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);
%
if j<20
str1=['data_xyz_pc.txt'];
fid1=fopen(str1,'a');
fprintf(fid1,'%f %f %f\n',temp)
fclose(fid1);
else if j<66
str2=['data_xyz_pc 2.txt'];
fid2=fopen(str2,'a');
fprintf(fid2,'%f %f %f\n',temp)
fclose(fid2);
%
else
str3=['data_xyz_pc 3.txt'];
fid3=fopen(str3,'a');
fprintf(fid3,'%f %f %f\n',temp)
fclose(fid3);
end
end
end