センリュウ Linuxプログラミング

＞お楽しみコーナー＞Linux＞科学プログラミング

最終更新日：2005/07/22

センリュウのホームページへようこそ！

科学プログラム２

　ここでは、簡単な科学プログラムを作ってみます。TurboLinux7ｗｏｒｋｓｔａｔｉｏｎ用に作ってありますので、それ以外で使用する場合は、ソースファイルの#define FONT1, 2 ……の部分をシステムに合わせたものに変更してmakeし直して使用してください（Redhat7.2はこのままで大丈夫。FedoraCore2もＯＫ）。

天体運動（３体問題）
　「われに初期値を与えよ。さればすべてを解いて見せよう」。ラプラスのデーモンというわけですが、じつはこんなことは現在の科学をもってしてもありえません。３体問題ですでに破産してしまうのです。人類が数式的に解けるのは２体運動までで、３体では、同一平面上にある場合など特殊な場合でしか解くことができません。さらにそれ以上となるとまったく無理となります。

　このプログラムは、計算上解けない３体以上の天体運動をシミュレーションするものです。ただし、ｘ－ｙ軸上での表示にしてあります。
　dataファイルにそれぞれの天体の初期値を入れておきます。300:-50:0:0:0:5:0 となっていれば、質量300、位置(-50,0,0)、初速(0,5,0)となります。最初の３体までは色つきで表示されますが、それ以上のものはすべて黒となります。また、最初の３体までは下の欄にもデータが表示されるようになっています（画面をクリックすると原点が変更になるとともに、それぞれの現在値が表示されます）。
　デフォルトでは３体なので、これを変更したい場合は、ソースの#define COUNT 3の部分を変更してmakeし直してください。さらにdataファイルを書き換えるのも忘れないように！

　近似計算には、独自のフィードバック方式と３点補間（２次式）による平均法を使っています。なお、天体同士が真心で衝突してしまうような異常事態では計算誤差は当然極端になってしまいます。

ダウンロード

プログラムのソース

//３体問題(増やしたいときは #define COUNT 3 を変更し、dataファイルも追記しておく。)
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/keysym.h> //functionKey

#include <pthread.h>
#include <stdlib.h>
#include <string.h>
#include <signal.h>
#include <math.h> //sin(), cos()
#include <stdio.h> //FILE

#define FONT1 "-adobe-courier-bold-r-normal-*-14-*-*-*-*-*-*-*"
#define FONT2 "-adobe-courier-bold-r-normal-*-12-*-*-*-*-*-*-*"
#define G 0.67 //6.67x10<-11>だが便宜上変えてある
#define DT 0.05 //近似計算の基本時間間隔(これが描画タイミングとなる)
#define N 200 //実際の近似計算の時間間隔はさらに細かくなり DT/N となる
#define BUFSIZE 1024
#define COUNT 3 ←天体の個数に合わせて変更してください。

Display *d;
Window wr, w1, w2, wb;
GC gc, gc1, gc2, gcb;
Pixmap pr, p1, p2, pb;
Font f1,f2;
unsigned long black,white;
XColor col[8][8][8];
int bwidth, bheight;
int xx0=0, yy0=0;
double t=0;
double px=0, py=0, pz=0, en, eno;

typedef struct {
double m;
double x;
double y;
double z;
double xo;
double yo;
double zo;
double vx;
double vy;
double vz;
double vxo;
double vyo;
double vzo;
double vxoo;
double vyoo;
double vzoo;
double ax;
double ay;
double az;
double axo;
double ayo;
double azo;
double axoo;
double ayoo;
double azoo;
} pl_element;

pl_element pl[COUNT];

int main()
{
pthread_t th1;
void th_func();
void color();
void draw_base();
void draw_w1base(unsigned long c);
void sig_fnc();
void draw_t(int r, int x, int y, int z);
void calc();

FILE *fp;
char buf[BUFSIZE+1];
int i, j;
double pld, en0=0;

system("clear");

//初期値データの読み込み
if((fp=fopen("data","r")) == NULL){
printf("can not open data file!\n");
exit(-1);
}
for(i=0;i<COUNT;i++){
if(fgets(buf,BUFSIZE,fp) != NULL){
pl[i].m=atof(strtok(buf,":"));
pl[i].xo=atof(strtok(NULL,":"));
pl[i].yo=atof(strtok(NULL,":"));
pl[i].zo=atof(strtok(NULL,":"));
pl[i].vxo=atof(strtok(NULL,":"));
pl[i].vyo=atof(strtok(NULL,":"));
pl[i].vzo=atof(strtok(NULL,":"));
}
}
fclose(fp);

d=XOpenDisplay(NULL);
color();

signal(SIGINT,sig_fnc);

pthread_create(&th1, NULL, (void *)th_func, NULL);

//エネルギーと運動量の初期値計算
for(i=0;i<COUNT;i++){
for(j=0;j<COUNT;j++){
if(j != i){
pld=sqrt(pow((pl[j].xo-pl[i].xo),2)+pow((pl[j].yo-pl[i].yo),2)+pow((pl[j].zo-pl[i].zo),2));
en0=en0-G*pl[i].m*pl[j].m/(pld*2);
}
}
px=px+pl[i].vxo*pl[i].m;
py=py+pl[i].vyo*pl[i].m;
pz=pz+pl[i].vzo*pl[i].m;
en0=en0+pl[i].m*(pow(pl[i].vxo,2)+pow(pl[i].vyo,2)+pow(pl[i].vzo,2))/2;
en=en0; eno=en0;
}

draw_base();

while(1)
{
calc(); //時間間隔DT毎の計算

//////計算の後処理
eno=en;
px=0; py=0; pz=0; en=0;
for(i=0;i<COUNT;i++){
for(j=0;j<COUNT;j++){
if(j != i){
pld=sqrt(pow((pl[j].x-pl[i].x),2)+pow((pl[j].y-pl[i].y),2)+pow((pl[j].z-pl[i].z),2));
en=en-G*pl[i].m*pl[j].m/(pld*2);
}
}
en=en+pl[i].m*(pow(pl[i].vx,2)+pow(pl[i].vy,2)+pow(pl[i].vz,2))/2;
px=px+pl[i].vx*pl[i].m;
py=py+pl[i].vy*pl[i].m;
pz=pz+pl[i].vz*pl[i].m;
}

draw_w1base(black);
for(i=0;i<COUNT;i++){
if(i == 0) XSetForeground(d,gc1,col[7][0][0].pixel);
else if(i == 1) XSetForeground(d,gc1,col[0][5][0].pixel);
else if(i == 2) XSetForeground(d,gc1,col[0][0][5].pixel);
else XSetForeground(d,gc1,col[0][0][0].pixel);
draw_t((int)(cbrt(pl[i].m)), (int)(pl[i].x), (int)(pl[i].y), (int)(pl[i].z));
}
XCopyArea(d,p1,w1,gc1,0,0,780,400,0,0);
XFlush(d);

usleep(10000); //これがないとthreadのイベント類が動かない
}

return(0);
}

//////////////////////////////////////////////////////////////////////////////////////////////////

void sig_fnc()
{
exit(0);
}

//////////////////////////////////////////////////////////////////////////////////////////////////

void color()
{
Colormap cmap;
int l,m,n;

black=BlackPixel(d,0);
white=WhitePixel(d,0);

cmap=DefaultColormap(d,0);
for(l=0;l<8;l++)
{
for(m=0;m<8;m++)
{
for(n=0;n<8;n++)
{
col[l][m][n].red=65535-(7-l)*65535/7;
col[l][m][n].green=65535-(7-m)*65535/7;
col[l][m][n].blue=65535-(7-n)*65535/7;
XAllocColor(d,cmap,&col[l][m][n]);
}
}
}
}

//////////////////////////////////////////////////////////////////////////////////////////////////

void draw_base()
{
int x_hot, y_hot;
void draw_w2base();

wr=XCreateSimpleWindow(d,RootWindow(d,0),100,50,800,530,1,black,white);
w1=XCreateSimpleWindow(d,wr,10,10,780,400,1,black,white);
w2=XCreateSimpleWindow(d,wr,10,410,780,80,1,black,black);
wb=XCreateSimpleWindow(d,wr,700,500,50,15,3,col[7][5][5].pixel,white);

XStoreName(d,wr,"test");

XSelectInput(d,wr,ExposureMask | KeyPressMask);
XSelectInput(d,wb,ButtonPressMask);
XSelectInput(d,w1,ButtonPressMask);

XMapWindow(d,wr);
XMapSubwindows(d,wr);

pr=XCreatePixmap(d,wr,800,530,DefaultDepth(d,0));
p1=XCreatePixmap(d,w1,780,400,DefaultDepth(d,0));
p2=XCreatePixmap(d,w2,780,80,DefaultDepth(d,0));
pb=XCreatePixmap(d,wb,50,15,DefaultDepth(d,0));

gc=XCreateGC(d,pr,0,0);
gc1=XCreateGC(d,p1,0,0);
gc2=XCreateGC(d,p2,0,0);
gcb=XCreateGC(d,pb,0,0);

f1=XLoadFont(d,FONT1);
f2=XLoadFont(d,FONT2);

XSetFont(d,gc,f1);
XSetFont(d,gc1,f1);
XSetFont(d,gc2,f2);

draw_w2base();

XSetForeground(d,gc,white);
XFillRectangle(d,pr,gc,0,0,800,530);
XSetForeground(d,gc,black);
XDrawString(d,pr,gc,500,515,"quit : press here =>",20);

XSetForeground(d,gcb,white);
XSetBackground(d,gcb,col[7][0][0].pixel);
XReadBitmapFile(d,wb,"./quit.bmp",&bwidth, &bheight, &pb, &x_hot, &y_hot);
}

//////////////////////////////////////////////////////////////////////////////////////////////////

void draw_w1base(unsigned long c)
{
XSetForeground(d,gc1,white);
XFillRectangle(d,p1,gc1,0,0,780,400);
XSetForeground(d,gc1,c);
XDrawLine(d,p1,gc1,0,200,780,200);
XDrawLine(d,p1,gc1,390,0,390,400);
//XCopyArea(d,p1,w1,gc1,0,0,780,400,0,0);
}

//////////////////////////////////////////////////////////////////////////////////////////////////

void draw_w2base()
{
char s[128];

XSetForeground(d,gc2,black);
XFillRectangle(d,p2,gc2,0,0,780,80);
XSetForeground(d,gc2,col[0][7][0].pixel);
sprintf(s,"1 m = %d p = (%f, %f, %f) v = (%f, %f, %f)",(int)(pl[0].m),pl[0].xo,pl[0].yo,pl[0].zo,pl[0].vxo,pl[0].vyo,pl[0].vzo);
XDrawString(d,p2,gc2,20,20,s,strlen(s));
sprintf(s,"2 m = %d p = (%f, %f, %f) v = (%f, %f, %f)",(int)(pl[1].m),pl[1].xo,pl[1].yo,pl[1].zo,pl[1].vxo,pl[1].vyo,pl[1].vzo);
XDrawString(d,p2,gc2,20,35,s,strlen(s));
sprintf(s,"3 m = %d p = (%f, %f, %f) v = (%f, %f, %f)",(int)(pl[2].m),pl[2].xo,pl[2].yo,pl[2].zo,pl[2].vxo,pl[2].vyo,pl[2].vzo);
XDrawString(d,p2,gc2,20,50,s,strlen(s));
sprintf(s,"xx0 = %d yy0 = %d t = %f MV = (%f, %f, %f) E = %f",xx0,-yy0,t,px,py,pz,en);
XDrawString(d,p2,gc2,20,65,s,strlen(s));
sprintf(s,"E-EO = %f",en-eno);
XDrawString(d,p2,gc2,20,80,s,strlen(s));
XCopyArea(d,p2,w2,gc2,0,0,780,80,0,0);
}

////////////////////////////////////////////////////////////////////////////////////////////////

void th_func()
{
XEvent event;
char ch[8];
KeySym key;
void draw_w2base();

while(1)
{
XNextEvent(d,&event);
switch(event.type)
{
case Expose:
XCopyArea(d,pr,wr,gc,0,0,800,530,0,0);
XCopyArea(d,p1,w1,gc1,0,0,780,400,0,0);
XCopyArea(d,p2,w2,gc2,0,0,780,80,0,0);
XCopyPlane(d,pb,wb,gcb,0,0,bwidth,bheight,0,0,1);
XFlush(d);

case ButtonPress:
if(event.xbutton.window == wb)
{
exit(EXIT_FAILURE);
}
if(event.xbutton.window == w1)
{
xx0=event.xbutton.x-390+xx0;
yy0=event.xbutton.y-200+yy0;
draw_w2base();
}

case KeyPress:
if(XLookupString(&event,ch,8,&key,NULL) == 1);
switch(key)
{
case XK_F1:
exit(EXIT_FAILURE);

case XK_F2:
exit(EXIT_FAILURE);
}
}
}
}

////////////////////////////////////////////////////////////////////////////////////////////////

void draw_t(int r, int x, int y, int z)
{
XFillArc(d,p1,gc1,390+x-xx0-r,200-y-yy0-r,2*r,2*r,0,360*64);
}

////////////////////////////////////////////////////////////////////////////////////////////////

void calc()
{
int i, j, k;
double pld;
double ddt = DT/N;

for(k=1;k<=N;k++){
//////仮計算
for(i=0;i<COUNT;i++){
pl[i].axo=0; pl[i].ayo=0; pl[i].azo=0;
for(j=0;j<COUNT;j++){
if(j != i){
pld=pow((pow((pl[j].xo-pl[i].xo),2)+pow((pl[j].yo-pl[i].yo),2)+pow((pl[j].zo-pl[i].zo),2)),1.5);
pl[i].axo=pl[i].axo+G*pl[j].m*(pl[j].xo-pl[i].xo)/pld;
pl[i].ayo=pl[i].ayo+G*pl[j].m*(pl[j].yo-pl[i].yo)/pld;
pl[i].azo=pl[i].azo+G*pl[j].m*(pl[j].zo-pl[i].zo)/pld;
}
}

pl[i].vx=pl[i].vxo+pl[i].axo*ddt;
pl[i].vy=pl[i].vyo+pl[i].ayo*ddt;
pl[i].vz=pl[i].vzo+pl[i].azo*ddt;

if(t == 0){
pl[i].x=pl[i].xo+(pl[i].vxo+pl[i].vx)*ddt/2;
pl[i].y=pl[i].yo+(pl[i].vyo+pl[i].vy)*ddt/2;
pl[i].z=pl[i].zo+(pl[i].vzo+pl[i].vz)*ddt/2;
}
else{
pl[i].x=pl[i].xo+(5*pl[i].vx+8*pl[i].vxo-pl[i].vxoo)*ddt/12;
pl[i].y=pl[i].yo+(5*pl[i].vy+8*pl[i].vyo-pl[i].vyoo)*ddt/12;
pl[i].z=pl[i].zo+(5*pl[i].vz+8*pl[i].vzo-pl[i].vzoo)*ddt/12;
}
}

//仮計算の結果による、ひとつ先の加速度の算出
for(i=0;i<COUNT;i++){
pl[i].ax=0; pl[i].ay=0; pl[i].az=0;
for(j=0;j<COUNT;j++){
if(j != i){
pld=pow((pow((pl[j].x-pl[i].x),2)+pow((pl[j].y-pl[i].y),2)+pow((pl[j].z-pl[i].z),2)),1.5);
pl[i].ax=pl[i].ax+G*pl[j].m*(pl[j].x-pl[i].x)/pld;
pl[i].ay=pl[i].ay+G*pl[j].m*(pl[j].y-pl[i].y)/pld;
pl[i].az=pl[i].az+G*pl[j].m*(pl[j].z-pl[i].z)/pld;
}
}
}

//////以上の結果を使った本計算
for(i=0;i<COUNT;i++){
if(t == 0){
pl[i].vx=pl[i].vxo+(pl[i].axo+pl[i].ax)*ddt/2;
pl[i].vy=pl[i].vyo+(pl[i].ayo+pl[i].ay)*ddt/2;
pl[i].vz=pl[i].vzo+(pl[i].azo+pl[i].az)*ddt/2;

pl[i].x=pl[i].xo+(pl[i].vxo+pl[i].vx)*ddt/2;
pl[i].y=pl[i].yo+(pl[i].vyo+pl[i].vy)*ddt/2;
pl[i].z=pl[i].zo+(pl[i].vzo+pl[i].vz)*ddt/2;
}
else{
pl[i].vx=pl[i].vxo+(5*pl[i].ax+8*pl[i].axo-pl[i].axoo)*ddt/12;
pl[i].vy=pl[i].vyo+(5*pl[i].ay+8*pl[i].ayo-pl[i].ayoo)*ddt/12;
pl[i].vz=pl[i].vzo+(5*pl[i].az+8*pl[i].azo-pl[i].azoo)*ddt/12;

pl[i].x=pl[i].xo+(5*pl[i].vx+8*pl[i].vxo-pl[i].vxoo)*ddt/12;
pl[i].y=pl[i].yo+(5*pl[i].vy+8*pl[i].vyo-pl[i].vyoo)*ddt/12;
pl[i].z=pl[i].zo+(5*pl[i].vz+8*pl[i].vzo-pl[i].vzoo)*ddt/12;
}
}

//////計算の後処理
t=t+ddt;
for(i=0;i<COUNT;i++){
pl[i].xo=pl[i].x; pl[i].yo=pl[i].y; pl[i].zo=pl[i].z;
pl[i].vxoo=pl[i].vxo; pl[i].vyoo=pl[i].vyo; pl[i].vzoo=pl[i].vzo;
pl[i].vxo=pl[i].vx; pl[i].vyo=pl[i].vy; pl[i].vzo=pl[i].vz;
pl[i].axoo=pl[i].axo; pl[i].ayoo=pl[i].ayo; pl[i].azoo=pl[i].azo;
}

}
}