Created

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist

Vpython Comet Dust model. In this model, dust particles are released at different times.

View comet_dust_time_r.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 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
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.