実践Linux                 TOPへ  Cプログラミング目次へ

サイエンス・プログラミング  2014年11月更新

●共通Makefile
Makefileは以下の内容(プログラムソースがtest.cの場合)
字下げはTabキーで行うこと。

INCLUDE = /usr/X11R6/include/X11
LIB = /usr/X11R6/lib
CC = gcc
CFLAGS = -lX11 -lpthread -lm
TARGETS = test
OBJ = ${SRC:.c=.o}
SRC = test.c

$(TARGETS):${OBJ}
  ${CC} -I ${INCLUDE} -L ${LIB} ${CFLAGS} ${OBJ} -o ${TARGETS}

${OBJ}:${SRC}
  ${CC} ${SRC} -c -o ${OBJ}

clean:
  rm -f *.o $(TARGETS)


プログラムのダウンロード(science.tar.gz)

science.tar.gz をダウンロードして、解凍して下さい。

ダウンロードしたら、まず、readmeファイルを読んでください。

ソースのフォントの部分をシステムに合わせて入れ替えて、make しなおす必要があるかもしれません。
 #define FONT1 "-urw-nimbus sans l-medium-i-condensed--16-0-0-0-p-0-iso10646-1"
 #define FONT2 "-urw-nimbus sans l-medium-i-condensed--50-0-0-0-p-0-iso10646-1"
の部分。
フォントの入れ替え方
xlsfontsコマンドでシステムフォントを一覧。 # xlsfonts -fn "*-iso10646-1"
iso10646-1は、ユニコードを意味する。
********--0-0-0-0-p-0-iso10646-1とかの最初の0を12などに変更すると文字サイズが指定できる。
********--12-0-0-0-p-0-iso10646-1

画像が不安定なときは、描画タイミングをとっているusleepの調整をしてくでださい。(readmeファイル参照)


●放物線運動      GTK+で作成したプログラムはもっとスマートです。



空気抵抗がない場合と速度の2乗に比例する空気抵抗がある場合の2つのシミュレーションをおこないます。
端末からプログラムを実行し、引数に初速度と打ち上げ角度を指定します。 (例 ./test 75 45  初速75 角度45)

プログラムのソース  ソースファイルの#define FONT1, 2の部分をシステムに合わせたものに変更して、makeし直して使用してください。
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/keysym.h>

#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define FONT1 "-urw-nimbus sans l-medium-i-condensed--16-0-0-0-p-0-iso10646-1"
#define FONT2 "-urw-nimbus sans l-medium-i-condensed--50-0-0-0-p-0-iso10646-1"
#define G 9.8
#define PI 3.14159
#define DT 0.01 //近似計算の時間間隔
#define K 0.001 //空気抵抗係数(物によってかわってくる)  ←空気抵抗の大きさを変えたい場合はここを変更してください。

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;
double v, q; //初速と角度


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

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,"parabola");

XSelectInput(d,wr,ExposureMask | KeyPressMask);
XSelectInput(d,wb,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);

///周辺ウィンドの描画
XSetForeground(d,gc,white);
XFillRectangle(d,pr,gc,0,0,800,530);
XSetForeground(d,gc,black);
XDrawString(d,pr,gc,570,513,"quit : press here =>",20);
XCopyArea(d,pr,wr,gc,0,0,800,530,0,0);

///上ウィンドの基本描画
XSetForeground(d,gc1,white);
XFillRectangle(d,p1,gc1,0,0,780,400);
XSetForeground(d,gc1,black);
XDrawString(d,p1,gc1,200,385,"200",3);
XDrawString(d,p1,gc1,400,385,"400",3);
XDrawString(d,p1,gc1,600,385,"600",3);
XDrawLine(d,p1,gc1,0,370,780,370);
XDrawLine(d,p1,gc1,0,170,780,170);
XDrawLine(d,p1,gc1,200,0,200,370);
XDrawLine(d,p1,gc1,400,0,400,370);
XDrawLine(d,p1,gc1,600,0,600,370);
XCopyArea(d,p1,w1,gc1,0,0,780,400,0,0);

///下ウィンドの基本描画
XSetForeground(d,gc2,black);
XFillRectangle(d,p2,gc2,0,0,780,80);
XCopyArea(d,p2,w2,gc2,0,0,780,80,0,0);

///quitウィンドの描画
XSetForeground(d,gcb,white);
XSetBackground(d,gcb,col[7][0][0].pixel);
XReadBitmapFile(d,wb,"./quit.bmp",&bwidth, &bheight, &pb, &x_hot, &y_hot);
XCopyPlane(d,pb,wb,gcb,0,0,bwidth,bheight,0,0,1);
}

////////////////////////////////////////////////////////////////////////////////////////////////
void draw_t(double t, int xo, int yo, int x, int y)
{
char s[16];

XDrawLine(d,p1,gc1,xo,370-yo,x,370-y);
XCopyArea(d,p1,w1,gc1,0,0,780,400,0,0);

sprintf(s,"%3dm %5.2fsec",x,t);
XSetForeground(d,gc2,black);
XFillRectangle(d,p2,gc2,0,0,780,80);
XSetForeground(d,gc2,col[0][7][0].pixel);
XDrawString(d,p2,gc2,100,60,s,strlen(s));
XCopyArea(d,p2,w2,gc2,0,0,780,80,0,0);

usleep(1000); //この待ち時間をうまく調整しないと(短いと)エラーが出る。イベントスレッドのusleep()との兼ね合いで決める。
}

////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char *argv[])
{
pthread_t th1, th2;
void th1_func();
void th2_func();
void color();
void draw_base();

system("clear");

if(argc < 3){
fputs("usage : program-name <speed> <angle>\n",stderr);
exit(1);
}
v=atof(argv[1]); //速度
q=atof(argv[2]); //角度

d=XOpenDisplay(NULL);

color();
draw_base();

pthread_create(&th1, NULL, (void *)th1_func, NULL);
pthread_create(&th2, NULL, (void *)th2_func, NULL);

while(1);

return(0);
}

/////////////////////////////////////////////////////////////////////////////////////////
void th1_func()
{
XEvent event;
char ch[8];
KeySym key;

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

case ButtonPress:
if(event.xbutton.window == wb)
{
exit(EXIT_FAILURE);
}

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

case XK_F2:
exit(EXIT_FAILURE);
}
}
usleep(1000); //これを入れるとdraw_t()のusleep()を短くできる。draw_t()と同じぐらいの長さにするとよい。
}
}

///////////////////////////////////////////////////////////////////
void th2_func()
{
double vx0, vy0;
double t=0, x=0, y=0, xo, yo, vx, vy, vxo, vyo;
void draw_t(double t, int xo, int yo, int x, int y);

vx0=v*cos(q*PI/180);
vy0=v*sin(q*PI/180);

XSetForeground(d,gc1,col[7][0][0].pixel);
while(1)
{
xo=x;
yo=y;
t=t+DT;
x=vx0*t;
y=vy0*t-G*t*t/2;
if(y<0)
{
x=vx0*2*vy0/G;
y=0;
t=2*vy0/G;
draw_t(t,(int)(xo),(int)(yo),(int)(x),(int)(y));
XFlush(d);
break;
}

draw_t(t,(int)(xo),(int)(yo),(int)(x),(int)(y));
XFlush(d); //これがないと動きがぎくしゃくする。
}

sleep(2);

XSetForeground(d,gc1,col[0][0][7].pixel);
t=0; x=0; y=0; vx=vx0; vy=vy0;
while(1)
{
xo=x;
yo=y;
vxo=vx;
vyo=vy;
vx=vxo-K*v*vxo*DT;
vy=vyo-(G+K*v*vyo)*DT;
x=xo+(vxo-K*v*vxo*DT/2)*DT;
y=yo+(vyo-(G+K*v*vyo)*DT/2)*DT;
v=sqrt(vx*vx+vy*vy);
t=t+DT;
if(y<0)
{
break;
}

draw_t(t,(int)(xo),(int)(yo),(int)(x),(int)(y));
XFlush(d); //これがないと動きがぎくしゃくする
}
}


●天体運動(3体問題)



このプログラムは、計算上解けない3体以上の天体運動をシミュレーションするものです。ただし、x−y軸上での表示にしてあります。
dataファイルにそれぞれの天体の初期値を入れておきます。300:-50:0:0:0:5:0 となっていれば、質量300、位置(-50,0,0)、初速(0,5,0)となります。最初の3体までは色つきで表示されますが、それ以上のものはすべて黒となります。また、最初の3体までは下の欄にもデータが表示されるようになっています(画面をクリックすると原点が変更になるとともに、それぞれの現在値が表示されます)。
デフォルトでは3体なので、これを変更したい場合は、ソースの#define COUNT 3の部分を変更してmakeし直してください。さらにdataファイルを書き換えるのも忘れないように!
近似計算には、独自のフィードバック方式と3点補間(2次式)による平均法を使っています。なお、天体同士が真心で衝突してしまうような異常事態では計算誤差は当然極端になってしまいます。


プログラムのソース  ソースファイルの#define FONT1, 2の部分をシステムに合わせたものに変更して、makeし直して使用してください。
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/keysym.h> //functionKey

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

#define FONT1 "-adobe-utopia-bold-r-normal--14-0-0-0-p-0-iso10646-1"
#define FONT2 "-adobe-utopia-bold-r-normal--12-0-0-0-p-0-iso10646-1"
#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];

//////////////////////////////////////////////////////////////////////////////////////////////////
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_w1base(unsigned long c);
void draw_w2();

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,"show planets");

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

///周辺ウィンドの描画
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);
XCopyArea(d,pr,wr,gc,0,0,800,530,0,0);

///上ウィンドの基本描画 線は黒で
draw_w1base(black);
XCopyArea(d,p1,w1,gc1,0,0,780,400,0,0);

///下ウィンドの描画
draw_w2();
XCopyArea(d,p2,w2,gc2,0,0,780,80,0,0);

///quitウィンドの描画
XSetForeground(d,gcb,white);
XSetBackground(d,gcb,col[7][0][0].pixel);
XReadBitmapFile(d,wb,"./quit.bmp",&bwidth, &bheight, &pb, &x_hot, &y_hot);
XCopyPlane(d,pb,wb,gcb,0,0,bwidth,bheight,0,0,1);
}

////////////////////////////////////////////上ウィンドの基本描画//////////////////////////////////////////////////////
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);
}

///////////////////////////////////////////下ウィンドの描画///////////////////////////////////////////////////////
void draw_w2()
{
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));
}

/////////////////////////////////////////////planetの描画///////////////////////////////////////////////////

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

}
}

////////////////////////////////////////////スレッド1(イベント)////////////////////////////////////////////////////
void th1_func()
{
XEvent event;
char ch[8];
KeySym key;
void draw_w2();

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

case ButtonPress:
if(event.xbutton.window == wb) //wbウィンドつまりquitボタンを押したとき終了
{
exit(EXIT_FAILURE);
}
if(event.xbutton.window == w1) //クリックした場所に中心をずらす
{
xx0=event.xbutton.x-390+xx0;
yy0=event.xbutton.y-200+yy0;
draw_w2();
}

case KeyPress: //F1かF2キーを押したとき終了
if(XLookupString((XKeyEvent *)&event,ch,8,&key,NULL) == 1);
switch(key)
{
case XK_F1:
exit(EXIT_FAILURE);

case XK_F2:
exit(EXIT_FAILURE);
}
}
usleep(6000); //これを入れるとth2_func()のusleep()を短くできる。th2_func()のより少し長めにするとよい。
}
}

////////////////////////////////////////////スレッド2////////////////////////////////////////////////////
void th2_func()
{
void draw_w1base(unsigned long c);
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;

while(1)
{
calc(); //計算 

//////計算の後処理
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); //上ウィンドの基本描画 線は黒で

//planetの色設定と描画
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)); //planetの描画
}
XFlush(d); //これがないとスムーズに動かない。
XCopyArea(d,p1,w1,gc1,0,0,780,400,0,0);

usleep(5000); //これをうまく調整しないと途中でクラッシュする。イベントスレッドのusleep()との兼ね合いに注意。
}
}

/////////////////////////////////////////////main//////////////////////////////////////////////////////
int main()
{
pthread_t th1, th2;
void th1_func();
void th2_func();
void color();
void draw_base();

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

system("clear"); //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);

//エネルギーと運動量の初期値計算
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;
}

//描画開始 ////////////////////////////////////
d=XOpenDisplay(NULL);
color(); //カラーの設定

draw_base(); //基本描画
XFlush(d);

//スレッド ////////////////////////////////////
pthread_create(&th1, NULL, (void *)th1_func, NULL);
pthread_create(&th2, NULL, (void *)th2_func, NULL);

while(1);

return(0);
}



TOPへ  Cプログラミング目次へ