(C 程序)视线通过一个倾斜的环:从这个环周围的源计算沿着我的视线的路径长度
(C program) Line of sight through a ring at an inclination : calculating the path length along my line of sight from sources around this ring
我有一个均匀的戒指,沿着我的视线倾斜了一个角度。戒指的尺寸是:
9.5 : inner radius
10.5 : outer radius
-3 to 3 : height of the ring
整个环上分布着许多源,我正在尝试计算从每个源绘制的一条线是否沿着我的方向穿过这个环。如果它穿过环,它会表现出穿过这个环的路径长度(例如穿过环两次或一次或从不)。
我有从每个源朝向我的视线的方向矢量,我正在使用一个简单的矢量加法来沿我的方向递增。
我的程序有问题:
它只检查路径长度是否小于环的外半径。如何检查我的源是否通过两个环?
任何帮助将不胜感激!!
我的程序:
/********************************************************************
xn, yn, zn : coordinates of the sources
ns_ux, ns_uy, ns_uz : unit vectors for these sources
ns : distance of a source from the sun
int_val : path length
********************************************************************/
int main(){
FILE *fp=NULL;
fp=fopen("Av_path.txt","w");
int k,number=9;
float step=0.02;
float ns_ux[number],ns_uy[number],ns_uz[number],xn[number], yn[number],zn[number],l[number],b[number],ns[number],int_val,x_comp,y_comp,z_comp,radial;
FILE* val= NULL;
val=fopen("novae_uniform_unitvectors.txt", "r");
for(k=0;k<=(number-1);k++){
fscanf(val,"%f %f %f %f %f %f %f %f %f", &xn[k], &yn[k], &zn[k], &ns_ux[k], &ns_uy[k], &ns_uz[k], &l[k], &b[k], &ns[k]);
float u=0.;
do {
u=u+step;
printf("%f\t%f\n" ,u,radial);
x_comp=xn[k]+u*ns_ux[k];
y_comp=yn[k]+u*ns_uy[k];
z_comp=zn[k]+u*ns_uz[k];
radial=pow((x_comp*x_comp+y_comp*y_comp+z_comp*z_comp),0.5);
}while (radial <10.5);
if(u >= 1.0 && z_comp >= -3.0 && z_comp <= 3.0){
int_val=1.0; // ring's width is only 1 unit
}
else if(u < 1.0 && z_comp >= -3.0 && z_comp <= 3.0) {
int_val=u-step;
}
else {
int_val=0.;
}
fprintf(fp, "%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",xn[k], yn[k], zn[k], ns[k], l[k], b[k], radial,z_comp,u,u/step, int_val);
}
return 0.;
}
我在将您的文本与代码结合起来时遇到了一些问题。查看代码,我看到 (xn, yn, zn)
处的一个点和从该点向 (ns_ux, ns_uy, ns_uz)
方向发出的光线。我想你想知道光线何时何地与你的环相交。
目前你做一些非常讨厌的循环来找到一个交点。让我们让它更优雅。您正在寻找 u
使得 sqrt(x_comp²+y_comp²+z_comp²)=10.5
。这是要求光线与球体相交的点,所以如果你想要一个圆柱体,你可能想要删除三个坐标之一。我现在会坚持使用球体,但这很容易适应。
首先要对等式求两边的平方:x_comp²+y_comp²+z_comp²=10.5²
。现在您可以将该平方和解释为向量 v
与自身的点积:v∙v = 10.5²
。但是那个向量 v
实际上是位置 n
加上 u
乘以某个向量 ns_u?
,我称之为 s
,所以你有 (n + u*s)∙(n + u*s)=10.5²
.现在您可以利用分配律可以应用于点积的事实,因为它是双线性的。所以你也可以这样写
(n∙n - 10.5²) + (2n∙s)u + (s∙s)u² = 0
这是u
中的一个二次方程,可以用二次公式求解。在代码中你会有
nn = xn*xn + yn*yn + zn*zn;
ns = xn*sn_ux + yn*sn_uy + zn*sn_uz;
ss = sn_ux*sn_ux + sn_uy*sn_uy + sn_uz*sn_uz;
r = 10.5;
a = ss;
b = 2*ns;
c = nn - r*r;
d = b*b - 4*a*c;
if (d < 0) {
// doesn't intersect sphere of radius r at all, but passes outside.
} else if (d == 0) {
// tangential to sphere, numerically extremely unlikely
// so you might want to merge this case with the one below.
} else {
// two points of intersection, parametrized by u1 and u2.
sqrtd = sqrt(d);
u1 = (-b+sqrtd)/(2*a);
u2 = (-b-sqrtd)/(2*a);
}
上面的 线 与球体相交。如果你想要射线,你会把自己限制在解决方案 u≥0
上。您可能希望对较小的半径执行相同的计算。一条射线在两个点上与较大的球体相交,但与较小的球体完全不相交,它只会通过这两个球体之间的区域一次。所有这些都假设是完整的领域;您仍然需要应用该厚度约束来查看光线在环的高度之外的位置。由于我还没有完全理解您在这方面的问题,因此我无法提供更多详细信息,但我希望以上内容能够帮助您自行解决该问题。如果您想要圆柱体而不是球体,只需从上述计算中删除所有 z
个组件。
我有一个均匀的戒指,沿着我的视线倾斜了一个角度。戒指的尺寸是:
9.5 : inner radius
10.5 : outer radius
-3 to 3 : height of the ring
整个环上分布着许多源,我正在尝试计算从每个源绘制的一条线是否沿着我的方向穿过这个环。如果它穿过环,它会表现出穿过这个环的路径长度(例如穿过环两次或一次或从不)。
我有从每个源朝向我的视线的方向矢量,我正在使用一个简单的矢量加法来沿我的方向递增。
我的程序有问题:
它只检查路径长度是否小于环的外半径。如何检查我的源是否通过两个环?
任何帮助将不胜感激!!
我的程序:
/********************************************************************
xn, yn, zn : coordinates of the sources
ns_ux, ns_uy, ns_uz : unit vectors for these sources
ns : distance of a source from the sun
int_val : path length
********************************************************************/
int main(){
FILE *fp=NULL;
fp=fopen("Av_path.txt","w");
int k,number=9;
float step=0.02;
float ns_ux[number],ns_uy[number],ns_uz[number],xn[number], yn[number],zn[number],l[number],b[number],ns[number],int_val,x_comp,y_comp,z_comp,radial;
FILE* val= NULL;
val=fopen("novae_uniform_unitvectors.txt", "r");
for(k=0;k<=(number-1);k++){
fscanf(val,"%f %f %f %f %f %f %f %f %f", &xn[k], &yn[k], &zn[k], &ns_ux[k], &ns_uy[k], &ns_uz[k], &l[k], &b[k], &ns[k]);
float u=0.;
do {
u=u+step;
printf("%f\t%f\n" ,u,radial);
x_comp=xn[k]+u*ns_ux[k];
y_comp=yn[k]+u*ns_uy[k];
z_comp=zn[k]+u*ns_uz[k];
radial=pow((x_comp*x_comp+y_comp*y_comp+z_comp*z_comp),0.5);
}while (radial <10.5);
if(u >= 1.0 && z_comp >= -3.0 && z_comp <= 3.0){
int_val=1.0; // ring's width is only 1 unit
}
else if(u < 1.0 && z_comp >= -3.0 && z_comp <= 3.0) {
int_val=u-step;
}
else {
int_val=0.;
}
fprintf(fp, "%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",xn[k], yn[k], zn[k], ns[k], l[k], b[k], radial,z_comp,u,u/step, int_val);
}
return 0.;
}
我在将您的文本与代码结合起来时遇到了一些问题。查看代码,我看到 (xn, yn, zn)
处的一个点和从该点向 (ns_ux, ns_uy, ns_uz)
方向发出的光线。我想你想知道光线何时何地与你的环相交。
目前你做一些非常讨厌的循环来找到一个交点。让我们让它更优雅。您正在寻找 u
使得 sqrt(x_comp²+y_comp²+z_comp²)=10.5
。这是要求光线与球体相交的点,所以如果你想要一个圆柱体,你可能想要删除三个坐标之一。我现在会坚持使用球体,但这很容易适应。
首先要对等式求两边的平方:x_comp²+y_comp²+z_comp²=10.5²
。现在您可以将该平方和解释为向量 v
与自身的点积:v∙v = 10.5²
。但是那个向量 v
实际上是位置 n
加上 u
乘以某个向量 ns_u?
,我称之为 s
,所以你有 (n + u*s)∙(n + u*s)=10.5²
.现在您可以利用分配律可以应用于点积的事实,因为它是双线性的。所以你也可以这样写
(n∙n - 10.5²) + (2n∙s)u + (s∙s)u² = 0
这是u
中的一个二次方程,可以用二次公式求解。在代码中你会有
nn = xn*xn + yn*yn + zn*zn;
ns = xn*sn_ux + yn*sn_uy + zn*sn_uz;
ss = sn_ux*sn_ux + sn_uy*sn_uy + sn_uz*sn_uz;
r = 10.5;
a = ss;
b = 2*ns;
c = nn - r*r;
d = b*b - 4*a*c;
if (d < 0) {
// doesn't intersect sphere of radius r at all, but passes outside.
} else if (d == 0) {
// tangential to sphere, numerically extremely unlikely
// so you might want to merge this case with the one below.
} else {
// two points of intersection, parametrized by u1 and u2.
sqrtd = sqrt(d);
u1 = (-b+sqrtd)/(2*a);
u2 = (-b-sqrtd)/(2*a);
}
上面的 线 与球体相交。如果你想要射线,你会把自己限制在解决方案 u≥0
上。您可能希望对较小的半径执行相同的计算。一条射线在两个点上与较大的球体相交,但与较小的球体完全不相交,它只会通过这两个球体之间的区域一次。所有这些都假设是完整的领域;您仍然需要应用该厚度约束来查看光线在环的高度之外的位置。由于我还没有完全理解您在这方面的问题,因此我无法提供更多详细信息,但我希望以上内容能够帮助您自行解决该问题。如果您想要圆柱体而不是球体,只需从上述计算中删除所有 z
个组件。