import java.applet.Applet;
import java.awt.Graphics;
import java.awt.Color;
 
/**惑星形成シミュレーション ver 0.3 
 * ワン・クラス・高速タイプ＋作用反作用＆エネルギー正規化
 */
//<applet code="PlanetarySystem3.java" width="300" height="300"></applet>
public class PlanetarySystem3 extends Applet implements Runnable{

	//汎用定数
	public double G = 6.672*Math.pow(10,-11); //万有引力定数
	public double mAU = 149597870;            //天文単位/1000

	//Applet用定数など
	public static final int APPLET_WIDTH =300;//アプレットウィンドウ幅
	public static final int APPLET_HEIGHT=300;//アプレットウィンドウ高さ
	public static final int APPLET_DEPTH =300;//アプレットウィンドウ仮想奥行き
	public static final Color bgColor   = Color.black;     //初期背景色
	public static final Color clearColor= Color.white;     //消去色
	public static final Color fillColor = Color.red;       //描画色
	public Thread thread;//スレッド
	public Graphics graphics;//画像
	
	//観測する系のフィールド
	public int minorbodies;                  //メンバ小天体数
	public boolean conflictionSwitch;        //衝突イベントスイッチ
	public boolean reflectionSwitch;         //壁反射イベントスイッチ
	public boolean energyNormalizationSwitch;//エネルギー正規化スイッチ

	//観測者のフィールド
	public double dt;                   //微少時間
	public long time;                   //累積時間
	public double RANDOM_PATTERNS=10000;//乱数パターン数

	//天体データ配列
	public double[] x;	//x座標
	public double[] y;	//y座標
	public double[] z;	//z座標
	public double[] px;	//運動量x成分
	public double[] py;	//運動量y成分
	public double[] pz;	//運動量z成分
	public double[] m;	//質量
	public double[] r;	//半径
	public int[] xx;	//2DApplet表示用x座標
	public int[] yy;	//2DApplet表示用y座標
	public int[] zz;	//2DApplet表示用z座標
	public boolean[] isSeen;//2DApplet表示フラグ
	public boolean[] isConflicted;//衝突フラグ


	//ウィンドウ変数
	public double scale;//表示される領域の幅(=高さ=奥行き)
	public double CX;//表示領域センターx座標
	public double CY;//表示領域センターy座標
	public double CZ;//表示領域センターz座標
	public double WX0;//表示領域始点x座標
	public double WY0;//表示領域始点y座標
	public double WZ0;//表示領域始点z座標
	public double WX1;//表示領域終点x座標
	public double WY1;//表示領域終点y座標
	public double WZ1;//表示領域終点z座標

	//Applet初期化
	public void init(){
		time=0;
		//コンピュータ時間によるMath.random()空回し
		for(int i=0;i<(int)(System.currentTimeMillis()%RANDOM_PATTERNS);i++){
			Math.random();
		}

		initialize();                 //フィールドの初期化
		initWindowCtoW();             //ウィンドウ初期化
		graphics = this.getGraphics();//画像クラス取得

	}

	//------------------INITIALIZATION--------------------------------------------
	public void initialize(){
		CX=0;//宇宙のセンターをフォーカス
		CY=0;//宇宙のセンターをフォーカス
		CZ=0;//宇宙のセンターをフォーカス

		scale=10*mAU;//１画面＝10ミリ天文単位
		dt=5000;     //微少時間
		conflictionSwitch         = true;//衝突スイッチ＝ON
		reflectionSwitch          = true;//壁反射スイッチ＝ON
		energyNormalizationSwitch = false;//エネルギー正規化スイッチ=OFF

		minorbodies = 300;//天体数300個
		x=new double[minorbodies];
		y=new double[minorbodies];
		z=new double[minorbodies];
		px=new double[minorbodies];
		py=new double[minorbodies];
		pz=new double[minorbodies];
		m=new double[minorbodies];
		r=new double[minorbodies];
		xx=new int[minorbodies];
		yy=new int[minorbodies];
		zz=new int[minorbodies];
		isSeen=new boolean[minorbodies];
		isConflicted=new boolean[minorbodies];
		
		for(int i=0;i<minorbodies;i++){
			//0.001天文単位^3に一様分布
			x[i]=Math.random()*10*mAU-10*mAU/2.0;
			y[i]=Math.random()*10*mAU-10*mAU/2.0;
			z[i]=Math.random()*10*mAU-10*mAU/2.0;
			//とりあえず初期速度０m/s
			px[i]=0.0;py[i]=0.0;pz[i]=0.0;
			//質量＝とりあえず小惑星マチルドの100〜200倍の値
			m[i]=(Math.random()+1.0)*1.033*Math.pow(10.0,20.0)*100;
			//半径＝とりあえず小惑星マチルドの500000〜1000000倍の値
			r[i]=(Math.random()+1.0)*9.368*500000;
			isConflicted[i]=false;
		}
	}
	
	//Applet描画---------------------------------------------
	//再描画
	public void paint(Graphics g){/*nop*/}
	//任意描画
	public void paintDirectly(Graphics g){
		for(int i=0;i<minorbodies;i++){

			//2d変換
			double xx_= (x[i]-WX0)/scale   ;//一時的にdouble
			double yy_= (y[i]-WY0)/scale   ;//一時的にdouble
			double zz_= (z[i]-WZ0)/scale   ;//一時的にdouble

			if(xx_>0.0 && yy_>0.0 && xx_<1.0 && yy_<1.0){//画面内なら
				
				g.setColor(clearColor);
				g.fillOval(xx[i]-(int)zz[i]/2,yy[i]-(int)zz[i]/2,(int)zz[i],(int)zz[i]);//消去

				xx[i]=(int)(xx_*APPLET_WIDTH );   //2d変換
				yy[i]=(int)(yy_*APPLET_HEIGHT);   //2d変換
				zz[i]=(int)(zz_*APPLET_DEPTH )/20;//2d変換
				if(zz_<0){
					zz[i]=1;
				}

				g.setColor(fillColor);
				g.fillOval(xx[i]-(int)zz[i]/2,yy[i]-(int)zz[i]/2,(int)zz[i],(int)zz[i]);//描画
			}
		}
		String displayedString = longToString((time         )/31557600,0) + "Y "//時間・天体数文字列計算
		                        +longToString((time%31557600)/   86400,3) + "D "
		                        +longToString((time%   86400)/    3600,2) + "°"
		                        +longToString((time%    3600)/      60,2) + "\' "
		                        +longToString((time%      60)/       1,2) + "\"  天体数:"
		                        +Integer.toString(minorbodies);

		//時間・天体数表示座標決定
		int countx = g.getFontMetrics().stringWidth(displayedString);
		int county = g.getFontMetrics().getHeight();

		g.setColor(clearColor);//時間・天体数消去
		g.fillRect(300-countx,300-county-5,countx,county+5);
		g.setColor(fillColor);//時間・天体数表示
		g.drawString(displayedString,300-countx,300-5);
	}
	//longのString化(前に0をつける) len=文字数
	public String longToString(long n,int len){
		if(len==0){
			return Long.toString(n);
		}else{
			String s = "0000" + Long.toString(n);
			return s.substring(s.length()-len,s.length());
		}
	}

	/**小天体メソッド*///------------------------------------

	/**aのbによる万有引力による運動量変化（dtによる近似）及び衝突
	   形式   :   beEffected(相手の小天体クラス)
   	   戻り値 : = 衝突したか                */
	public boolean beEffectedAndMove(int a,int b){
		double dxx=x[b]-x[a];//２天体距離x成分
		double dyy=y[b]-y[a];//２天体距離y成分
		double dzz=z[b]-z[a];//２天体距離z成分
		double l2= dxx*dxx + dyy*dyy + dzz*dzz;//２乗距離
		double dpx,dpy,dpz;                    //運動量微小変化
		boolean ans=false;                     //戻り値
		double E0,Er;     //前エネルギー・エネルギー増加率
		
	 	//衝突判定
	 	if(conflictionSwitch){
			if(l2<(r[b]+r[a])*(r[b]+r[a])){//衝突したとき
				px[a] = px[a]+px[b];
				py[a] = py[a]+py[b];
				pz[a] = pz[a]+pz[b];
				x[a] = (x[a]*m[a]+x[b]*m[b])/(m[a]+m[b]);//重心の位置を新天体の位置に
				y[a] = (y[a]*m[a]+y[b]*m[b])/(m[a]+m[b]);//重心の位置を新天体の位置に
				z[a] = (y[a]*m[a]+z[b]*m[b])/(m[a]+m[b]);//重心の位置を新天体の位置に
				m[a] = m[a]+m[b];
				r[a] = Math.pow(r[a]*r[a]*r[a]+r[b]*r[b]*r[b] , 1.0/3.0);//体積から計算
				isConflicted[b] = true;//相手の衝突フラグを立てる
				return true;//---return
			}
		}
		//エネルギー計算
		E0=getEnergy(a,b,l2);

		//運動量変化
		dpx = G*m[a]*m[b]/l2 * dxx/Math.abs(dxx) * dt/2.0;
		dpy = G*m[a]*m[b]/l2 * dyy/Math.abs(dyy) * dt/2.0;
		dpz = G*m[a]*m[b]/l2 * dzz/Math.abs(dzz) * dt/2.0;
		px[a] += dpx; py[a] += dpy; pz[a] += dpz;
		px[b] -= dpx; py[b] -= dpy; pz[b] -= dpz;
		//天体移動
		x[a]+=px[a]/m[a]*dt/2.0;
		y[a]+=py[a]/m[a]*dt/2.0;
		z[a]+=pz[a]/m[a]*dt/2.0;
		x[b]+=px[b]/m[b]*dt/2.0;
		y[b]+=py[b]/m[b]*dt/2.0;
		z[b]+=pz[b]/m[b]*dt/2.0;

		//壁反射
		if(reflectionSwitch){
			if(x[a]<WX0){x[a]=WX0;px[a]= Math.abs(px[a])*0;}
			if(y[a]<WX0){y[a]=WY0;py[a]= Math.abs(py[a])*0;}
			if(z[a]<WX0){z[a]=WZ0;pz[a]= Math.abs(pz[a])*0;}
			if(x[a]>WX1){x[a]=WX1;px[a]=-Math.abs(px[a])*0;}
			if(y[a]>WX1){y[a]=WY1;py[a]=-Math.abs(py[a])*0;}
			if(z[a]>WX1){z[a]=WZ1;pz[a]=-Math.abs(pz[a])*0;}
			if(x[b]<WX0){x[b]=WX0;px[b]= Math.abs(px[b])*0;}
			if(y[b]<WX0){y[b]=WY0;py[b]= Math.abs(py[b])*0;}
			if(z[b]<WX0){z[b]=WZ0;pz[b]= Math.abs(pz[b])*0;}
			if(x[b]>WX1){x[b]=WX1;px[b]=-Math.abs(px[b])*0;}
			if(y[b]>WX1){y[b]=WY1;py[b]=-Math.abs(py[b])*0;}
			if(z[b]>WX1){z[b]=WZ1;pz[b]=-Math.abs(pz[b])*0;}
		}

		//エネルギー再計算
		if(energyNormalizationSwitch){
			l2 = (x[b]-x[a])*(x[b]-x[a]) + (y[b]-y[a])*(y[b]-y[a]) + (z[b]-z[a])*(z[b]-z[a]);
			Er = Math.abs(getEnergy(a,b,l2)/E0);
			//正規化
			px[a]=px[a]/Math.sqrt(Er*0.5);
			py[a]=py[a]/Math.sqrt(Er*0.5);
			pz[a]=pz[a]/Math.sqrt(Er*0.5);
			px[b]=px[b]/Math.sqrt(Er*0.5);
			py[b]=py[b]/Math.sqrt(Er*0.5);
			pz[b]=pz[b]/Math.sqrt(Er*0.5);
		}
		return false;
	}

	//エネルギー計算
	public double getEnergy(int a,int b,double l2){
		return (px[a]*px[a] + py[a]*py[a] + pz[a]*pz[a])/(2*m[a])//運動エネルギー
		      +(px[b]*px[b] + py[b]*py[b] + pz[b]*pz[b])/(2*m[b])//運動エネルギー
		      -G*m[a]*m[b]/Math.sqrt(l2);                        //万有引力ポテンシャル
	}
	
	//計算-----------------------------------------------
	public void calc(){

		//運動量変化・衝突計算
		for(int i=0;i<minorbodies;i++){
			for(int j=0;j<minorbodies;j++){
				//自分を除外
				if(i!=j){
					//計算
					if(beEffectedAndMove(i,j)){
						//衝突が起こったら天体[j]を消す
						for(int k=j;k<minorbodies-1;k++){

							 x[k]= x[k+1]; y[k]=y[k+1]; z[k]=z[k+1];
							px[k]=px[k+1];py[k]=py[k+1];pz[k]=pz[k+1];
							 m[k]= m[k+1]; r[k]=r[k+1];
							 xx[k]= xx[k+1]; yy[k]=yy[k+1];zz[k]=zz[k+1];
							isConflicted[k]=isConflicted[k+1];
							isSeen[k]=isSeen[k+1];

						}
						minorbodies--;//天体数減少
					}
				}
			}
		}
	}

	//スレッド-------------------------------
	public void run(){
		graphics.setColor(bgColor);//背景塗り
		graphics.drawRect(0,0,APPLET_WIDTH,APPLET_HEIGHT);
		while(true){
			paintDirectly(graphics);//描画
			calc();                 //計算
			time+=dt;               //時間累積
		}
	}

	public void start(){
		if(thread == null){
			thread = new Thread(this);
			thread.start();
		}
	}

	//ウィンドウ変数の初期化------------------
	public void initWindowCtoW(){
		WX0=CX-scale/2.0;
		WY0=CY-scale/2.0;
		WZ0=CZ-scale/2.0;
		WX1=CX+scale/2.0;
		WY1=CY+scale/2.0;
		WZ1=CZ+scale/2.0;
	}

	public void initWindowWtoC(){
		if((WX1-WX0)!=(WY1-WY0)||(WX1-WX0)!=(WZ1-WZ0)){//正方形でなければWY1を正方形に正規化
			WY1=WY0+(WX1-WX0);
			WZ1=WZ0+(WX1-WX0);
		}
		scale=WX0-WX1;
		CX=(WX0+WX1)/2.0;
		CY=(WY0+WY1)/2.0;
		CZ=(WZ0+WZ1)/2.0;
	}

}