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

comet_dust_time_r.py

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)