产生在球面上均匀分布的点

如何产生在球面上均匀分布的点呢？这里提供若干种思路。

1. 生成两个随机分布的角度，并且按照球坐标的形式产生。

缺点：产生的点是按照角度均匀分布的，但是注意这并不是球面上均匀分布的点，后面可以看到两极的地方角度明显密集。并且，由于在计算过程中有大量的三角函数的计算，程序的效率不高。

要求不高时可以采用这种方法。

思路是，采用10807 方法产生一组两个随机数，分别作为角度使用。将产生的随机数输出到plt 文件中，使用tecplot绘图。（tecplot是非常好用的CFD后处理可视化的软件，强烈推荐使用）。关于10807 方法产生随机数，有空就另外再开一篇帖子。

下面附上 C++ 代码：

 1 #include <iostream>
 2 #include<cmath>
 3 #include<fstream>
 4 
 5 using namespace std;
 6 
 7 class cRandom
 8 {
 9     public:
10         cRandom(int x,double y):seed(x),random(y){};
11         cRandom():seed(0),random(0){};
12 
13         int seed;
14         double random;
15 };
16 
17 cRandom my_random(int z)
18 // 16807 way to create random numbers
19 // z is the seed number, num is the total random number to create
20 {
21     //z(n+1)=(a*z(n)+b) mod m
22     //describe m=a*q+r to avoid that the number is large than the computer can bear
23     const int m=pow(2,31)-1;
24     const int a=16807;
25     const int q=127773;
26     const int r=2836;
27 
28     int temp=a*(z%q)-r*(z/q);
29 
30     if(temp<0)
31     {
32         temp=m+temp;
33     }
34     //z is the seed number
35     z=temp;
36     double t = z*1.0/m;
37 
38     cRandom cr;
39     cr.random=t;
40     cr.seed=z;
41 
42     return cr;
43 }
44 
45 int main()
46 {
47     cout << "Hello world!" << endl;
48     const int num = pow(10,5);
49     int z1 = 20;
50     int z2=112;
51 
52     ofstream fcout;
53     fcout.open("result.plt");
54     fcout<<" title=random  "<<endl;
55     fcout<<"  VARIABLES = \"X\",\"Y \",\"Z \" "<<endl;
56     fcout<<"zone I= "<<num<<",datapacking=POINT"<<endl;
57     cRandom sita(z1,0.0);
58     cRandom pesi(z2,0.0);
59     const double pi = 3.141592;
60 
61     for(int i=0;i!=num;++i)
62     {
63         sita=my_random(pesi.seed);
64         pesi=my_random(sita.seed);
65 
66         double x=1.0*sin(pi*sita.random)*cos(2*pi*pesi.random);
67         double y=1.0*sin(pi*sita.random)*sin(2*pi*pesi.random);
68         double z=1.0*cos(pi*sita.random);
69 
70         fcout<<x<<" "<<y<<" " <<z<<endl;
71     }
72 
73 
74     fcout.close();
75     fcout<<"this program is finished"<<endl;
76 
77     return 0;
78 }

View Code