Main

package nObjectProblem2;
 
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
 
public class Main {
	public static void main(String[] args) throws IOException{
		Comet comet=new Comet();
 
 
		/**ユリウス日によるシミュレーション開始時刻[day]*/
		double start=2443550;
		/**積分間隔[day]*/
		double dt=1;
		/**表示間隔*/
		double step=1;
		/**シミュレーション終了時刻[day]*/
		double end=2443650;
		/*彗星の軌道要素*/
		/**近日点通過時[day]*/
		comet.t=2443542;
		/**近日点距離[au]*/
		comet.q=1.441685;
		/**離心率*/
		comet.e=0.599545;
		/**近地点黄経[deg]*/
		comet.w=328.9867;
		/**昇交点黄経[deg]*/
		comet.o=121.5246;
		/**軌道傾斜角[deg]*/
		comet.i=17.8559;
 
 
 
 
		/**天体*/
		Planet planet[]=new Planet[8];
		/**天文単位[m]*/
		final double AU=1.49597870691*Math.pow(10,11);
		/**万有引力定数[AU^3/day^2/kg]*/
		final double G=6.67259*Math.pow(10,-11)*3600*24*3600*24/AU/AU/AU;
		/**地軸の傾き[degree]*/
		double eps=23.43792638;
		/**現在時刻[day]*/
		double now=start;
		BufferedReader br=new BufferedReader(new FileReader("data.csv"));
		for(int i=0;i<8;i++){
			planet[i]=new Planet();
			String str=br.readLine();
			String strs[]=str.split(",");
			planet[i].l[0]=Double.parseDouble(strs[0]);
			planet[i].l[1]=Double.parseDouble(strs[1]);
			planet[i].l[2]=Double.parseDouble(strs[2]);
			planet[i].w[0]=Double.parseDouble(strs[3]);
			planet[i].w[1]=Double.parseDouble(strs[4]);
			planet[i].w[2]=Double.parseDouble(strs[5]);
			planet[i].o[0]=Double.parseDouble(strs[6]);
			planet[i].o[1]=Double.parseDouble(strs[7]);
			planet[i].o[2]=Double.parseDouble(strs[8]);
			planet[i].i[0]=Double.parseDouble(strs[9]);
			planet[i].i[1]=Double.parseDouble(strs[10]);
			planet[i].e[0]=Double.parseDouble(strs[11]);
			planet[i].e[1]=Double.parseDouble(strs[12]);
			planet[i].a=Double.parseDouble(strs[13]);
			planet[i].m=Double.parseDouble(strs[14])*5.974*Math.pow(10,24);
		}
		br.close();
 
		PrintWriter pw;
		pw=new PrintWriter(new BufferedWriter(new FileWriter("output.csv")));
 
		comet.p=comet.start(now,eps);
		double pn[]=new double[3];
		pn=comet.start(now+1,eps);
		comet.v[0]=(pn[0]-comet.p[0])/dt;
		comet.v[1]=(pn[1]-comet.p[1])/dt;
		comet.v[2]=(pn[2]-comet.p[2]-pn[2])/dt;
 
		double a[]=new double[3];
		a[0]=0;
		a[1]=0;
		a[2]=0;
		while(now<end){
 
			for(int i=0;i<8;i++){
				planet[i].p=planet[i].nowposv(now);
				double p=Math.sqrt(Math.pow(comet.p[0]-planet[i].p[0],2)+Math.pow(comet.p[1]-planet[i].p[1],2)+Math.pow(comet.p[2]-planet[i].p[2],2));
				//a=加速度
				a[0]=a[0]-G*planet[i].m*(comet.p[0]-planet[i].p[0])/Math.pow(p,3);
				a[1]=a[1]-G*planet[i].m*(comet.p[1]-planet[i].p[1])/Math.pow(p,3);
				a[2]=a[2]-G*planet[i].m*(comet.p[2]-planet[i].p[2])/Math.pow(p,3);
			}
			comet.dv[0]=dt*a[0];
			comet.dv[1]=dt*a[1];
			comet.dv[2]=dt*a[2];
			comet.v[0]=comet.v[0]+comet.dv[0];
			comet.v[1]=comet.v[1]+comet.dv[1];
			comet.v[2]=comet.v[2]+comet.dv[2];
			comet.p[0]=comet.p[0]+comet.v[0]*dt;
			comet.p[1]=comet.p[1]+comet.v[1]*dt;
			comet.p[2]=comet.p[2]+comet.v[2]*dt;
 
			double se[]=new double[6];
			se[0]=comet.p[0]-planet[3].p[0];
			se[1]=comet.p[1]-planet[3].p[1];
			se[2]=comet.p[2]-planet[3].p[2];
			se[3]=Math.sqrt(se[0]*se[0]+se[1]*se[1]+se[2]*se[2]);
			se[4]=Math.toDegrees(Math.asin(se[2]/se[3]));//delta
			se[5]=Math.toDegrees(Math.atan(se[1]/se[0]));//alpha
			if(se[1]>0&&se[0]<0){
				se[5]=180-se[5];
			}else if(se[1]<0&&se[0]>0){
				se[5]=360-se[5];
			}else if(se[1]<0&&se[0]<0){
				se[5]=180+se[5];
			}
			if(se[5]>360){
				while(true){
					se[5]=se[5]-360;
					if(se[5]-360<0){
						break;
					}
				}
			}
 
 
			if(now==(int)(now/step)*step){
				Comet c=comet;
				pw.println(now+","+c.p[0]+","+c.p[1]+","+c.p[2]);
				//pw.print(now+","+c.p[0]+","+c.p[1]+","+c.p[2]+",");
				/*for(int i=0;i<8;i++){
					pw.print(now+","+planet[i].p[0]+","+planet[i].p[1]+","+planet[i].p[2]+",");
				}
				pw.println();*/
				System.out.println(now+",α="+(int)(se[5]/15)+"h"+(se[5]-(int)(se[5]/15)*15)*4+",δ="+se[4]+",r="+Math.sqrt(comet.p[0]*comet.p[0]+comet.p[1]*comet.p[1]+comet.p[2]*comet.p[2])+",Δ="+se[3]);
			}
			now++;
		}
 
		pw.close();
 
	}
}