//===========================ファイル名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();
        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);
   }
}