/**惑星形成シミュレーション ver 0.5
 * ワン・クラス・高速タイプ＋作用反作用＆エネルギー正規化＋２面バッファ
 * ＆パラメータ設定可
 */
import java.applet.Applet;
import java.awt.Graphics;
import java.awt.Graphics2D;            // rotation
import java.awt.Color;
import java.util.Random;
import java.util.StringTokenizer;
import java.awt.Image;   
import java.awt.image.BufferedImage;   // double buffering
//<applet code="PlanetarySystem5.java" width="300" height="300"></applet>
public class PlanetarySystem5 extends Applet implements Runnable{

	//汎用定数
	public double AU = 149597870000.0;            //天文単位

	//Appletに関するフィールド
	public static boolean test=true;
	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 BG_COLOR   = Color.black;     //初期背景色
	public static final Color CLEAR_COLOR= Color.white;     //消去色
	public static final Color FILL_COLOR = Color.red;       //描画色
	public Graphics2D appGraph;    //アプレット画像インスタンス
  public BufferedImage memImg; // ２面画像バッファ
  public Graphics2D memGraph;  // ２面画像バッファの描画対象
	public Thread thread;//スレッド
	public boolean clearGraphicsFlag;//画面消去用フラグ

	//再生制御に関するフィールド
	public boolean running=false; //シミュレーション開始フラグ
	public int stepRunning=0;     //ステップ再生スタック
	public int stepRunningValue;  //ステップ再生時のステップ数
	public double dt;//微少時間
	public long time;//累積時間
	
	//系全体に関するフィールド
	public double G;                         //万有引力定数
	public int minorbodies;                  //メンバ小天体数
	public boolean energyNormalizationSwitch;//エネルギー正規化スイッチ
	public int edgeMode;                     //端処理モード
	public int conflictionMode;              //天体同士衝突モード
	public int randomSeedMode;               //乱数シードモード
	public double reflectionCV;              //垂直反射係数
	public double reflectionCH;              //水平反射係数
	public Random random;					 //乱数クラス
	public long randomSeed;                  //乱数シード値
	public double msum;						 //質量合計
	public double rsum;                      //半径合計

	public double distributionSize;//天体初期分布範囲
	public double mmin;//天体最小質量
	public double mmax;//天体最大質量
	public double pmin;//天体最小運動量
	public double pmax;//天体最大運動量
	public double rmin;//天体最小半径
	public double rmax;//天体最大半径

	//天体データ配列
	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 double[] xx;//2DApplet消去用x座標
	public double[] yy;//2DApplet消去用y座標
	public double[] zz;//2DApplet消去用z座標
	public boolean[] isConflicted;//衝突フラグ


	//ウィンドウ変数
	public double scale;//表示される領域の幅(=高さ=奥行き)
	public double scale0;//表示最小座標
	public double scale1;//表示最大座標
	public int radiusMode;//描画半径モード
	public int radiusMax;//描画半径最大
	public boolean perspectiveSwitch;//遠近法スイッチ


	//Applet初期化
	public void init(){
		appGraph = (Graphics2D)this.getGraphics();//画像クラス取得
    // double buffering
    memImg   = new BufferedImage((int)APPLET_WIDTH,
                                 (int)APPLET_HEIGHT,
                                 BufferedImage.TYPE_INT_BGR);
    memGraph = memImg.createGraphics();
	}

	//Applet描画---------------------------------------------
	
	//天体描画メソッド
	public void drawMinorbodies(int i,Color color){
		if((x[i]>=scale0)&&(x[i]<=scale1)&&(y[i]>=scale0)&&(y[i]<=scale1)&&(z[i]>=scale0)&&(z[i]<=scale1)){//表示領域内なら
			double x_,y_,r_;
			x_=x[i]/scale*APPLET_WIDTH;//x
			y_=y[i]/scale*APPLET_HEIGHT;//y
			r_=radiusMax;//z
			if(radiusMode==0){r_ = r_ * r[i]/rsum;}//半径に比例
			if(radiusMode==1){r_ = r_ * m[i]/msum;}//質量に比例
			if(perspectiveSwitch){//遠近法ONなら
				r_ = r_ * (z[i]/scale+0.5);
			}

			if(r_<3){r_=3;}
			memGraph.setColor(color);
			memGraph.fillOval((int)(x_-r_/2.0)+APPLET_WIDTH/2,(int)(y_-r_/2)+APPLET_HEIGHT/2,(int)r_,(int)r_);//消去
		}
	}
	//１ステップ前天体描画メソッド
	public void drawPreviousMinorbodies(int i,Color color){
		if((xx[i]>=scale0)&&(xx[i]<=scale1)&&(yy[i]>=scale0)&&(yy[i]<=scale1)&&(zz[i]>=scale0)&&(zz[i]<=scale1)){//表示領域内なら
			double x_,y_,r_;
			x_=xx[i]/scale*APPLET_WIDTH;//x
			y_=yy[i]/scale*APPLET_HEIGHT;//y
			r_=radiusMax;//z
			if(radiusMode==0){r_ = r_ * r[i]/rsum;}//半径に比例
			if(radiusMode==1){r_ = r_ * m[i]/msum;}//質量に比例
			if(perspectiveSwitch){//遠近法ONなら
				r_ = r_ * (zz[i]/scale+0.5);
			}

			if(r_<3){r_=3;}
			memGraph.setColor(color);
			memGraph.fillOval((int)(x_-r_/2.0)+APPLET_WIDTH/2,(int)(y_-r_/2)+APPLET_HEIGHT/2,(int)r_,(int)r_);//消去
		}
	}
	//再描画
	public void paint(Graphics g){/*nop*/}
	//任意描画
	public void paintDirectly(){
		//天体描画
			memGraph.setColor(CLEAR_COLOR);
			memGraph.fillRect(0,0,APPLET_WIDTH,APPLET_HEIGHT);
		for(int i=0;i<minorbodies;i++){
			drawMinorbodies(i,FILL_COLOR);
			xx[i]=x[i];
			yy[i]=y[i];
			zz[i]=z[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 = memGraph.getFontMetrics().stringWidth(displayedString);
		int county = memGraph.getFontMetrics().getHeight();

		memGraph.setColor(FILL_COLOR);//時間・天体数表示
		memGraph.drawString(displayedString,APPLET_WIDTH-countx,APPLET_HEIGHT-5);

    appGraph.drawImage(memImg,0,0,(int)APPLET_WIDTH, (int)APPLET_HEIGHT,this);

	}
	//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=0,Er=0;     //前エネルギー・エネルギー増加率
		
	 	//衝突判定
		if(l2<(r[b]+r[a])*(r[b]+r[a])){//衝突したとき
			if(conflictionMode==1){//合体モード
				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
			}else if(conflictionMode==2){//弾性反射
			}
		}
		//エネルギー計算
		if(energyNormalizationSwitch){
			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(edgeMode==2){
			if(x[a]<scale0){x[a]=scale0;px[a]=+Math.abs(px[a])*reflectionCV; py[a]=py[a] * reflectionCH; pz[a]= pz[a] * reflectionCH;}
			if(y[a]<scale0){y[a]=scale0;py[a]=+Math.abs(py[a])*reflectionCV; pz[a]=pz[a] * reflectionCH; px[a]= px[a] * reflectionCH;}
			if(z[a]<scale0){z[a]=scale0;pz[a]=+Math.abs(pz[a])*reflectionCV; px[a]=px[a] * reflectionCH; py[a]= py[a] * reflectionCH;}
			if(x[a]>scale1){x[a]=scale1;px[a]=-Math.abs(px[a])*reflectionCV; py[a]=px[a] * reflectionCH; pz[a]= pz[a] * reflectionCH;}
			if(y[a]>scale1){y[a]=scale1;py[a]=-Math.abs(py[a])*reflectionCV; pz[a]=px[a] * reflectionCH; px[a]= px[a] * reflectionCH;}
			if(z[a]>scale1){z[a]=scale1;pz[a]=-Math.abs(pz[a])*reflectionCV; px[a]=px[a] * reflectionCH; py[a]= py[a] * reflectionCH;}
			if(x[b]<scale0){x[b]=scale0;px[b]=+Math.abs(px[b])*reflectionCV; py[b]=py[a] * reflectionCH; pz[b]= pz[b] * reflectionCH;}
			if(y[b]<scale0){y[b]=scale0;py[b]=+Math.abs(py[b])*reflectionCV; pz[b]=pz[a] * reflectionCH; px[b]= px[b] * reflectionCH;}
			if(z[b]<scale0){z[b]=scale0;pz[b]=+Math.abs(pz[b])*reflectionCV; px[b]=px[a] * reflectionCH; py[b]= py[b] * reflectionCH;}
			if(x[b]>scale1){x[b]=scale1;px[b]=-Math.abs(px[b])*reflectionCV; py[b]=py[a] * reflectionCH; pz[b]= pz[b] * reflectionCH;}
			if(y[b]>scale1){y[b]=scale1;py[b]=-Math.abs(py[b])*reflectionCV; pz[b]=pz[a] * reflectionCH; px[b]= px[b] * reflectionCH;}
			if(z[b]>scale1){z[b]=scale1;pz[b]=-Math.abs(pz[b])*reflectionCV; px[b]=px[a] * reflectionCH; py[b]= py[b] * reflectionCH;}
		}else if(edgeMode==1){
			if(x[a]<scale0 && px[a]<0){x[a]=scale1;}
			if(y[a]<scale0 && py[a]<0){y[a]=scale1;}
			if(z[a]<scale0 && pz[a]<0){z[a]=scale1;}
			if(x[a]>scale1 && px[a]>0){x[a]=scale0;}
			if(y[a]>scale1 && py[a]>0){y[a]=scale0;}
			if(z[a]>scale1 && pz[a]>0){z[a]=scale0;}
			if(x[b]<scale0 && px[b]<0){x[b]=scale1;}
			if(y[b]<scale0 && py[b]<0){y[b]=scale1;}
			if(z[b]<scale0 && pz[b]<0){z[b]=scale1;}
			if(x[b]>scale1 && px[b]>0){x[b]=scale0;}
			if(y[b]>scale1 && py[b]>0){y[b]=scale0;}
			if(z[b]>scale1 && pz[b]>0){z[b]=scale0;}			
		}
		//エネルギー再計算
		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)){
						if(j<i){i--;}
						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];
						}
						minorbodies--;//天体数減少
					}
				}
			}
		}
	}

	//スレッド-------------------------------
	public void run(){
		while(true){
			if(running||stepRunning>0){
				paintDirectly();//描画
				calc();                 //計算
				time+=dt;               //時間累積
				stepRunning--;
			}
		}
	}

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

	//シミュレーション制御-------------------
	public void startSim(){
		clearGraphicsFlag=true;
		running = true;
	}
	public void stopSim(){
		running = false;
	}

	public void StepSim(String s){
		clearGraphicsFlag=true;
		stepRunning = stepRunningValue;
	}
	//javascriptからパラメータ読み込み---------------------
	public void recieveParameters(String s){
		StringTokenizer st=new StringTokenizer(s,":");
		st.nextToken();
		stepRunningValue = Integer.parseInt(st.nextToken());
		dt = precode(st.nextToken());
		scale = precode(st.nextToken())*AU;
		radiusMax = Integer.parseInt(st.nextToken());
		perspectiveSwitch = st.nextToken().equals("1");
		radiusMode = Integer.parseInt(st.nextToken());
		edgeMode = Integer.parseInt(st.nextToken());
		reflectionCV = precode(st.nextToken());
		reflectionCH = precode(st.nextToken());
		randomSeedMode = Integer.parseInt(st.nextToken());
		G = precode(st.nextToken());
		randomSeed = Long.parseLong(st.nextToken());
		conflictionMode = Integer.parseInt(st.nextToken());
		energyNormalizationSwitch = st.nextToken().equals("1");
		minorbodies = Integer.parseInt(st.nextToken());
		distributionSize = precode(st.nextToken()) * AU;
		mmin = precode(st.nextToken());
		mmax = precode(st.nextToken());
		rmin = precode(st.nextToken());
		rmax = precode(st.nextToken());
		pmin = precode(st.nextToken());
		pmax = precode(st.nextToken());
		
		scale0= -scale/2.0;
		scale1= +scale/2.0;
		initialize();

	}

	//小天体関連配列設定＆乱数使用
	public void initialize(){
		time=0;
		if(randomSeedMode==0){
			randomSeed=System.currentTimeMillis();
		}
		random = new Random(randomSeed);

		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 double[minorbodies];
		yy=new double[minorbodies];
		zz=new double[minorbodies];
		isConflicted=new boolean[minorbodies];
		rsum=0.0;
		msum=0.0;
		for(int i=0;i<minorbodies;i++){
			x[i] = random.nextDouble() * distributionSize - distributionSize/2.0;
			y[i] = random.nextDouble() * distributionSize - distributionSize/2.0;
			z[i] = random.nextDouble() * distributionSize - distributionSize/2.0;
			px[i] = random.nextDouble() * (pmax-pmin) + pmin;
			py[i] = random.nextDouble() * (pmax-pmin) + pmin;
			pz[i] = random.nextDouble() * (pmax-pmin) + pmin;
			m[i] = random.nextDouble() * (mmax-mmin) + mmin;
			r[i] = random.nextDouble() * (rmax-rmin) + rmin;
			isConflicted[i]=false;
			rsum += r[i]*r[i]*r[i];
			msum += m[i];
		}
		rsum = Math.pow(rsum , 1.0/3.0);
	}
	public static double precode(String s){
		int i;
		while((i=s.indexOf(" "))!=-1){//空白の除去
			s=s.substring(0,i)+s.substring(i+1);
		}
		i=s.indexOf("*10^");
		if(i!=-1){//10のべき乗があれば計算
			return Double.parseDouble(s.substring(0,i)) * Math.pow(10.0,Double.parseDouble(s.substring(i+4)));
		}else{
			return Double.parseDouble(s);
		}
	}

}