fishies.pyde

# particle physics simulation of charged particles in a toroidal field with electrostatic repoulsion
# (basically all the fish are attracted to an invisible rotating donut, and repelled by each other
# assumes all particles of equal mass and charge

# constant simulation parameters:

# screen size
SS = 800
# number of particles
NP = 16

# sphere radius
# the radius of the main sphere of each particle
SR = 20

# toroidal radius
# radius of the donut that the fish are attracted to
TR = 200
# toroidal force
TG = 100

# inverse square threshold
# if the distance between particles is less than ist, their inverse square force
# is bounded to prevent simulation blowing up.
IST = 0.00001

class DVec(object):
    """Dev's PVector methods"""

    @classmethod
    def random3D(cls):
        return PVector(random(-1, 1), random(-1, 1), random(-1, 1)).normalize()

    @classmethod
    def rotateX(cls, vec, t):
        """return this vec rotated around the X axis by t radians"""
        return PVector(
            vec.x,
            vec.y * cos(t) - vec.z * sin(t),
            vec.y * sin(t) + vec.z * cos(t)
        )

    @classmethod
    def inv_sq(cls, vec, other):
        """
        inverse square force experienced towards `other`, bounded using `ist`
        """
        global IST
        d = other - vec
        ds = sq(d.mag())
        # if the square of the distance is too small, don't divide by it
        if(ds > IST):
            d.mult(1/ds)
        return d
    
    @classmethod
    def xyHeading(cls, vec):
        """The angle of the vector projected onto the xy plane"""
        return atan2(vec.y, vec.x)
    
    @classmethod
    def zHeading(cls, vec):
        """The angles between the vector and the Z axis"""
        return PVector.angleBetween(vec, PVector(0, 0, 1))
        
        
class Particle(object):
    """A particle with position, velocity, mass and charge"""

    def __init__(self, p, v, m=1, q=1):
        super(Particle, self).__init__()
        self.p = p
        self.v = v
        self.m = m
        self.q = q
        self.h = []

    def step(self, a):
        """compute next step in simulation, given particle's acceleration"""
        # increment position by velocity
        self.p = self.p + self.v
        # increment velocity by accelleration
        self.v = self.v + a
        # slow down the velocity with a drag factor
        self.v.mult(0.99)

    def draw(self):
        # draw a sphere at the particle's current position
        pushMatrix()
        noStroke()
        fill(255, 0, 0)
        translate(self.p.x, self.p.y, self.p.z)
        tail = -1*self.v
        # orient the fish so that its tail is pointing towards x
        rotateY(DVec.zHeading(tail)-QUARTER_PI)
        rotateZ(DVec.xyHeading(tail))
        pushMatrix()
        scale(1.2, 0.4, 0.8)
        sphere(SR)
        popMatrix()
        # draw the tail fin
        translate(SR, 0, 0)
        scale(1, 0.4, 1)
        sphere(SR/1.5)
        popMatrix()

    def calculate_toroidal(self):
        """
        The toroidal acceleration experienced by any particle at point p

        The torus is a donut on the x-y plane that all particles are attracted to
        """
        global ts  # toroidal spin
        
        # phi is the angle the particle makes with the torus
        ph = DVec.xyHeading(self.p)

        # tp is the closest point to the torus
        tp = DVec.rotateX(PVector(TR*sin(ph),  TR*cos(ph), 0), ts)
        # tn is on the opposite side of the torus
        tn = DVec.rotateX(PVector(-TR*sin(ph),  -TR*cos(ph), 0), ts)

        tf = TG * (DVec.inv_sq(self.p, tp)
                    + DVec.inv_sq(self.p, tn)
                    )

        return tf

    def calculate_electrostatic(self):
        """The electrostatic acceleration for a particle"""
        global ps  # list of particles
        global eg  # electrostatic force
        assert(isinstance(self, Particle))
        tf = PVector(0, 0, 0)
        for q in ps:
            # print "tf: ", tf
            tf = tf + DVec.inv_sq(self.p, q.p)
        tf.mult(eg)
        return tf

    def __repr__(self):
        return "Particle({}, {})".format(
            repr(self.p),
            repr(self.v)
        )


def setup():
    print("Starting setup")
    # energy previously
    # how much energy the system had last iteration
    global ep
    ep = 0
    # electrostatic force
    # how much particles are attracted to each other: the toroidal force
    # divided by the number of particles
    global eg
    eg = - TG / NP
    # toroidal spin
    # how much the torus is rotated in the x axis, is adjusted to give the particles
    # an extra kick of energy when they start to slow down
    global ts
    ts = 0
    # list of particles
    global ps
    ps = []

    for _ in xrange(NP):
        p = Particle(DVec.random3D() * TR/3, DVec.random3D() * TR/50)
        print(p)
        ps.append(p)
    size(SS, SS, P3D)
    
    sphereDetail(10)


def draw():
    # print "Starting draw"
    global ps  # list of particles
    global ep  # energy previously
    global ts  # toroidal spin
    background(color(54, 70, 93, 10))
    
    camera(0, 0, SS, 0, 0, 0, 1, 0, 0)
    
    # lighting stuff
    directionalLight(102, 102, 102, 0, 0, -1)
    lightSpecular(204, 204, 204)
    directionalLight(102, 102, 102, 0, 1, -1)
    lightSpecular(102, 102, 102)
    specular(51, 51, 51)

    noStroke()
    fill(255, 127, 0)

    for p in ps:
        # print(p)
        p.draw()

    # move particles to next frame, calculate the total energy
    ec = 0
    for p in ps:
        p.step(p.calculate_toroidal() + p.calculate_electrostatic())
        ec += p.v.mag()
    # if we're losing energy, spin the torus!
    if(ep != 0):
        ed = (ec - ep) / NP
        if(ed < 0):  # system is losing energy
            ts += 0.005
            print("boosting spin to {}".format(ts))
    ep = ec