character(len=4)::num(0:1990)
integer,parameter::n= 400000,s=50000
integer,parameter::star=3                          ! 星の数
real,parameter:: h = 0.003                         ! 時間間隔
real,parameter:: pi=3.1415926536
real :: m(star),vx(star),vy(star),vz(star),k1,k2,k3,k4,l1,l2,l3,l4
real :: a(star),b(star),c(star),next_a(star),next_b(star),next_c(star),at(star,s),bt(star,s),ct(star,s)

! 質量
m(1)=1000.0
m(2)=1000.0
m(3)=1000.0

! 初期位置
a(1)=50.0 ;b(1)=0 ;c(1)=0
a(2)=-25.0 ;b(2)=25.0*sqrt(3.0) ;c(2)=0
a(3)=-25.0 ;b(3)=-25.0*sqrt(3.0) ;c(3)=0

! 初期速度
vx(1)=0 ;vy(1)=2.0 ;vz(1)=1.0
vx(2)=-2.0*sqrt(3.0)/2.0 ;vy(2)=-2.0/2.0 ;vz(2)=-1.0
vx(3)=2.0*sqrt(3.0)/2.0 ;vy(3)=-2.0/2.0 ;vz(3)=0


open(21,file='3objs.dat')

do j = 0, n-1
   do k=1,star
      k1 = vx(k)
      l1=0
      do l=1,star
         if(k==l)cycle
         l1 = l1 - m(l)/((a(k)-a(l))**2+(b(k)-b(l))**2+(c(k)-c(l))**2)**1.5*(a(k)-a(l))
      enddo
      k2= vx(k) + l1*h/2.0
      l2=0
      do l=1,star
         if(k==l)cycle
         l2 = l2 - m(l)/((a(k)+k1*h/2.0-a(l))**2+(b(k)-b(l))**2+(c(k)-c(l))**2)**1.5*(a(k)+k1*h/2.0-a(l))
      enddo
      k3 = vx(k) + l2*h/2.0
      l3=0
      do l=1,star
         if(k==l)cycle
         l3 = l3 - m(l)/((a(k)+k2*h/2.0-a(l))**2+(b(k)-b(l))**2+(c(k)-c(l))**2)**1.5*(a(k)+k2*h/2.0-a(l))
      enddo
      k4 = vx(k) + l3*h
      l4=0
      do l=1,star
         if(k==l)cycle
         l4 = l4 - m(l)/((a(k)+k3*h-a(l))**2+(b(k)-b(l))**2+(c(k)-c(l))**2)**1.5*(a(k)+k3*h-a(l))
      enddo

      next_a(k) = a(k) + ( k1 + 2.0*k2 + 2.0*k3 + k4 ) * h/6.0
      vx(k) = vx(k) + ( l1 + 2.0*l2 + 2.0*l3 + l4 ) * h/6.0


      k1 = vy(k)
      l1=0
      do l=1,star
         if(k==l)cycle
         l1 = l1 - m(l)/((a(k)-a(l))**2+(b(k)-b(l))**2+(c(k)-c(l))**2)**1.5*(b(k)-b(l))
      enddo
      k2= vy(k) + l1*h/2.0
      l2=0
      do l=1,star
         if(k==l)cycle
         l2 = l2 - m(l)/((a(k)-a(l))**2+(b(k)+k1*h/2.0-b(l))**2+(c(k)-c(l))**2)**1.5*(b(k)+k1*h/2.0-b(l))
      enddo
      k3 = vy(k) + l2*h/2.0
      l3=0
      do l=1,star
         if(k==l)cycle
         l3 = l3 - m(l)/((a(k)-a(l))**2+(b(k)+k2*h/2.0-b(l))**2+(c(k)-c(l))**2)**1.5*(b(k)+k2*h/2.0-b(l))
      enddo
      k4 = vy(k) + l3*h
      l4=0
      do l=1,star
         if(k==l)cycle
         l4 = l4 - m(l)/((a(k)-a(l))**2+(b(k)+k3*h-b(l))**2+(c(k)-c(l))**2)**1.5*(b(k)+k3*h-b(l))
      enddo

      next_b(k) = b(k) + ( k1 + 2.0*k2 + 2.0*k3 + k4 ) * h/6.0
      vy(k) = vy(k) + ( l1 + 2.0*l2 + 2.0*l3 + l4 ) * h/6.0


      k1 = vz(k)
      l1=0
      do l=1,star
         if(k==l)cycle
         l1 = l1 - m(l)/((a(k)-a(l))**2+(b(k)-b(l))**2+(c(k)-c(l))**2)**1.5*(c(k)-c(l))
      enddo
      k2= vz(k) + l1*h/2.0
      l2=0
      do l=1,star
         if(k==l)cycle
         l2 = l2 - m(l)/((a(k)-a(l))**2+(b(k)-b(l))**2+(c(k)+k1*h/2.0-c(l))**2)**1.5*(c(k)+k1*h/2.0-c(l))
      enddo
      k3 = vz(k) + l2*h/2.0
      l3=0
      do l=1,star
         if(k==l)cycle
         l3 = l3 - m(l)/((a(k)-a(l))**2+(b(k)-b(l))**2+(c(k)+k2*h/2.0-c(l))**2)**1.5*(c(k)+k2*h/2.0-c(l))
      enddo
      k4 = vz(k) + l3*h
      l4=0
      do l=1,star
         if(k==l)cycle
         l4 = l4 - m(l)/((a(k)-a(l))**2+(b(k)-b(l))**2+(c(k)+k3*h-c(l))**2)**1.5*(c(k)+k3*h-c(l))
      enddo

      next_c(k) = c(k) + ( k1 + 2.0*k2 + 2.0*k3 + k4 ) * h/6.0
      vz(k) = vz(k) + ( l1 + 2.0*l2 + 2.0*l3 + l4 ) * h/6.0

      write(21,*)next_a(k),next_b(k),next_c(k)

   enddo

   do k=1,star
     a(k)=next_a(k)
     b(k)=next_b(k)
	 c(k)=next_c(k)
   enddo
enddo

close(21)

open(22,file='3objs.dat')

do i=1,s
   do k=1,star
      read(22,*)at(k,i),bt(k,i),ct(k,i)
   enddo

   do j=1,star*7
      read(22,*)
   enddo

enddo

close(22)


open(11,file='number')
do i=0,1990
   write(11,*)i
enddo
close(11)

open(12,file='number')
do i=0,1990
   read(12,*)num(i)
enddo
close(12)

do i=0,1990
   open(14,file=trim(num(i)))
      do j=0,49
         do k=1,star
            write(14,*)at(k,25*i+1+j),bt(k,25*i+1+j),ct(k,25*i+1+j)
         enddo
      enddo
   close(14)
enddo

end