Okay I from a past class in Astronomy. I grabbed this program and tweaked it a bit. I managed to add initial momentum into a given system, and I were attempting to simulate something of a galactic rotational plane (without expansion). Several variables worth noting for changing the simulator. L changes the distance between stars in general, while expanding the field of view. Its worth noting that you'd likely in order to see the stars better when say expanding out by or in by powers of ten correspondingly changing the
variable pertaining the visible radius set for stars.
This is found on the line
While I've tried using something of a centripetal acceleration for tangential motion of stars related to a given center of mass, the effort is off or even in this case, setting a larger scale mass (with zero initial momentum),
relative to a given star mass, I've still encountered some issues, namely, an expansion in the overall system. I've managed to tweak, however, initial momentum of stars in general so that rotation is about as stable as I've seen. My system is a bit slow in processing, so I couldn't remotely get close to something like our own galaxy in the number of stars, and likely one would be dealing with much greater length scales here, relative to the given length scales presented in L.
Adjusting the parameter cmass/6e10: vmags.append([vmag(cmass/6e10,rmag1)])
variable pertaining the visible radius set for stars.
This is found on the line
Stars = Stars+[sphere(pos=(x,y,z), radius=10*r, color=col0)]
While I've tried using something of a centripetal acceleration for tangential motion of stars related to a given center of mass, the effort is off or even in this case, setting a larger scale mass (with zero initial momentum),
relative to a given star mass, I've still encountered some issues, namely, an expansion in the overall system. I've managed to tweak, however, initial momentum of stars in general so that rotation is about as stable as I've seen. My system is a bit slow in processing, so I couldn't remotely get close to something like our own galaxy in the number of stars, and likely one would be dealing with much greater length scales here, relative to the given length scales presented in L.
Adjusting the parameter cmass/6e10: vmags.append([vmag(cmass/6e10,rmag1)])
from visual import *
from time import clock
from random import random,uniform
import numpy
import math
# Sets up a bunch of masses and lets them move under the influence of their
# mutual gravity.
# Masses are put in random positions within a sphere originally. An additional
# force is applied which simulates an equal density distributed over the
# rest of the universe outside this initial sphere.
def randomdirection():
# Generates random direction on sky.
# RA (in radians)
# Dec (in radians)
ra = 2.0*math.pi*random()
dec = math.acos(2.0*random()-1.0)-0.5*math.pi
return ra,dec
def ranvec():
# Generates a randomly orientated unit vector.
theta, phi = randomdirection()
z = math.sin(phi)
x = math.cos(phi)*math.sin(theta)
y = math.cos(phi)*math.cos(theta)
vec = numpy.array([x,y,z])
return vec
# Stars interacting gravitationally
# Program uses numpy arrays for high speed computations
Nstars = 100 # change this to have more or fewer stars
G = 6.7e-11 # Universal gravitational constant
# Typical values
Msun = 2E30
Rsun = 6.96e8
Rtrail = 2e8
L = 4e11
vsun = 0.9*sqrt(G*Msun/Rsun)
h0 = 2.5e-6 # Hubble's constant - expansion rate 2.5e-6 is good.
scene = display(title="Stars", width=1300, height=740,
range=L, forward=(-1,-1,-1))
Stars = []
poslist = []
plist = []
mlist = []
rlist = []
p0 = 0.0*Msun*100000.0
for i in range(Nstars):
vec = L*ranvec()*(random()**0.3333)
x = vec[0]
y = vec[1]
z = vec[2]/1e4
r = Rsun
col0 = (uniform(0.7,1.0),uniform(0.7,1.0),
uniform(0.7,1.0))
Stars = Stars+[sphere(pos=(x,y,z), radius=10*r, color=col0)]
mass = Msun
px = p0*uniform(-1,1)
py = p0*uniform(-1,1)
pz = p0*uniform(-1,1)
poslist.append((x,y,z))
plist.append((px,py,pz))
mlist.append(mass)
rlist.append(r)
pos = array(poslist)
p = array(plist)
m = array(mlist)
m.shape = (Nstars,1) # Numeric Python: (1 by Nstars) vs. (Nstars by 1)
radius = array(rlist)
vcm = sum(p)/sum(m) # velocity of center of mass
def centerofmass(m, pos):
# find center of mass
## print(m.shape)
## print(pos.shape)
## print(pos*m)
##
## print(sum(pos*m,axis=0))
R = sum(pos*m,axis=0)/sum(m)
return R
def rfromRvec(pos,R):
return pos-R
def unitnorm(r2):
v = []
for x in r2:
v.append([numpy.linalg.norm(x)])
v = array(v)
return r2/v
def tangentfromr2(r2):
def orthogonalangle(x,y):
return math.pi/2 + math.atan(y/x)
def findtancoords(rxpyp,oalpha):
return [rxpyp*math.cos(oalpha),
rxpyp*math.sin(oalpha)]
def findphi(rxy,z):
return math.atan(z/rxy)
def findrxy(x,y):
return (x**2+y**2)**.5
def findrxyz(x,y,z):
return (x**2+y**2+z**2)**.5
def rotationphi(phi):
m = [(math.cos(phi),math.sin(phi)),
(-math.sin(phi),math.cos(phi))]
return numpy.array(m)
def revrotationphi(phi):
m = [(math.cos(phi),-math.sin(phi)),
(math.sin(phi),math.cos(phi))]
return numpy.array(m)
def rotatecoords(x,z,R):
xzcrd = [x,z]
XZ = numpy.array(xzcrd)
rXZ = numpy.dot(R,XZ)
return rXZ
# to find tangent one solution, involves using rotated coordinate system
# on xz plane by phi where phi is represented by the angle of r2.
# Here we then have under rotated coordinates a tangent vector given
# by a perpendicular on the 2 dimensional plane by pi/2, then covert
# computed unit tangent vector back to our original coordinate system.
# a little easier despite the fuss of coordinate transformations since
# computing a tangent on a 2 dimensional plane is a little less cumbersome.
#first find r for xy plane
x=r2[0];y=r2[1];z=r2[2]
rxy = findrxy(x,y)
theta = findphi(x,y)
if theta < -math.pi:
print(theta)
phi = findphi(rxy,z)
if phi <0:
print(phi)
Rxy = rotationphi(theta)
Rxz = rotationphi(phi)
revRxz = revrotationphi(phi)
revRxy = revrotationphi(theta)
#phi positive is 1rst or 3rd quad
#theta positive is 1rst or 3rd quad
XpYp = rotatecoords(x,y,Rxy)
#z = zp
xp = XpYp[0]; yp = XpYp[1]; zp = z;
XppZpp = rotatecoords(xp,zp,Rxz)
xpp = XppZpp[0]; zpp = XppZpp[1]; ypp = yp
# under coordinate rotations ypp = yp,
## orthetap = 1*math.pi/2
#here a twice rotated coordinate system puts the zpp at zero, and likewise ypp should be zero)
rxppypp = findrxy(xpp,ypp)
typp, txpp = (xpp,ypp)
## tangentcoords = findtancoords(rxppypp, orthetap)
## txpp,typp = tangentcoords
if typp < 0:
print(typp)
## typp = abs(typp); txpp = abs(txpp)
#convert txp, zp back to original coordinates, typ = ty
Txpzp = rotatecoords(txpp,zpp,revRxz)
txp = Txpzp[0]; typ = typp; tzp = Txpzp[1]
Txy = rotatecoords(txp,typ,revRxy)
tx = Txy[0]; ty = Txy[1]; tz=tzp
return [tx,ty,tz]
def vmag(M,rmag):
# M is the total mass, and rmag is the distance from center of total mass M to star
G = 6.7e-11
return (G*M/rmag)**.5
def angles(pos):
#find spherical coordinate angles
c = pos[:,1]/pos[:,0]
d = (pos[:,0]**2 + pos[:,1]**2)**.5
e = pos[:,2]/d
f = []; g = [];
i = 0
for x in c:
f.append([math.atan(x),math.atan(e[i])])
i += 1
return numpy.array(f)
def finddcoords(rxy,oalpha):
r2 = []
for oalphai in oalpha:
theta,phi = oalphai
rx = rxy*math.sin(theta)*math.cos(phi)
ry = rxy*math.sin(theta)*math.sin(phi)
rz = rxy*math.cos(theta)
r2.append([rx,ry,rz])
return numpy.array(r2)
def getCMlist(r2,m):
magnr2 = (r2[:,0]**2+r2[:,1]**2+r2[:,2]**2)**.5
i = 0;Mlist =[];
for m2 in magnr2:
n=0;M=0
for m2a in magnr2:
if not i == n:
if m2a < m2 or m2a == m2:
M += m[n]
n+=1
Mlist.append(M)
i+=1
return numpy.array(Mlist)
m[0] = sum(m)/1
cmass = sum(m)
R = centerofmass(m, pos)
print('center of mass:')
print(R)
pos[0] = R
R = centerofmass(m, pos)
M = sum(m)
#Rsun in the direction of r
ang = angles(pos)
rSunarr = finddcoords(Rsun,ang)
r2 = rfromRvec(pos,R) + rSunarr
Mlist = getCMlist(r2,m)
print(r2)
r2n = unitnorm(r2)
print(r2n)
tpos = []
vmags = []
itr = 0
print('mlist:')
print(Mlist)
for r2i in r2n:
txyz = tangentfromr2(r2i)
tpos.append(txyz)
x = r2[itr,0]; y=r2[itr,1]; z=r2[itr,2];
rmag1 = (x**2+y**2+z**2)**.5
## vmags.append([vmag(cmass/6e10,rmag1)])
vmags.append([vmag(Mlist[itr],rmag1)])
itr+=1
tposa = array(tpos)
print(tposa)
vmagsa = array(vmags)
print(vmags)
vs = tposa*vmagsa
print(vs)
p = p-m*vcm # make total initial momentum equal zero
p = vs*Msun
p[0] = [0,0,0]
print(p)
dt = 100.0
Nsteps = 0
pos = pos+(p/m)*(dt/2.) # initial half-step
time = clock()
Nhits = 0
while 1:
rate(10000)
L *= 1.0+h0*dt
con = 0.9*G*Nstars*Msun/(L*L*L)# strength of force to allow for external mass
# Compute all forces on all stars
try: # numpy
r = pos-pos[:,newaxis] # all pairs of star-to-star vectors
for n in range(Nstars):
r[n,n] = 1e6 # otherwise the self-forces are infinite
rmag = sqrt(add.reduce(r*r,-1)) # star-to-star scalar distances
hit = less_equal(rmag,radius+radius[:,newaxis])-identity(Nstars)
hitlist = sort(nonzero(hit.flat)[0]).tolist() # 1,2 encoded as 1*Nstars+2
F = G*m*m[:,newaxis]*r/rmag[:,:,newaxis]**3 # all force pairs
except: # old Numeric
r = pos-pos[:,newaxis] # all pairs of star-to-star vectors
for n in range(Nstars):
r[n,n] = 1e6 # otherwise the self-forces are infinite
rmag = sqrt(add.reduce(r*r,-1)) # star-to-star scalar distances
hit = less_equal(rmag,radius+radius[:,newaxis])-identity(Nstars)
hitlist = sort(nonzero(hit.flat)) # 1,2 encoded as 1*Nstars+2
F = G*m*m[:,newaxis]*r/rmag[:,:,newaxis]**3 # all force pairs
for n in range(Nstars):
F[n,n] = 0 # no self-forces
p = p+sum(F,1)*dt+pos*con*dt*m
# Having updated all momenta, now update all positions
pos = pos+(p/m)*dt
# Expand universe
#pos *= 1.0+h0*dt
# Update positions of display objects; add trail
for i in range(Nstars):
Stars[i].pos = pos[i]
# If any collisions took place, merge those stars
for ij in hitlist:
i, j = divmod(ij,Nstars) # decode star pair
if not Stars[i].visible: continue
if not Stars[j].visible: continue
# m[i] is a one-element list, e.g. [6e30]
# m[i,0] is an ordinary number, e.g. 6e30
newpos = (pos[i]*m[i,0]+pos[j]*m[j,0])/(m[i,0]+m[j,0])
newmass = m[i,0]+m[j,0]
newp = p[i]+p[j]
newradius = Rsun*((newmass/Msun)**(1./3.))
iset, jset = i, j
if radius[j] > radius[i]:
iset, jset = j, i
Stars[iset].radius = newradius
m[iset,0] = newmass
pos[iset] = newpos
p[iset] = newp
Stars[jset].visible = 0
p[jset] = vector(0,0,0)
m[jset,0] = Msun*1E-30 # give it a tiny mass
Nhits = Nhits+1
pos[jset] = (10.*L*Nhits, 0, 0) # put it far away
Nsteps += 1
No comments:
Post a Comment