c - (C 程序) 倾斜穿过环的视线 : calculating the path length along my line of sight from sources around this ring

Question

我有一个均匀的环，沿着我的视线倾斜一定角度。戒指的尺寸为:

9.5 : inner radius
10.5 : outer radius
-3 to 3 : height of the ring 

环上分布着许多源，我正在尝试计算从每个源绘制的线是否沿着我的方向穿过该环。如果它穿过环，它会表现出穿过该环的路径长度是多少(例如，穿过环两次或一次或从未穿过)。

我有从每个源到我的视线的方向 vector ，并且我使用简单的 vector 加法沿着我的方向递增。

我的程序有问题:

它仅检查路径长度是否小于环的外半径。如何检查我的源是否通过了两个环？

任何帮助将不胜感激!!

我的程序:

/********************************************************************
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.;
}

Answer 1

我在将您的文本与代码结合起来时遇到一些问题。查看代码，我看到 (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²。现在，您可以将该平方和解释为 vector v 与其自身的点积:v∙v = 10.5²。但该 vector v 实际上是位置 n 加上 u 乘以某个 vector 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 分量即可。