Skip to content

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Vpython Comet Dust model. In this model, dust particles are released at different times.
from visual import *
K=1.975e17 #this is the radiation pressure constant for the sun
G=6.67e-11
Ms=1.9891e30 #mass of sun
Rs=6.955e8 #radius of sun
m=1000 #mass of the comet (it doesn't matter)
s=1#s is the scale so I can see stuff
h = 1.8e9 #closest approach
vc=4e5 #just a guess here
c=1.5 #reflectivity of dust
sun=sphere(pos=(0,0,0), radius = s*Rs, color=color.yellow)
ison=sphere(pos=(-h,0,0), radius = .1*s*Rs, make_trail=True)
dust1=sphere(pos=ison.pos, radius=ison.radius*2, color=color.cyan)
dust2=sphere(pos=ison.pos, radius=ison.radius*2, color=color.cyan)
dust3=sphere(pos=ison.pos, radius=ison.radius*2, color=color.cyan)
dust4=sphere(pos=ison.pos, radius=ison.radius*2, color=color.cyan)
#The lines are just connecting the dust particles together for visual effect
line1=cylinder(pos=ison.pos, axis=(dust4.pos-ison.pos), radius=dust4.radius, opacity=.5)
line2=cylinder(pos=dust4.pos, axis=(dust3.pos-dust4.pos), radius=dust4.radius, opacity=.5)
line3=cylinder(pos=dust3.pos, axis=(dust2.pos-dust3.pos), radius=dust4.radius, opacity=.5)
line4=cylinder(pos=dust2.pos, axis=(dust1.pos-dust2.pos), radius=dust4.radius, opacity=.5)
ion=cylinder(pos=ison.pos, axis=mag(dust2.pos-ison.pos)*norm(ison.pos), radius=dust4.radius, color=color.cyan, opacity=.5)
#these are the different radius of dust
#R=.5e-6
R1=.5e-6
R2=1e-6
R3=2.5e-6
R4=5e-6
R2=R1
R3=R1
R4=R1
rho=3000 #density of dust
#mass = Volume of sphere * density
dust1.m=rho*(4./3)*pi*R1**3
dust2.m=rho*(4./3)*pi*R2**3
dust3.m=rho*(4./3)*pi*R3**3
dust4.m=rho*(4./3)*pi*R4**3
ison.m=m
#initial momentum
ison.p=vector(0,-vc,0)*ison.m
dust1.p=vector(0,-vc,0)*dust1.m
dust2.p=vector(0,-vc,0)*dust2.m
dust3.p=vector(0,-vc,0)*dust3.m
dust4.p=vector(0,-vc,0)*dust4.m
dt0=100
t=0
te=90300 #this is the time at the end for a previous run
tl=te/4. #this is the interval that dust is released
dr0=vc*dt0 #this is a constant to create a variable time interval
dt=dr0/vc #this is the calculation time interval
tmax=38000 #end of simulation time
rt = 38000/6 #the time between releases of dust
while mag(ison.pos)<5*h:
rate(30)
#gravitaitonal force on the comet
F=-G*Ms*ison.m*norm(ison.pos)/mag(ison.pos)**2
#the first dust has both gravity and light pressure
Fd1=-G*Ms*dust1.m*norm(dust1.pos)/mag(dust1.pos)**2+c*pi*R1**2*K*norm(dust1.pos)/mag(dust1.pos)**2
#For the other dusts, they just have gravity unless it is after a time interval
Fd2=-G*Ms*dust2.m*norm(dust2.pos)/mag(dust2.pos)**2
if t>rt:
Fd2=-G*Ms*dust2.m*norm(dust2.pos)/mag(dust2.pos)**2+c*pi*R2**2*K*norm(dust2.pos)/mag(dust2.pos)**2
Fd3=-G*Ms*dust3.m*norm(dust3.pos)/mag(dust3.pos)**2
if t>2*rt:
Fd3=-G*Ms*dust3.m*norm(dust3.pos)/mag(dust3.pos)**2+c*pi*R3**2*K*norm(dust3.pos)/mag(dust3.pos)**2
Fd4=-G*Ms*dust4.m*norm(dust4.pos)/mag(dust4.pos)**2
if t>3*rt:
Fd4=-G*Ms*dust4.m*norm(dust4.pos)/mag(dust4.pos)**2+c*pi*R4**2*K*norm(dust4.pos)/mag(dust4.pos)**2
#update momentum
ison.p=ison.p+F*dt
dust1.p=dust1.p+Fd1*dt
dust2.p=dust2.p+Fd2*dt
dust3.p=dust3.p+Fd3*dt
dust4.p=dust4.p+Fd4*dt
#update position
ison.pos=ison.pos+ison.p*dt/ison.m
dust1.pos=dust1.pos+dust1.p*dt/dust1.m
dust2.pos=dust2.pos+dust2.p*dt/dust2.m
dust3.pos=dust3.pos+dust3.p*dt/dust3.m
dust4.pos=dust4.pos+dust4.p*dt/dust4.m
#draw the lines update position and axis
line1.pos=ison.pos
line1.axis=dust4.pos-ison.pos
line2.pos=dust4.pos
line2.axis=dust3.pos-dust4.pos
line3.pos=dust3.pos
line3.axis=dust2.pos-dust3.pos
line4.pos=dust2.pos
line4.axis=dust1.pos-dust2.pos
ion.pos=ison.pos
ion.axis=mag(dust2.pos-ison.pos)*norm(ison.pos)
t=t+dt
#here I recalculate the time interval size based on the speed of the comet
dt=dr0*ison.m/mag(ison.p)
print(ison.pos)
print(ison.p)
#print(mag(ison.pos-dust.pos)/1e9)
print(t)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.