Skip to content

@rhettallain /comet_dust_diff_size.py secret
Created

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
This vpython model shows the motion of different sized dust particles released from a comet.
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
R1=.5e-6
R2=1e-6
R3=2.5e-6
R4=5e-6
rho=3000 #density of dust
dust1.m=rho*(4./3)*pi*R1**3 #mass = density*volume of sphere
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
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
tl=te/4.
dr0=vc*dt0 #this is the distance in one time interval (needed later)
dt=dr0/vc #calculate the time step based on the speed
while mag(ison.pos)<5*h:
rate(30)
#calculate the force
F=-G*Ms*ison.m*norm(ison.pos)/mag(ison.pos)**2
Fd1=-G*Ms*dust1.m*norm(dust1.pos)/mag(dust1.pos)**2+c*pi*R1**2*K*norm(dust1.pos)/mag(dust1.pos)**2
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+c*pi*R3**2*K*norm(dust3.pos)/mag(dust3.pos)**2
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
#update lines (for drawing)
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
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.