Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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