2

假设我在一个球体上有 10 个点(随机分布),我想旋转整个系统以确保一个点位于北极。我将如何使用 c++ 做到这一点?

我通过查看 3D 旋转矩阵来解决这个问题:

http://en.wikipedia.org/wiki/Rotation_matrix

我围绕 x 轴旋转我的点,直到 y 分量为零,然后围绕 y 轴旋转,直到 x 分量为零。这应该将问题点留在北极或南极,对吗?

我的代码如下所示:

#include <stdio.h> 
#include <string.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <time.h>
#include <stdlib.h>
#include <sstream>
using namespace std;
#define PI 3.14159265358979323846

int main()
{

  int a,b,c,f,i,j,k,m,n,s;
  double p,Time,Averagetime,Energy,energy,Distance,Length,DotProdForce,Forcemagnitude,
         ForceMagnitude[101],Force[101][4],E[1000001],En[501],x[101][4],y[101][4],
         Dist[101][101],Epsilon,z[101],theta,phi;

    /*  set the number of points */
    n=10;

    /* check that there are no more than 100 points */
    if(n>100){
      cout << n << " is too many points for me :-( \n";
      exit(0);
    }

    /* reset the random number generator */
    srand((unsigned)time(0));  

    for (i=1;i<=n;i++){
      x[i][1]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][2]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;
      x[i][3]=((rand()*1.0)/(1.0*RAND_MAX)-0.5)*2.0;

      Length=sqrt(pow(x[i][1],2)+pow(x[i][2],2)+pow(x[i][3],2));

      for (k=1;k<=3;k++){
        x[i][k]=x[i][k]/Length;
      }
    }

    /* calculate the energy */
    Energy=0.0;

    for(i=1;i<=n;i++){
      for(j=i+1;j<=n;j++){
        Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                    +pow(x[i][3]-x[j][3],2));

        Energy=Energy+1.0/Distance;
      }
    }

   cout << fixed << setprecision(10) << "energy=" << Energy << "\n";  

  /* Save Values so far */

  for(i=1;i<=n;i++){
    for(j=1;j<=3;j++){
      y[i][j]=x[i][j];
    }
  }

  /* Choose each point in turn and make it the north pole note this is what the while loop os for, but have set it to a<2 to just look at first point */

  a=1;
  b=0;
  c=0;

  while(a<2){

  /* Find theta and phi to rotate points by */

  cout << fixed << setprecision(5) << "x[" << a << "][1]=" << x[a][1] << 
  " x[" << a << "][2]=" << x[a][2] << " x[" << a << "][3]=" << x[a][3] << "\n";

  theta=x[a][2]/x[a][3];
  theta=b*PI+atan(theta);

  /* Rotate Points by theta around x axis and then by phi around y axis */

  for(i=1;i<=n;i++){
    x[i][1]=x[i][1];
    x[i][2]=x[i][2]*cos(theta)-x[i][3]*sin(theta);
    x[i][3]=x[i][2]*sin(theta)+x[i][3]*cos(theta);
  }

  phi=x[a][1]/x[a][3];
  phi=c*PI+atan(phi);

  for(i=1;i<=n;i++){
    x[i][1]=x[i][1]*cos(phi)-x[i][3]*sin(phi);
    x[i][2]=x[i][2];
    x[i][3]=x[i][1]*sin(phi)+x[i][3]*cos(phi);
  }

  cout << fixed << setprecision(5) << "x[" << a << "][1]=" << x[a][1] << 
  " x[" << a << "][2]=" << x[a][2] << " x[" << a << "][3]=" << x[a][3] << "\n";

   if(x[a][3]==1.0 && x[a][2]==x[a][3]==0)
    a=a+1;
  else if(b==0 && c==0)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        c=1;
      }
    }
  else if(b==0 && c==1)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        b=1;
        c=0;
      }
    }
  else if(b==1 && c==0)
    for(i=1;i<=n;i++){
      for(j=1;j<=3;j++){
        x[i][j]=y[i][j];
        c=1;
      }
    }
  else if(b==1 && c==1)
    break;

  }

  energy=0.0;

  for(i=1;i<=n;i++){
    for(j=i+1;j<=n;j++){

      Distance=sqrt(pow(x[i][1]-x[j][1],2)+pow(x[i][2]-x[j][2],2)
                    +pow(x[i][3]-x[j][3],2));
      energy=energy+1.0/Distance;
    }
  }  

  cout << fixed << setprecision(10) << "ENERGY=" << energy << "\n";  
  cout << fixed << setprecision(5) << "x[" << a << "][1]=" << x[a][1] << 
  " x[" << a << "][2]=" << x[a][2] << " x[" << a << "][3]=" << x[a][3] << "\n";

  /* Output to help with gnuin.txt */
  ofstream File4 ("mypoints");
  for(i=1;i<=n;i++){
    File4 << x[i][1] << " " <<   x[i][2] << " " << x[i][3] << "\n";
  }
  File4.close(); 

  return 0;

}

好的,所以这里有很多问题,比如 if 语句(第 103 行)不应该真的有等于 double 的条件,因为它永远不会起作用,但我可以稍后使用间接比较来解决这个问题(一些 epsilon 的东西)。我真正的疑问是为什么旋转即使作用于所有点也会使点脱离球体?(如您所见,这些值已被标准化以使它们都在第 38 行的单位球体上)。

注意:b、c 的东西是检查点是在北极还是南极。

4

1 回答 1

4

您的轮换代码有问题。例如:

x[i][1]=x[i][1];
x[i][2]=x[i][2]*cos(theta)-x[i][3]*sin(theta);
x[i][3]=x[i][2]*sin(theta)+x[i][3]*cos(theta);

x[i][2]在第 2 行进行修改,然后在第 3 行使用它。您应该对中间结果使用临时存储,以避免在完成引用之前修改值。

第一行是相当多余的,其余的应该看起来更像:

double new_y, new_z;
new_y=x[i][2]*cos(theta)-x[i][3]*sin(theta);
new_z=x[i][2]*sin(theta)+x[i][3]*cos(theta);
x[i][2] = new_y;
x[i][3] = new_z;

(显然在你执行这样的计算的地方这样做)

定位球体的更好方法可能是以与“查看”矩阵相同的方式计算变换矩阵。在查看矩阵中,框架被旋转以将某个向量与 z 轴对齐。在您的情况下,您可能希望沿 y 轴对齐,但原理是相同的。

我还要评论您似乎忽略了数组中的第 0 个元素......恕我直言,这是一个坏习惯 - 您应该习惯数组从零开始的事实。迟早你会得到错误的索引,或者你将不得不与其他人的代码交互。

于 2013-01-22T20:33:10.247 回答