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

comet_dust_diff_size.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
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)