1 NBabel.java
2 Date: April 2010
3 Author: Philip Breen (University of Edinburgh)
4
5 Integration scheme: Predictor-corrector leapfrog
6 Compiler: javac
7 Operating system: Scientific Linux
8 Hardware: Intel(R) Pentium(R) 4 CPU 3.20GHz
9
10 Initial conditions: 32,64 and 128 particles Plummer distribution of equal mass
11
12 Time step (dt): constant and shared 0.001 N-body unit
13 Total time (tEnd): 1.0
14
15 Performance:
16
17 N CPU time dE/E
18 32 2.778 7.34E-7
19 64 9.882 3.60E-6
20 128 38.906 1.72E-6
21
22 **This code is submitted free of license
23
24 How to compile:
25 >javac NBable.java
26
27 How to run:
28 >java NBable < inputfile
1 import java.lang.Math;
2 import java.util.*;
3 import java.io.*;
4 import java.lang.Double;
5
6 public class NBable{
7 public static void main(String[] args)
8 {
9 int count = 0;
10 double t = 0.0;
11 double dt = 0.001;
12 double tend = 1.0;
13
14 cluster cl = new cluster();
15 //Scanner in = new Scanner(System.in);
16 BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
17 String Input = null;
18 String dummy=null;
19
20 try{
21 while((Input = br.readLine()) != null){
22
23 StringTokenizer numbers = new StringTokenizer(Input);
24 double mass;
25 double[] pos = new double[3];
26 double[] vel = new double[3];
27 while (numbers.hasMoreTokens()) {
28 dummy=numbers.nextToken();
29 mass = Double.parseDouble(numbers.nextToken());
30 pos[0] = Double.parseDouble(numbers.nextToken());
31 pos[1] = Double.parseDouble(numbers.nextToken());
32 pos[2] = Double.parseDouble(numbers.nextToken());
33 vel[0] = Double.parseDouble(numbers.nextToken());
34 vel[1] = Double.parseDouble(numbers.nextToken());
35 vel[2] = Double.parseDouble(numbers.nextToken());
36 cl.N.add(new Star(mass,pos,vel));
37 }
38
39 }
40 } catch (IOException ioe) {
41 System.out.println("IO error: check input");
42 System.exit(1);
43 }
44
45 double[] E=cl.energies();
46 double[] Eold=E;
47 System.out.println(E[0]+" "+E[1]+" "+E[2] );
48 cl.acc();
49 while(t<tend){
50 cl.update_positions(dt);
51 cl.update_acc();
52 cl.acc();
53 cl.update_velocities(dt);
54
55 count++;
56 t+=dt;
57 if(count%10==0){
58 Eold=E;
59 E=cl.energies();
60 System.out.print(" t="+t+" E="+E[0]+" KE="+E[1]+" PE="+E[2]+" de=");
61 System.out.println((E[0]-Eold[0])/Eold[0]);
62 }
63 System.out.println((0.25+E[0])/0.25);
64 }
65 }
66 }
67
68 //cluster contains Vector of star and methods to calculate acceleration, update
69 acceleration, calculate energy
70 class cluster{
71 public Vector<Star> N = new Vector<Star>();
72 public cluster(){}
73
74 public void addstar(double mass,double[] pos, double[] vel){
75 N.add(new Star(mass,pos,vel));
76 }
77
78 public void acc(){
79 for (int kk =0;kk<N.size();kk++) {
80 N.get(kk).reset_a(); //reset star acceleration (a) to zero
81 }
82
83 double rd; //rd used to store distance
84 for (int ii = 0; ii != N.size(); ii++){
85 for (int j = ii+1; j != N.size(); j++){
86 rd=0;
87 for (int i = 0; i != 3; i++){
88 rd += (N.get(ii).r[i]-N.get(j).r[i])*(N.get(ii).r[i]-N.get(j).r[
89 i]);
90 }
91 double coff=Math.pow(Math.sqrt(rd),3);
92 for (int i = 0; i != 3; i++){
93 N.get(ii).a[i] +=(N.get(j).m/coff)*(N.get(j).r[i]-N.get(ii).r[i]);
94 N.get(j).a[i] += (N.get(ii).m/coff)*(N.get(ii).r[i]-N.get(j).r[i]);
95 }
96 }
97 }
98 }
99
100 public void update_positions(double dt){ //updates the position of each star
101 for (int kk =0;kk<N.size();kk++) {
102 N.get(kk).update_position(dt);
103 }
104 }
105
106 public void update_velocities(double dt){ // Updates the velocity of
107 each star and
108 for (int kk =0;kk<N.size();kk++) { // then updata aold (Sets aold equal to
109 a)
110 N.get(kk).update_velocity(dt);
111 }
112 }
113 public void update_acc(){ // Updates aold (Sets aold equal to a)
114 for (int kk =0;kk<N.size();kk++) {
115 N.get(kk).aold_equal_a();//aold[0]=N.get(kk).a[0];
116 }
117 }
118
119
120 public double[] energies(){
121 double[] E = {0,0,0};
122 double rd =0;
123
124 //Kinetic energy
125 for (int i = 0; i != N.size(); i++){
126 for (int j = 0; j != 3; j++){
127 E[1] += N.get(i).m*(N.get(i).v[j]*N.get(i).v[j]);
128 }
129 }
130
131 //Potential energy
132 for (int ii = 0; ii != N.size()-1; ii++){
133 for (int j = ii+1; j != N.size(); j++){
134 rd=0;
135 for (int i = 0; i != 3; i++){
136 rd += (N.get(ii).r[i]-N.get(j).r[i])*(N.get(ii).r[i]-N.get(j).r[
137 i]);
138 }
139 E[2] -= (N.get(ii).m*N.get(j).m)/Math.sqrt(rd);
140 }
141 }
142 E[1]=0.5*E[1];
143 E[0] = E[1] + E[2];
144 return E;
145 }
146
147 }
148
149 /*Star class contains information about star (mass,position and velocity) and
150 methods to update position, velocity and acceleration */
151 class Star
152 {
153 public double m;
154 public double[] r = new double[3];
155 public double[] v = new double[3];
156 public double[] a = new double[3];
157 public double[] aold = new double[3];
158 public Star()
159 {
160 m = 0;
161 double[] ze= {0,0,0};
162 r=ze;
163 v=ze;
164 }
165
166 public Star(double mass,double[] pos, double[] vel)
167 {
168 m=mass;
169 r=pos;
170 v=vel;
171 }
172 public void reset_a(){
173 for (int i = 0; i != 3; i++){
174 a[i]=0.0;
175 }
176 }
177 public void aold_equal_a(){
178 for (int i = 0; i != 3; i++){
179 aold[i]=a[i];
180 }
181 }
182 public void update_position(double dt){
183 for (int i = 0; i != 3; i++){
184 double dt2=dt*dt;
185 r[i]+=v[i]*dt+0.5*a[i]*dt2;
186 }
187 }
188 public void update_velocity(double dt){
189 for (int i = 0; i != 3; i++){
190 v[i]+=0.5*dt*(a[i]+aold[i]);
191 }
192 }
193 }
194