은하 시뮬레이션

 

 

 

program galaxy

implicit none

double precision vx(10000,2000),vy(10000,2000) !+++++++++++++

double precision  x(10000,2000),y(10000,2000) !+++++++++++

double precision vx1(10000,2000), vy1(10000,2000) !++++++++++ 

double precision m(2000) !==== maass

integer i,j,k,q !!!!!!-------------step variable

double precision dt !time step

double precision g !gravitioal constant

double precision r,dx,dy, dx1, dy1 ! radius, sub varible

double precision ra,ka!random varible

double precision ro !====

call random_seed()



g=50

dt=0.01



do i=1,2000,1

call random_number(ra)


call random_number(ka)


m(i)=(1-(1/4000)*ra)*10

x(1,i)=300*exp(abs(ra))*sin(ra*i)

y(1,i)=300*exp(abs(ra))*cos(ra*i)

vx(1,i)=sin(ka*i)*ra*50

vy(1,i)=cos(ka*i)*ra*50


end do

do i=1,500,1

call random_number(ra)


call random_number(ka)


x(1,i)=30*exp(abs(ra))*sin(ra*i)

y(1,i)=30*exp(abs(ra))*cos(ra*i)

vx(1,i)=sin(ka*i)*ra*30

vy(1,i)=cos(ka*i)*ra*30



end do




!++++++++++++++++++++++++++++++++++++++++++++++++++++


!++++++++++++++++++++++++++++++++++++++++++++++++++++++
do i=1,500,1

dx=0
dy=0
dx1=0
dy1=0

do j=1,2000,1

vx1(i,j)=vx(i,j)+dx*dt/2
vy1(i,j)=vy(i,j)+dy*dt/2

x(i+1,j)=x(i,j)+vx1(i,j)*dt
y(i+1,j)=y(i,j)+vy1(i,j)*dt

vx(i+1,j)=vx1(i,j)+dx1
vy(i+1,j)=vy1(i,j)+dy1

do k=1,2000,1


!!++++++++++++++++++++++++++++++++++++++++++
if(k<j .and. k>j) then
dx=dx-dt*g*m(k)*(x(i,j)-x(i,k))/(((x(i,j)-x(i,k))**2+(y(i,j)-y(i,k))**2+0.001)**(3/2))
dy=dy-dt*g*m(k)*(y(i,j)-y(i,k))/(((x(i,j)-x(i,k))**2+(y(i,j)-y(i,k))**2+0.001)**(3/2))

else if(k==i) then
dx=dx
dy=dy
end if


!+++++++++++++++++++++++++++++++++++++++++++++++++
if(k<j .and. k>j)then
dx=dx-dt*m(k)*g*(x(i+1,j)-x(i+1,k))/(((x(i+1,j)-x(i+1,k))**2+(y(i+1,j)-y(i+1,k))**2+0.001)**(3/2))
dy=dy-dt*m(k)*g*(y(i+1,j)-y(i+1,k))/(((x(i+1,j)-x(i+1,k))**2+(y(i+1,j)-y(i+1,k))**2+0.001)**(3/2))

else if(k==i)then
dx=dx
dy=dy
end if
!+++++++++++++++++++++++++++++++++++++++++++++++

end do


end do



end do

!!!!!!!+++++++++++++++++++++++++++


open(unit=3,file="galaxy2.txt",status="replace")


do i=1,2000,1

write(3,*) x(500,i), y(500,i)

end do



close(3)

open(unit=4,file="galaxy3.txt",status="replace")


do i=1,2000,1

write(4,*) x(50,i), y(50,i)

end do



close(4)



end program galaxy