NBabel

 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