Created

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist

This vpython model shows the motion of different sized dust particles released from a comet.

View comet_dust_diff_size.py
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
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.