1 NBabel.f95
2 Date: April 2010
3 Author: Philip Breen (University of Edinburgh)
4
5 Integration scheme: Predictor-corrector leapfrog
6 Compiler: gfortran
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(s) dE/E
18 32 0.189 7.339E-07
19 64 0.694 3.596E-06
20 128 2.978 1.727E-06
21
22 **This code is submitted free of license
23
24 How to compile:
25 >gfortran -o NBable NBable.f95
26 This executable NBable
27
28 How to run:
29 >./NBable
30
31 The program then will ask for the input file name and number of stars
32 (e.g input128 128)
33
34 Alternativity create a new input file by
35 >echo "input128 128" > newinput
36
37 and run by
38 >./NBable < newinput
39
1 PROGRAM NBABLE
2 Implicit none
3 INTEGER :: COUNTER=0,DUMMY,ios,i,k,N
4 Double precision :: t,dt,dt2,tend,E,Eold,PE,KE
5 Double precision, dimension(:),allocatable :: m
6 Double precision, dimension(:,:),allocatable :: X,V,A,Aold
7 Character(len=24) :: filename
8 Intrinsic :: mod
9
10 WRITE(6,'(a)',advance="no") "Enter name of input file and number of
11 stars: "
12 Read(*,*) filename,N !reads in input file name of number of stars
13 allocate(X(N,3),V(N,3),m(N)) !sets array size
14 allocate(A(N,3),Aold(N,3))
15
16 Open(unit=10,file=filename,STATUS='OLD') !opens input file and reads in data
17 Do i=1,N
18 READ(10,*) Dummy,m(i),X(i,1),X(i,2),X(i,3),V(i,1),V(i,2),V(i,3)
19 end do
20
21 t=0.0d0
22 dt=1.0d-3
23 tend=1.0d0
24 dt2=dt*dt
25 call acc()
26 call K_energy()
27 Eold=PE+KE
28 write(6,*)PE,KE
29 do while(t<tend)
30
31 X=X + dt*V+0.5d0*dt2*A !leapfrog algorithm
32 Aold=A
33 call acc()
34 V=V+ 0.5d0*dt*(Aold+A) !leapfrog algorithm
35 t=t+dt
36 counter=counter+1
37 call K_energy()
38 e=PE+KE
39 if(mod(counter,10)==0) write(6,*)" t=",t,"E=",E,"dE/E=",(E-Eold)/Eold
40 Eold=E
41 end do
42 write(6,*) "total dE/E = ",(E+0.25)/-0.25 !uncomment to get total energy
43 change
44 Contains
45 !subroutine acc() calculates acceleration and potential energy
46 subroutine acc()
47 Integer :: ii,jj,kk
48 Double precision :: distance,distance2,coff
49 Double precision, dimension(3) :: distance_vector
50 Intrinsic :: dsqrt, sum
51 A(:,:)=0.0d0
52 PE=0.0d0
53 do ii=1,N-1
54 do jj=ii+1,N
55 distance_vector=X(ii,:)-X(jj,:)
56 distance2=sum(distance_vector*distance_vector)
57 distance=DSQRT(distance2)
58 coff=distance*distance*distance
59 PE=PE-M(II)*M(JJ)/distance
60 A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
61 A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
62 end do
63 end do
64 end subroutine acc
65 !subroutine acc() kinetic energy
66 subroutine K_energy()
67 Integer :: ii
68 KE=0.0d0
69 do ii=1,N
70 KE=KE+M(II)*(v(ii,1)*v(ii,1)+v(ii,2)*v(ii,2)+v(ii,3)*v(
71 ii,3))
72 end do
73 KE=0.5d0*KE
74 end subroutine K_energy
75 END PROGRAM NBABLE
76