2013-01-22 88 views
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)真的不應該有平等的條件雙重因爲它會從來沒有工作,但我可以稍後使用間接比較(一些epsilon的東西)。我真正的疑問是爲什麼輪轉,即使它在所有點上作用,都會將球體上的點取下來? (正如你所看到的,這些值已經過標準化,使它們全部在第38行的單位球面上)。

注意:b,c的東西是檢查點是在北極還是南極。

+0

在這一點上,你的問題實際上是低估的:自由度不受約束導致整個類的最終狀態。這可能或可能不是一個問題。 – dmckee

+0

您必須添加一些約束條件。只需選取一個隨機點並沿着與該點相交的平面旋轉整個系統點,即球體的中心和球體頂部的虛擬點。 – Blender

回答

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); 

您在第2行修改x[i][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個元素......恕我直言,這是一個壞習慣 - 你應該習慣數組從零開始的事實。遲早你會得到錯誤的索引,或者你將不得不接口到別人的代碼。

+0

是的陣列的東西是一個壞習慣。代碼絕不美觀或體面! – adrem7