movement_equation.py

#!/usr/bin/env python

#  movement_equation_2.0.py
#
#  Copyright 2016 Nagy Gergely
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  :Author: Nagy Gergely
#  :Adress: nagygeri11@gmail.com
#  :Version: 0.2.1 beta
#  :Status: Prototype
#  :Date: 2016.01.05

import math
from PIL import Image, ImageDraw

#dimensions for vectors
x = 0
y = 1


"Non-physics functions"

def plot_base(xmin, xmax, ymin, ymax):
    """
    Creates a base for the plot: a white picture with two orthogonal lines, the x
    and y axis, according to the given minimum and maximum coordinates in pixels.
    """
    base_color = 'white'
    axis_color = 'grey'

    xsize = abs(xmax - xmin)
    ysize = abs(ymax - ymin)
    plot = Image.new('RGB', (xsize, ysize), base_color)
    draw = ImageDraw.Draw(plot)

    #draw x axis if shown:
    if ymin < 0 and ymax > 0:
        draw.line((0, ymax, xsize, ymax), axis_color)
    elif ymax > 0:
        draw.line((0, ysize, xsize, ysize), axis_color)
    else:
        draw.line((0, 0, xsize, 0), axis_color)

    #draw y axis if shown:
    if xmin < 0 and xmax > 0:
        draw.line((0-xmin, 0, 0-xmin, ysize), axis_color)
    elif xmax > 0:
        draw.line((0, 0, 0, ysize),axis_color)
    else:
        draw.line((xsize, 0, xsize, ysize),axis_color)

    return plot


def draw_path(m, r0, v0, force_function, time, dt, plot_size, resolution):
    """
    Draws the path of an object with weight m starting from r0 with velocity v0,
    according to the force function given in the force_function method.

    :param m: the mass of the object
    :param r0: the coordinates of the object at t=0 in meters (2-tuple)
    :param v0: the coordinates of the initial velocity vector (2-tuple)
    :param force_function: the function of the force applied to the object, described below
    :param time: the time of the movement
    :param dt: time elapsed ed between two calculated point: the bigger it is, the faster but less accurate is the result
    :param plot_size: the minimum and the maximum x and y coordinates as a 4-tuple: (xmin, xmax, ymin, ymax)
    :param resolution: resolution of the result plot image (meter/pixel)
    """
    trace_color = (255, 0, 0)
    plot_pixsize = [int(i/resolution) for i in plot_size]
    plot = plot_base(*plot_pixsize)
    draw = ImageDraw.Draw(plot)
    m = float(m)
    r = (float(r0[x]), float(r0[y]))
    v = (float(v0[x]), float(v0[y]))

    i = 0
    pos, prev_pos = None, None
    while float(i)*dt < time:

        #physics calculations
        F = force_function(m=m, r=r, v=v)
        a = (F[x]/m                        , F[y]/m                       )
        v = (v[x] + a[x]*dt                , v[y] + a[y]*dt               )
        r = (r[x] + v[x]*dt + a[x]/2*dt**2 , r[y] + v[y]*dt + a[y]/2*dt**2)

        pix_pos_x = round(r[x]/resolution) - plot_pixsize[0]
        pix_pos_y = round(r[y]/resolution) *-1 + plot_pixsize[3]
        if ((0 < pix_pos_x < plot_pixsize[1]-plot_pixsize[0]) and
            (0 < pix_pos_y < plot_pixsize[3]-plot_pixsize[2])):
            pos = pix_pos_x, pix_pos_y
        else:
            pos = None

        if pos is not None and prev_pos is not None:
            draw.line((prev_pos, pos), fill=trace_color)
            #plot.putpixel((pix_pos_x, pix_pos_y), trace_color)     #'dotty' but may visualize speed

        prev_pos=pos
        i += 1

    return plot



"""
*******************************************************************************
FORCE FUNCTIONS
*******************************************************************************
"""



"physics constants"

_c_  = 299792458 # m/s                  speed of light
_y_  = 6.67384 * (10**-11) # Nm²/kg²    gravitational constant
_h_  = 6.62606957 * (10**-34) # Js      Planck's constant
_E0_ = 8.854187817 * (10**-12) # C²/Nm² electric constant
_u0_ = 4.0 * math.pi # Tm/A             magnetic constant
_g_  = 9.80665 # m/s²                   standard gravity
_e_  = 1.602176565 * (10**-19) # C      elementary charge
_me_ = 9.10938291 * 10**-31 # kg        mass of electron
_mp_ = 1.672621777 * 10**-27 #kg        mass of proton
_Na_ = 6.02214129 * 10**23 #1/mol       Avogadro's constant

_Me_  = 5.972 * 10**24 # kg             mass of the Earth
_Ms_  = 1.989 * 10**30 # kg             mass of the Sun
_Mm_  = 7.34767309 * 10**22 #kg         mass of the Moon


"""
Force functions for the movement plotter. They should take the paramaters as keyword
arguments (use **kwargs to be compatible with more parameters in the future), and
return the coordinates of the force vector at these parameters as a 2-tuple.
Possible variablesat the moment: m, r, v
"""

def r_xy_dep(**variables):
    """force depends on the x and the y coordinates"""
    r = variables['r']
    F = [0.0 , -1.0]

    F[x] = 0
    F[y] = _g_*F[y]

    return F


def v_xy_dep(**variables):
    """force depends on the x and the y velocity"""
    v = variables('v')
    F = [0.0 , 1.0]
    F[x] = -(v[x]**2) * v[x]/abs(v[x])
    F[y] = -(v[y]**2) * v[y]/abs(v[y])

    return F


def central(**variables):
    """central force field, force depends on the
       vector from the centrum to the object"""

    r0 = variables['r']
    c = (0.0 , 0.0)               # centrum
    r = (r0[x]-c[x] , r0[y]-c[y]) # vector from centrum to r0
    r_ = math.hypot(r[x], r[y])   # length of r
    fi = math.atan2(r[y], r[x])   # angle of r

    F_ = -r_                      #force dependency
    dfi = 0                       #the angle between the force vector and r

    F_x = F_*math.cos(fi+dfi)
    F_y = F_*math.sin(fi+dfi)
    return (F_x, F_y)

def gravitational(**variables):
    #(special type of central dependency)

    central_mass = _Me_   # mass of the scource object
    m = variables['m']
    r0 = variables['r']
    c = (0.0 , 0.0)
    r = (r0[x]-c[x] , r0[y]-c[y])
    r_ = math.hypot(r[x], r[y])
    fi = math.atan2(r[y], r[x])
    F_ = -r_
    F_x = F_*math.cos(fi+dfi)
    F_y = F_*math.sin(fi+dfi)
    return (F_x, F_y)

def v_dep(**variables):
    """force depends on the velocity vector
       (i.e. charged particle in magnetic field)"""
    v_xy = variables['v']
    v_ = math.hypot(v_xy[x], v_xy[y])
    fi = math.atan2(v_xy[y], v_xy[x])

    F_ = v_
    dfi = math.pi/2

    F_x = F_*math.cos(fi+dfi)
    F_y = F_*math.sin(fi+dfi)
    return (F_x, F_y)



"""
******************************************************************************
MAIN MOVEMENT DRAWER FUNCTION
******************************************************************************
"""

def main(arg):

    m = 0.1               #mass of the object in kilograms
    r0 = (0, 0)           #initial coordinates in meters
    v0 = (30, 0)          #initial velocity vector in m/s
    dependency = v_dep    #force function from above
    time = 10             #time of movement to draw in secs
    time_res = 1/70000    #time between steps in secs, determines accuracy and running time
    plot_size = (-10, 10, -10, 10)  #minimum and maximum coordinates in meters (xmin, xmax, ymin, ymax)
    plot_res = 0.01       #meters per pixel on the plot

    plot = draw_path(m, r0, v0, dependency, time, time_res, plot_size, plot_res)
    plot.show()



if __name__ == '__main__':
    import sys
    sys.exit(main(sys.argv))