//===========================ファイル名Planet===========
//４次のルンゲ-クッタ法を用いた惑星運動シミュレーションVer.0.9
// 作成：加藤徳善 1997.3.18
// 作成にあたっては長嶋さん北村さんによるpenduram6
//を参考にしています。
// ルンゲ-クッタ法はm次元n体について一般的に作ってありますが，
//このアプレットではその一部だけを用いています。
// BASICで自作したものをJavaに変換したため，変数の整理ができ
//ていません。悪しからず。
//====================================================
import java.applet.*;
import java.awt.*;
//import java.lang.*;
//import java.util.*;
public class Planet extends Applet implements Runnable
{
Graphics g_this;
Image offScreen;
Graphics offGraphics;
Dimension dim1;
Thread th;
Color col=new Color(75,85,130);
//ルンゲ-クッタ関連配列
double x[][]=new double[10][3];
double v[][]=new double[10][3];
double dx[][][]=new double[10][10][3];
double nrm[][][]=new double[10][10][3];
double pow[][]=new double[10][10];
double ik[][][]=new double[4][10][3];
double im[][][]=new double[4][10][3];
double r[][]=new double[10][10];
double aa[][]=new double[10][10];
double al[][]=new double[10][10];
int ap[][]=new int[10][10];
double b[]=new double[10];
double c[]=new double[10];
double d[][]=new double[10][3];
double x0[][]=new double[10][3];
double mass[]=new double[10];
double x1[][]=new double[10][3];
double v1[][]=new double[10][3];
double v2[][]=new double[10][3];
double vn[][]=new double[10][3];
double w[][]=new double[10][3];
double w1[][]=new double[10][3];
double va[]=new double[10];
double f[][]=new double[10][3];
double coef[]={1.,2.,2.,1.};
//力のスイッチ
boolean swa=true;//距離に関係する力（万有引力，弾性力など）
boolean swb=true;//速さに比例する力（粘性抵抗）
boolean swc=true;//速さの二乗に比例する力（慣性抵抗）
boolean swd=true;//一定の力（地上の重力など）
int dim;//次元-1
int n;//物体数-1
double dt;//時間間隔
double time;//時間
int FONTSIZE1=12;
int oldx[]=new int[3];
int oldnv[]=new int[4];
Point pt0=new Point(220,200);
int drawmode=0;//スレッドによる 1:描画 2:描画停止
int drawinitflag=1;// 1:画面初期化 0:何もしない
int stepflag=0;//1:ストロボモード 0:何もしない
int drawvflag=1;//速度ベクトルを表示する
int DD=15;
int dd=5;
public void init(){
int i,k;
Font f=new Font("TimesRoman",Font.PLAIN,FONTSIZE1);//フォント作
setFont(f);//フォント選択
setLayout(new BorderLayout());//ボーダーレイアウトに
Panel p = new Panel();//パネル作製
p.setLayout(new GridLayout(4,1));
p.add(new Button("Start"));//部品をパネルに貼り付け
p.add(new Button("Stop"));
p.add(new Button("Strobo"));
p.add(new Button("Reset"));
add("East", p);//パネルを上に貼り付け
time=0;
initrk();//ルンゲ-クッタ関係の初期化
//calc();
drawinitflag=1;
resize(500,400);
setBackground(col);
dim1=this.size();//現在の表示画面の大きさを調べる
if (offScreen == null) {
offScreen = createImage(dim1.width,dim1.height);
offGraphics = offScreen.getGraphics();
}
}
//加速度の計算
public void acce(){
int i,j,k;
double fw,fw1;
if(swa||swc){
for(i=0;i<=n;i++){
if(swa){
for(j=0;j<=i;j++){
fw=0.;
for(k=0;k<=dim;k++){
if(i!=j){
dx[i][j][k]=x1[i][k]-x1[j][k];
dx[j][i][k]=x1[j][k]-x1[i][k];
}
else
dx[i][j][k]=x1[i][k]-x0[j][k];
fw+=dx[i][j][k]*dx[i][j][k];
}
r[i][j]=Math.sqrt(fw);
for(k=0;k<=dim;k++){
if(r[i][j]!=0)
nrm[i][j][k]=dx[i][j][k]/r[i][j];
else
nrm[i][j][0]=1.;
if(i!=j)
nrm[j][i][k]=-nrm[i][j][k];
}
r[i][j]-=al[i][j];
r[j][i]=r[i][j];
if(ap[i][j]==0)
pow[i][j]=1.;
if(ap[i][j]==-2)
pow[i][j]=1./r[i][j]/r[i][j];
if(ap[i][j]==1)
pow[i][j]=r[i][j];
pow[i][j]=pow[i][j]*aa[i][j];
pow[j][i]=pow[i][j];
}
}
fw1=0;
for(k=0;k<=dim;k++)
fw1+=v1[i][k]*v1[i][k];
va[i]=Math.sqrt(fw1);
for(k=0;k<=dim;k++){
if(va[i]!=0.)
vn[i][k]=v1[i][k]/va[i];
}
}
}
for(i=0;i<=n;i++){
for(k=0;k<=dim;k++)
f[i][k]=0.;//clear
}
for(i=0;i<=n;i++){
for(k=0;k<=dim;k++){
if(swa){
for(j=0;j<=n;j++)
f[i][k]+=pow[i][j]*nrm[i][j][k];
}
if(swb)
f[i][k]+=b[i]*v1[i][k];
if(swc)
f[i][k]+=c[i]*va[i]*va[i]*vn[i][k];
if(swd)
f[i][k]+=d[i][k];
f[i][k]/=mass[i];
}
}
}
//ルンゲ-クッタ法
public void rk(){
int i,j,k,l;
for(i=0;i<=n;i++){
for(k=0;k<=dim;k++){
x1[i][k]=x[i][k];
v1[i][k]=v[i][k];
w[i][k]=0.;
w1[i][k]=0.;
}
}
acce();
for(l=0;l<=2;l++){
for(i=0;i<=n;i++){
for(k=0;k<=dim;k++){
ik[l][i][k]=dt*v1[i][k];
im[l][i][k]=dt*f[i][k];
x1[i][k]=x[i][k]+ik[l][i][k]/coef[l+1];
v1[i][k]=v[i][k]+im[l][i][k]/coef[l+1];
w[i][k]+=ik[l][i][k]*coef[l];
w1[i][k]+=im[l][i][k]*coef[l];
}
}
acce();
}
for(i=0;i<=n;i++){
for(k=0;k<=dim;k++){
ik[l][i][k]=dt*v1[i][k];
im[l][i][k]=dt*f[i][k];
x1[i][k]=x[i][k]+(w[i][k]+ik[3][i][k])/6.;
v1[i][k]=v[i][k]+(w1[i][k]+im[3][i][k])/6.;
}
}
for(i=0;i<=n;i++){
for(k=0;k<=dim;k++){
x[i][k]=x1[i][k];
v[i][k]=v1[i][k];
}
}
}
public void initrk(){
int i,j;
dim=1;//次元-1
n=0;//物体数-1
dt=.03;//計算での時間間隔
//物体の初めの位置x[i][n] iは物体番号,nは座標軸番号
// x y z
x[0][0]=120 ; x[0][1]=0. ;
//x[1][0]=0. ; x[1][1]=0. ;
//x[2][0]=0. ; x[2][1]=0. ; x[2][2]=0.;
//固定点の位置（物体ごとに力を受ける固定点を設定できる）
// x0 y0 z0
x0[0][0]=0. ; x0[0][1]=0. ;
//x0[1][0]=0. ; x0[1][1]=0. ;
//x0[2][0]=0. ; x0[2][1]=0. ; x0[2][2]=0.;
//物体の初速度v[i][n] iは物体番号,nは座標軸番号
// v x v y vz
v[0][0]=0. ; v[0][1]=-18. ;
//v[1][0]=0. ; v[1][1]=0. ; v[1][2]=0.;
//v[2][0]=0. ; v[2][1]=0. ; v[2][2]=0.;
//質量
mass[0]=1.;
//mass[1]=10000.;
//mass[2]=1.;
//---力の設定-------
swa=true;//２物体間の相互作用
//強さ
aa[0][0]=-50000;
//aa[1][0]=-50000; aa[1][1]=0.;
//aa[2][0]=100000.; aa[2][1]=0.; aa[2][2]=0.;
//力のべき（-2は万有引力やクーロン力，1はバネの弾性力）
ap[0][0]=-2;
//ap[1][0]=-2; ap[1][1]=0;
//ap[2][0]=100000.; ap[2][1]=0.; ap[2][2]=0.;
//バネの長さ
al[0][0]=0.;
//al[1][0]=0.; al[1][1]=0.;
//al[2][0]=100000.; al[2][1]=0.; al[2][2]=0.;
if(n>0){//対称化（作用・反作用の法則）
for(i=1;i<=n;i++){
for(j=0;j<i;j++){
aa[j][i]=aa[i][j];
ap[j][i]=ap[i][j];
al[j][i]=al[i][j];
}
}
}
swb=false;//速さに比例する力（粘性抵抗）
//b[0]=0.; b[1]=0.; b[2]=0.;
swc=false;//速さの二乗に比例する抵抗（慣性抵抗）
//c[0]=0.; c[1]=0.; c[2]=-0.001.;
swd=false;//一定の力（重力など）
//d[0][0]=0.; d[0][1]=9.8*mass[0];
//d[1][0]=0.; d[1][1]=9.8*mass[1];
//d[2][0]=0.; d[2][1]=9.8*mass[2];
}
public void calc(){
rk();//ルンゲ-クッタ呼び出し
}
public boolean mouseDrag(Event e,int xd,int yd){
int i=0;
int k;
repaint();
v[i][0]=(double)(xd-x[i][0]-pt0.x)/2.;
v[i][1]=(double)(yd-x[i][1]-pt0.y)/2.;
return true;
}
public boolean action(Event e, Object o) {
if(e.target instanceof Button) {
if ("Start".equals(o)) {
drawinitflag=1;
stepflag=0;
drawmode = 2;
calc();
drawmode = 1;
}
if ("Stop".equals(o)) {
stepflag=0;
drawmode = 2;
return true;
}
if ("Strobo".equals(o)) {
drawinitflag=1;
stepflag=1;
drawmode = 1;
return true;
}
if ("Reset".equals(o)) {
drawinitflag=1;
stepflag=0;
drawmode = 2;
initrk();
drawvflag=1;
repaint();
return true;
}
return true;
}
return super.handleEvent(e);
}
public void start(){
if(th==null){
th=new Thread(this);
th.start();
}
}
public void run(){
while(th!=null){
if(drawmode==1){
try{
calc();
calc();
calc();
calc();
calc();
calc();
calc();
calc();
calc();
calc();
repaint();
th.sleep(50);
}catch(InterruptedException e){
;
}
}
}
}
public void stop(){
if(th !=null){
th.stop();
th=null;
}
}
public void update(Graphics g){
paint(g);//updateを無効にし、paintで置き換え
}
public void paint(Graphics g){
int i,j,k;
if(drawinitflag==1){
offGraphics.setColor(getBackground());
offGraphics.fillRect(0,0,dim1.width,dim1.height);
//背景色で画面全体を塗りつぶす（初期化）
drawinitflag=0;
}
if(stepflag !=1){
offGraphics.setColor(getBackground());
offGraphics.fillRect(0,0,dim1.width,dim1.height);
//for(i=0;i<=n;i++){
// offGraphics.fillOval(oldx[0],oldx[1],(int)dd,(int)dd);
// offGraphics.drawLine(oldnv[0],oldnv[1],oldnv[2],oldnv[3]);
// }
}
offGraphics.setColor(getForeground());
offGraphics.setColor(Color.yellow);
offGraphics.fillOval(pt0.x-(int)(DD*.5),pt0.y-(int)(DD*.5),(int)DD,(int)DD);
for(i=0;i<=n;i++){
offGraphics.setColor(Color.cyan);
oldx[0]=pt0.x+(int)(x[i][0]-dd*.5);
oldx[1]=pt0.y+(int)(x[i][1]-dd*.5);
offGraphics.fillOval(oldx[0],oldx[1],(int)dd,(int)dd);
if(drawvflag==1){
offGraphics.setColor(Color.pink);
oldnv[0]=pt0.x+(int)(x[i][0]);
oldnv[1]=pt0.y+(int)(x[i][1]);
oldnv[2]=pt0.x+(int)(x[i][0]+v[i][0]*2.);
oldnv[3]=pt0.y+(int)(x[i][1]+v[i][1]*2.);
offGraphics.drawLine(oldnv[0],oldnv[1],oldnv[2],oldnv[3]);
}
}
g.drawImage(offScreen,0,0,this);
}
}