debug prints tweak

[physics.git] / particles.py

 #!/bin/usr/python
# Copyright 2009 Daniel Cardenas All rights reserved.
# Sub atomic particles as standing electromagnetic waves.
#   Strong and weak nuclear force are manifestations of electromagnetic forces.
# Algorithm:
#    Set up:
#        Initial time quantum
#        Number of particles
#            For each particle:
#                Position
#                Momentum
#                Force on other particles? (nah, just calulate it twice)
#        Trail
#
# Energy problem with this simulation:
#   When parts are very close to each other the electric force is very strong.  Approaching infinity as distance approaches zero.
#   Simple simulation with discrete time quantum can't handle this problem.
#   When electron is approaching proton its speed is increasing.  When passing proton it has higher speed.
#   Because of higher speed it is further away and therefore won't lose all of the energy gain.
#   Lets call this digital energy gain.
#
#    Simulation has error, energy gain, because wl.dt is not infinitely small.
#    When electron approaches proton there are more samples taken as electron approaches proton then
#    when leaving proton, because aproaching proton speed is slower and when leaving speed is higher.
#    Gaining or Losing Energy: http://www.schlitt.net/xstar/n-body.pdf
#    Simulated orbits spiral outwards: http://burtleburtle.net/bob/math/multistep.html
#    http://www.amara.com/papers/nbody.html
#    http://www.physicsforums.com/showpost.php?p=677453&postcount=5
#
# What are options for compensating for digital energy gain?
#  1. Only lose energy when in escape mode
#  2. Variable wl.dt, variable time quantum.  Change the time quantum so it loses energy leaving close particle.
#  3. Lose energy after gained a lot?  Track energy of the system.  Problem how to calculate energy of system with infinite energy.
#  4. Figure out the path and velocity that the electron would really take and have.  Difficult since electric fields are changing.
#      Possible if use 3rd order form where acceleration is 2nd order when constant.  3rd order form involves jerk.
#      Difficult when several protons and one electron involved.
#  5. Learn about skip method
#  6. Search for and study other methods on the internet
#  7. Collect data when it occurs and thus make smarter decision
#  8. Use special relativity.  Didn't help.
#  9. Use integration algorithm for calculating velocity rather than dt method.
#      Requires equation for velocity over the path
#  10. If particle is outside of the normal size of an atom then assume it is escaping due to energy gain.
#       reset kinetic energy to zero.
#
# Color scheme:
#   Need colors for:
#       1.  x axis           torquoise, greenish, blueish
#       2.  y axis           yellow
#       3.  electron         redish,  a bit of random color added to distinguish individual electrons
#       4.  proton           blueish, a bit of random color added to distinguish individual protons
#       5.  direction        green, equivalent to velocity 
#       6.  total    force   white
#       7.  electric force   red
#       8.  magnetic force   blue  ,  Direction that external magnetic field is pushing us
#       9.  magnetic moment  yellow (.5,.5,0)    Intrinsic magnetic field or spin as popular known.
#       10. velocity         green
#       11.
#
# Variable comments:
#   intrinsic_B_field_plus_qv - This is a vector using the right hand rule.
#                                       It points in the direction of v in qv for example.
#                                       The B field rotates around the direction that this vector points.
#
#
from visual import *    # See vpython.org for more info.
import math
import random           # http://effbot.org/pyfaq/how-do-i-generate-random-numbers-in-python.htm
#import threading               # Have computations on one thread and drawing on another
import pprint
import sys              # for argv, flush
import getopt

class simple_class:         # Used to hold constants for electron and proton and other global variables.
    """An empty class"""

wl = simple_class()            # workload.  Holds global variables related to our set of particles or system of particles

wl.particle_count = 2       # hydrogen
#wl.particle_count = 3       # Neutron and proton
#wl.particle_count = 4  # heavy hydrogen, Deuterium
#wl.particle_count = 8       # Helium, 2 neutrons
#wl.particle_count = 14      # lithum
wl.size_adjust    = wl.particle_count / 2   # Larger particles when there are more of them because they are further apart.
                                        # See setup_globals().  Just eye candy improvement.
debug_verbosity = 20
debug_sum_forces_on_one_particle = 21
random.seed(0)          # Can comment out for real randomness.
c     = 299792458       # speed of light, http://en.wikipedia.org/wiki/Speed_of_light
c2    = c**2            # Used for Lorentz factor       http://en.wikipedia.org/wiki/Lorentz_factor
h_eVs = 4.135667e-15    # http://en.wikipedia.org/wiki/Planck_constant
h_Js  = 6.62606896e-34  # http://en.wikipedia.org/wiki/Planck_constant
Ke    = 8.987551787e9   # Proportionality constant, called Coulomb's constant.  http://en.wikipedia.org/wiki/Electric_force
u0    = 1.2566e-6       # Henry's per meter or Tesla*meters/Ampere.             http://en.wikipedia.org/wiki/Magnetic_constant
coulomb = 6.241509629e18 # http://en.wikipedia.org/wiki/Coulomb
size_of_hydrogen_nucleus = 1.6e-15              # In meters.   http://en.wikipedia.org/wiki/Atomic_nucleus
B_field_in_space   =  1.0e-11         # Teslas
earth_B_field      = vector(0,0,40.0e-06)   # http://en.wikipedia.org/wiki/Earth%27s_magnetic_field#Field_characteristics.  Arbitrarily assigned in the B field direction
background_magnetic_field = earth_B_field
elementary_charge  =  1.602176487e-19 # coulombs
base_electric_pot_energy = Ke * ( elementary_charge ** 2 ) / size_of_hydrogen_nucleus
K_m     = u0 /(4*pi)     # magnetic constant.  http://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law


proton   = simple_class()                    # Just defines some constants
electron = simple_class()
proton  .radius   = .875e-15                 # http://en.wikipedia.org/wiki/Proton.  Somewhat arbitrary.  Doesn't effect experiment.  Just for visual candy.
proton  .visual_radius = proton.radius
electron.visual_radius = proton.visual_radius * .5     # Arbitrary, just for visual candy
proton  .rest_mass_kg  = 1.672621636e-27         # kg
proton  .mass_MeV = 938.272013
proton  .mass_eV  = proton.mass_MeV * 1e6
electron.rest_mass_kg  = 9.10938215e-31      # http://en.wikipedia.org/wiki/Electron
electron.mass_MeV = 0.51099891               # http://en.wikipedia.org/wiki/Electron and http://en.wikipedia.org/wiki/Electronvolt
electron.mass_eV  = electron.mass_MeV * 1e6
electron.frequency = ( 2 * electron.mass_eV ) / h_eVs   # e = h*f  , f = e /h, 2 because that is what is required from pair production.
proton  .frequency = ( 2 *   proton.mass_eV ) / h_eVs   # e = h*f  , f = e /h
electron.wavelength = h_Js / ( electron.rest_mass_kg * c )   # http://en.wikipedia.org/wiki/Compton_wavelength
electron.period   = 1 / electron.frequency
proton  .period   = 1 /   proton.frequency
print "\t electron frequency ", electron.frequency, "  period ", electron.period, electron.wavelength
print "\t   proton frequency ",   proton.frequency, "  period ",   proton.period
proton  .charge   = elementary_charge
electron.charge   = elementary_charge * -1
electron.magnetic_moment = -928.476377e-26     # http://en.wikipedia.org/wiki/Electron_magnetic_dipole_moment
proton  .magnetic_moment =   14.106067e-27     # http://en.wikipedia.org/wiki/Magnetic_moment#Elementary_particles
                                               # Doesn't make sense that a proton would have a small magnetic moment.
                                               # Joules / Teslas.  Is that the unit of measure we need it in?
electron.magnetic_moment = elementary_charge   # Experimenting
proton  .magnetic_moment = elementary_charge   # Experiment

el_pr_weight_diff = proton.rest_mass_kg / electron.rest_mass_kg   # Used for arbitrary initial set up.
#print " el_pr_weight_diff %d" % el_pr_weight_diff
electron.radius   = proton.radius / el_pr_weight_diff       # Just a wild guess.  Maybe the electron true radius is larger than that of the proton.


scene.title="Sub atomic particles as standing electromagnetic waves"
scene.width =800        # go ahead and change these to suit your preference.
scene.height=600
scene.x     =0
scene.y     =0

axis_length = size_of_hydrogen_nucleus * 4 * wl.size_adjust
y_axis_label = label(pos=(0,axis_length,0), text="Y", opacity=.3, box=0, line=0)
y_axis = arrow(pos=(0,0,0), axis=(0,size_of_hydrogen_nucleus,0), length=axis_length, opacity=0.7, color=(.5,.5,0))
y_axis.shaftwidth  = y_axis.shaftwidth  / 4
y_axis. headwidth  = y_axis. headwidth  / 4
y_axis. headlength = y_axis. headlength / 4

x_axis_label = label(pos=(axis_length,0,0), text="X", opacity=.3, box=0, line=0)
x_axis       = y_axis.__copy__()
x_axis.axis  = (size_of_hydrogen_nucleus,0,0)
x_axis.color = (0,.5,0.5)
x_axis.length= axis_length

initial_velocity_max = 3e4                 # arbitrary

max_ptr_len = size_of_hydrogen_nucleus * wl.size_adjust
min_ptr_len = max_ptr_len / 2

#scene.show_rendertime = True

kb_input_label = label(visible=0)

def setup_workload() :
    wl.number_of_dt_per_proton_period = 16
    #wl.number_of_dt_per_proton_period = 32
    #wl.number_of_dt_per_proton_period = 64
    #wl.number_of_dt_per_proton_period = 128
    wl.desired_dt = proton.period / wl.number_of_dt_per_proton_period  # wl.dt = time quantum or delta time.
    wl.dt = wl.desired_dt  # wl.dt = time quantum or delta time.  Adjusts to smaller values when particles are close.
    wl.old_dt = wl.dt
    print "\t wl.dt %1.1e" % wl.dt
    wl.last_dt_change = wl.dt

    #def update_dt():
    # If the electrical force is strong I want to have a smaller dt so that we can model orbits better
    # Also have smaller digital energy gain issue.
    #new_dt = large_dt / wl.largest_ele_force
    # I want to lower dt when electrical force is greater than 50.
    #new_dt = large_dt / 32
    # large_dt = desired_dt * 32
    wl.large_dt = wl.desired_dt * 32     # From running simulation I know that electron gets thrown when electrical force reaches this value.
    wl.pars = []               # wl.pars = sub-atomic particles, a list of vpython sphere objects with added attributes.
    wl.total_time = 0
    wl.main_loop_count = 0      # Just a counter for number of iterations we've done.
    wl.curr_KE_system    = 0.0    # KE = Kinetic Energy
    wl.orig_KE_system    = 0.0    # KE = Kinetic Energy
    wl.prev_KE_system    = 0.0
    wl.part_with_most_energy = 0;
    wl.play_mode = True             # if s == 's' or s == 'p' or s == ' ':   #stop, pause or play
    wl.auto_zoom = True

def print_unit_vector(string, vec) :     # Vector should already be in normalized form
    print "%s%2d" % (string, int(vec.x*10)),
    print "%2d" % int(vec.y*10),
    print "%2d" % int(vec.z*10),

def find_perpendicular(vel,to_proton_vec):
    # http://local.wasp.uwa.edu.au/~pbourke/geometry/disk/
    count = 0
    while True:
        # http://local.wasp.uwa.edu.au/~pbourke/geometry/disk/
        # Step 1.  Choose any point P randomly which doesn't lie on the line through P1 and P2
        # Step 2. Calculate the vector R as the cross product between the vectors
        #         P - P1 and P2 - P1. This vector R is now perpendicular to P2 - P1.
        #         (If R is 0 then step 1 wasn't satisfied)
        r = cross( vel, to_proton_vec )
        r_norm = norm(r) * -1
        #print_unit_vector("vel", norm(vel)); print_unit_vector(" r", r_norm),
        if r.x != 0 or r.y != 0 or r.z != 0:
            break
            print "\t Cross product is not zero.  Success."
        else:
            print "\t Our points lie on the same line.",
            print vel, to_proton_vec,
            if count > 3 :
                print "Can't find random points that are not on the same line."
                exit(1)
        random_vec = vector(random.random(),random.random(),random.random());
        if count%2 == 0 :
            vel = random_vec
        else :
            to_proton_vec = random_vec
        count += 1
        print "random vec ", random_vec
    return r_norm


def setup_particles():
    initial_spacing = size_of_hydrogen_nucleus * 16 * ( wl.particle_count / 2 ) # arbitrary
    electron_count = 0
    proton_count = 0
    num_particles_of_same_type = int ( (wl.particle_count / 2) + 0.5 )
    for i in range(0,wl.particle_count):
        # For describing sub atomic particles particles we start with visual python sphere object then add to it.
        wl.pars.append(sphere(radius=electron.visual_radius))   # Specify a radius that gets reset so that we don't get flicker with autoscale.
        part = wl.pars[i]           # part = particle.  Create a pointer to the sphere object which holds most of our data.
        part.index = i              # Set one of the many variables in our particle object.
        part.num   = i
        print "num %d " % i,
        #part.is_proton   = (i%2)==0;
        part.is_proton   = i < ( (wl.particle_count+1) / 2 );
        part.is_electron = not part.is_proton
        # For info on vector see: http://vpython.org/webdoc/visual/vector.html
        spacing = initial_spacing
        need_random_loc_assigned = True
        if part.is_electron and i >= 2 :                # Put electron between protons.
            need_random_loc_assigned = False
            #if part.is_electron and ((i+1)%4) == 0 :                # Make a neutron.  Put every other electron between previous protons.
            #part.pos = mid_point(previous_proton, prev_prv_proton)
            proton_index = int((i-1)/2)
            p  = wl.pars[proton_index]
            np = wl.pars[proton_index+1]        # np = next proton
            if np.is_electron :
                if i > 4 :
                    np = wl.pars[0]         # put electron between first proton and last proton.
                else :
                    need_random_loc_assigned = True
            if not need_random_loc_assigned :
                part.pos = vector(p.x-((p.x-np.x)/2),p.y-((p.y-np.y)/2),p.z-((p.z-np.z)/2))
                print "Putting electron in the middle of two protons"
        if need_random_loc_assigned :
            if part.is_proton :
                spacing = spacing / 2
            part.pos = vector((random.random()-.5)*spacing,    # Arbitrary values assigned
                              (random.random()-.5)*spacing,
                              (random.random()-.5)*spacing)
        part.opacity = 0.3                      # Arbitrary.  Just chose what I think looks good.
        #print i, 'is electron', part.is_electron
        part.velocity = - norm(part.pos)
        part.vel_norm = norm(part.velocity)
        part.vel_pointer = arrow(pos=part.pos, axis=(1,0,0), length=max_ptr_len, color=(0,1,0), opacity=.2)  # direction of velocity
        if part.is_electron :
            part.ptype_name  = "electron"
            part.charge_avg  = electron.charge
            part.frequency   = electron.frequency
            part.rest_mass_kg= electron.rest_mass_kg
            part.period      = electron.period
            part.radius      = electron.visual_radius * wl.size_adjust
            part.color       = (1.0,random.random()**2,random.random()**2)   # Redish color, ^2 the other colors makes them closer to zero
            part.EF_pointer = arrow(pos=part.pos, length=max_ptr_len, color=(1,0,0), opacity=.2, axis=(1,0,0))  # Where is the electric force pushing us?
            part.MF_pointer = arrow(pos=part.pos, length=max_ptr_len, color=(0,0,1), opacity=.2, axis=(1,0,0))  # Where is the magnetic force pushing us?
            part.tot_pointer = arrow(pos=part.pos, length=max_ptr_len, color=(1,1,1), opacity=.2, axis=(1,0,0))  # Where is the total force pushing us?
            part.mag_moment_unit_vec = vector(0,0,-1)        # initial z direction is arbitrary
            part.      avg_magnetic_moment = electron.magnetic_moment;
            part.intrinsic_magnetic_moment = part.mag_moment_unit_vec * part.avg_magnetic_moment;  # http://en.wikipedia.org/wiki/Spin_%28physics%29#Magnetic_moments
            part.this_type_count = electron_count
            electron_count += 1;
        else:
            part.ptype_name  = "proton  "
            part.charge_avg  = proton.charge * 1
            part.frequency   = proton.frequency
            part.rest_mass_kg= proton.rest_mass_kg
            part.period      = proton.period
            part.radius      = proton.visual_radius * wl.size_adjust
            part.color       = (random.random()**2,random.random()**2,1.0)   # Blue-ish color, ^2 the other colors makes them closer to zero
            part.mag_moment_unit_vec = vector(0,0,1)        # initial z direction is arbitrary
            part.      avg_magnetic_moment = proton.magnetic_moment;
            part.intrinsic_magnetic_moment = part.mag_moment_unit_vec * part.avg_magnetic_moment   # http://en.wikipedia.org/wiki/Spin_%28physics%29#Magnetic_moments
            part.this_type_count = proton_count
            proton_count += 1;
        part.previous_magnetic_moment = vector(0,0,-1)

        #part.type = part.ptype_name
            #print "part.radius", part.radius
            # magnetic Dipole moment pointer.  http://en.wikipedia.org/wiki/Magnetic_dipole_moment
        part. mm_pointer = arrow(pos=part.pos, length=part.radius*2, color=(0.5,0.5,0.0), opacity=.2, axis=part.mag_moment_unit_vec)  # What direction is the particle oriented to.  Where is the magnetic moment vector pointing.
        #part.mag_di_moment_arrow = arrow(pos=part.pos, axis=(0,1,0), length=part.radius*2, color=(1,1,0.25), opacity=.2)  # Where is the total force pushing us?
        part.angular_velocity = 0.0

        if part.this_type_count == 0 :
            #part.initial_wave_height_in_terms_of_dt = random.random() * part.period
            part.initial_wave_height_in_terms_of_dt = 0
        else :
            # handle pauli exclusion principle
            # Create particles at different times in their cycles.
            previous_part_of_same_type = wl.pars[i-1]
            part.initial_wave_height_in_terms_of_dt = ( part.period / num_particles_of_same_type ) + previous_part_of_same_type.initial_wave_height_in_terms_of_dt
        part.lorentz_factor = 1.0               # http://en.wikipedia.org/wiki/Lorentz_transformation
        part.rela_mass = part.rest_mass_kg      # Set relativistic mass
        part.mag_velocity = mag(part.velocity)
        part.old_dist_apart = 1.0               # Initialize to huge value.
        print " ", part.ptype_name,
        print "position ", part.pos
        print "\t velocity %1.1e" % mag(part.velocity), part.velocity
        print "\t wl.dt offset random ", part.initial_wave_height_in_terms_of_dt
        part.trail = curve(color=part.color, radius = part.radius / 8, pos=part.pos  )
        #part.trail.append(pos=part.pos,retain=50)         # 50 = arbitrary
        part.old_vel      = vector(0,0,0)
        part.old_pos      = vector(0,0,0)
        part.old_f_unit_v = vector(0,0,0)
        part.getting_further_apart = False
        prev_part = part
        part.getting_closer     = False
        part.old_getting_closer = False
        part.previous_force_was_3_digits = False
        part.label = label(pos=part.pos,text='%d' % part.index)
        part.moment_of_inertia = ( 2.0 / 5.0 ) * part.rest_mass_kg * ( part.radius ** 2 )   # http://en.wikipedia.org/wiki/List_of_moments_of_inertia
                                                                    # I modeled it as a sphere but perhaps it is close to a disc.
        print " part.moment_of_inertia", part.moment_of_inertia
    #calc_energy_of_system()
    #wl.orig_total_energy = wl.total_energy
    #wl.energy_limit = wl.orig_total_energy * 1.1
    #print "\t total energy of system %5g" % wl.orig_total_energy
    #wl.prev_total_energy = wl.orig_total_energy
    #wl.energy_ratio = 1.0           # Energy ratio from original
    # end of setup()

def display_particle_labels() :
    for part in wl.pars :
        part.label.visible = True
        part.label.pos = part.pos
    wl.display_particle_labels_bool = True
    num_loops_to_display_part_index = 25 * wl.particle_count
    wl.hide_number_label_count = wl.main_loop_count + num_loops_to_display_part_index

def hide_particle_labels() :
    for part in wl.pars :
        part.label.visible = False
    kb_input_label.visible = 0
    wl.display_particle_labels_bool = False


def calc_energy_of_system():
    wl.KE_system = 0
    wl.PE_system = 0
    for part in wl.pars :
        part.kinetic_energy = .5 * part.rest_mass_kg * mag2(part.velocity);         # mag2 = magnitude ^ 2
        wl.KE_system += part.kinetic_energy
        # Electric potential energy          http://en.wikipedia.org/wiki/Electric_potential_energy
        # Simplifying assumptions:
        # 1. The system has zero net charge
        # 2. The systems center of mass and center of charge is at the origin of the axis.  (0,0,0)
        # 3. If a particle is away from the center than potential energy is the amount of energy from current position to center, for single charge.
        #     Because charge at center is zero when all particles are there.
        # 4.
        # gravitational potential energy = U = mgh
        # For comparison purposes
        # pe = k*(q1*q2)/r
        mag_pos = mag(part.pos)
        if mag_pos == 0 :
            part.electric_energy = 0
        else :
            # Doesn't work.  What happens when part is inside the base_electric_pot_energy?  Value is negative.
            part.electric_energy = base_electric_pot_energy - ( Ke * ( elementary_charge ** 2 ) / mag_pos )
        wl.PE_system += part.electric_energy
    wl.total_energy = wl.KE_system + wl.PE_system
    print " PE", wl.PE_system, ", KE", wl.KE_system, ", tot", wl.total_energy

def update_particle_em_wave_heights(part):
    cycle            = ( wl.total_time + part.initial_wave_height_in_terms_of_dt ) * part.frequency * 2 * pi
    part.charge_cur  = ( part.charge_avg          * sin(cycle) ) + part.charge_avg    # Theory is that it oscillates around the average
    part.B_field_mag = ( part.avg_magnetic_moment * cos(cycle) ) + part.avg_magnetic_moment
    part.B_from_qv          = K_m * part.charge_cur * part.velocity                      # http://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law#Point_charge_at_constant_velocity
    if mag2(part.intrinsic_magnetic_moment) == 0 :
        if mag2(part.previous_magnetic_moment) == 0 :
            part.mag_moment_unit_vec = vector(0,0,-1)        # arbitrary
        else :
            part.mag_moment_unit_vec = norm(part.previous_magnetic_moment)
    else :
        part.mag_moment_unit_vec  = norm(part.intrinsic_magnetic_moment)
    part. previous_magnetic_moment = part.intrinsic_magnetic_moment                      # http://en.wikipedia.org/wiki/Spin_%28physics%29#Magnetic_moments
    part.intrinsic_magnetic_moment = part.mag_moment_unit_vec * part.B_field_mag         # Just using random direction for this case.
    part.intrinsic_B_field_plus_qv = part.intrinsic_magnetic_moment + part.B_from_qv
    if debug_verbosity > 25 :
        print " n ", part.num, part.ptype_name, 
        #print "b", part.B_field_mag, " uv", part.mag_moment_unit_vec,
        print " qv ", part.B_from_qv,
        print " intrinsic B field", part.intrinsic_magnetic_moment
        print " intrinsic_B_field_plus_qv", part.intrinsic_B_field_plus_qv

    #if (part.index != 0 and wl.main_loop_count&0x0f) or debug_sum_forces_on_one_particle > 25:
    #if (wl.main_loop_count&0xff == 0) or debug_sum_forces_on_one_particle > 25:
    if debug_sum_forces_on_one_particle > 25:
        percent_diff = int((part.charge_cur / part.charge_avg)*100)
        print " %%%3d" % percent_diff,
        #    print "\t charge_cur %+1.1e" % part.charge_cur, "  delta %%%d" % percent_diff, "  index ", part.index
    #else
    #    print "\t charge_cur %+5e" % part.charge_cur, "  delta %%%5.3f" % (part.charge_cur / part.charge_avg),   \
    #        "  index ", part.index

    #        if part.index != 0 and debug_sum_forces_on_one_particle > 50:
    #            print "charge_cur %+5e" % part.charge_cur, "  delta %5.3f" % (part.charge_cur / part.charge_avg), "  index ", part.index
    #    print "charge_cur delta from part average ", (part.charge_avg - part.charge_cur), "  index ", part.index


    #################################################################################
    # Might be documented at this web page:
    # http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/index.htm
    # http://www.wikihow.com/Find-the-Angle-Between-Two-Vectors
    # There are just the first two web pages that came when searching on this topic
def angle_between_two_norm_vectors(n1,n2,v1,v2,caller) :
    dot_product = dot(n1,n2)
    #print " dot ", dot_product
    if dot_product > 1.0 or dot_product < -1.0 :
        result = 0
        if wl.main_loop_count > 0 :
            #print " ref_pos ", ref_pos
            #print " pos1 ", pos1
            #print " pos2 ", pos2
            print caller
            print " v1 ", v1, n1
            print " v2 ", v2, n2
            print " dot ", dot_product
            #print " zero angle "
    else :
        # print "dot product", dot_product
        result = math.acos(dot_product)
        #print " angle ", result
    return result

def angle_between_two_vectors(v1,v2,caller) :
    return angle_between_two_norm_vectors(norm(v1),norm(v2),v1,v2,caller)


def angle_change_between_two_particles(pos1,pos2,ref_pos) :
    if pos1 == pos2 :
        result = 0.0
    else :
        v1 = pos1 - ref_pos
        v2 = pos2 - ref_pos
        result = angle_between_two_vectors(v1,v2,"Angle change")
    return result

def sum_forces_on_one_particle(part):
    # part = particle.  A pointer to the sphere object which holds most of our data.
    wl.largest_ele_force = 1e-15
    part.e_force_sum  = vector(0,0,0)       # Sum of electric forces acting on the particle
    part.B_force_sum  = vector(0,0,0)       # Sum of magnetic forces acting on the particle
    part  .force_sum  = vector(0,0,0)       # Sum of all forces acting on the particle
    part.closest_part = part
    part.closest_dist = 1                   # 1 meter apart
    has_printed = False
    #if part.is_electron :
    if ( debug_sum_forces_on_one_particle > 20 ):
        print " #%d" % part.index,
    for j in range(wl.particle_count):          # Calculate the force each particle exerts on this one
        if ( part.index == j ):
            continue                        # No need to calculate force that particle exerts on itself.
        if ( debug_sum_forces_on_one_particle > 60 ):
            print "part.index %d, j %d" % (part.index, j)
        other_part = wl.pars[j]
        dist_vec = part.pos - other_part.pos                    # dist is a vector
        mag2_dist       = mag2(dist_vec)                        # mag2 = magnitude squared
        part.dist_apart = mag (dist_vec)
        if part.dist_apart < part.closest_dist :
            part.closest_part = other_part                  # used to calculate B field direction
            part.closest_dist = part.dist_apart
        if j == 0 :
            if ( debug_sum_forces_on_one_particle > 25 ):
                print " dis a %1.0e" % part.dist_apart,
            if ( debug_sum_forces_on_one_particle > 35 ):
                print " delta %1.0e" % (part.dist_apart-part.old_dist_apart),
        # The below part is used to track when an electron is tossed by a proton.
        # Tossed because of digitization simulation error.
        # The digitization problem is corrected in implement_particle_calcs()
        if part.closest_part == other_part :
            part.getting_further_apart = (part.dist_apart > part.old_dist_apart)    # Don't know if these attributes are used.
            part.old_getting_closer    = part.getting_closer
            part.    getting_closer    = not part.getting_further_apart
            if part.getting_closer :
                if part.old_getting_closer :
                   part.getting_closer_count += 1
                else :
                   part.getting_closer_count  = 1

            # http://en.wikipedia.org/wiki/Electric_force
            # positive force is repulsive, negative value is attraction
            # For info on mag see: http://vpython.org/webdoc/visual/VisualIntro.html or http://vpython.org/webdoc/visual/vector.html
        dist_unit_vector = norm(dist_vec)
        e_force = dist_unit_vector * (Ke * part.charge_cur * other_part.charge_cur / mag2_dist)
        if ( debug_sum_forces_on_one_particle > 20 ):
            print "ele_force %3.0f " % mag(e_force),
        part.e_force_sum += e_force
        if mag(e_force) > wl.largest_ele_force :
            wl.largest_ele_force = mag(e_force)

        #################################################################
        # Calculate the magnetic force due to the some of the B fields. #
        #################################################################

        # Want to sum all of the forces on the moving charge
        # Force from a point charge: http://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law#Point_charge_at_constant_velocity
        # http://maxwell.ucdavis.edu/~electro/magnetic_field/pointcharge.html
        #
        # The B field of the other part relative to this one:
        # Parts have two sources of magnetic field
        # 1. Magnetic field created from moving charge
        # 2. Magnetic field from standing EM wave.  In this simulation also called intrinsic B field.
        #
        # What do the field lines look for the different magnetic fields?
        # 1. Magnetic field do to a moving charge are circular similar to those from a current in a wire:
        #       http://en.wikipedia.org/wiki/File:Manoderecha.svg
        #       http://www.youtube.com/watch?v=eWaA9RiLsno
        #       http://web.mit.edu/jbelcher/www/java/part_biot/part_biot.html
        #       http://www.physics.sjsu.edu/becker/physics51/mag_field.htm
        # The magnetic field created by the other part is circular
        #   so need to calculate what its field is at part location.
        # Using the right hand rule with the direction of motion.
        #

        #B_force = K_m * cross(other_part.intrinsic_B_field_plus_qv,-dist_unit_vector/mag2_dist)
        # Below is a form of the biot-savart law: http://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law#Point_charge_at_constant_velocity
        b_field_from_other_part_at_this_point_in_space = K_m * cross(
            other_part.intrinsic_B_field_plus_qv, -dist_unit_vector/mag2_dist )
        B_force = cross(b_field_from_other_part_at_this_point_in_space, part.intrinsic_B_field_plus_qv)
        if ( debug_sum_forces_on_one_particle > 25 ):
            if part.index >= 0 :
                print "\t type", part.ptype_name, " dist apart", part.dist_apart
                print "\t intrinsic_B_field_plus_qv ", part.intrinsic_B_field_plus_qv, " other ", other_part.intrinsic_B_field_plus_qv
                #print "\t b_field_from_other_part_at_this_point_in_space %1.0e %1.0e %1.0e mag %1.1e" % (
                #        b_field_from_other_part_at_this_point_in_space.x,
                #        b_field_from_other_part_at_this_point_in_space.y,
                #        b_field_from_other_part_at_this_point_in_space.z,
                #        b_field_from_other_part_at_this_point_in_space.mag )
                print "\t qv ", part.B_from_qv, " mag", mag(part.B_from_qv)
                print "\t other_part.intrinsic_magnetic_moment", other_part.intrinsic_magnetic_moment
                print "\t B x%1.0e y%1.0e z%1.0e %g" % (B_force.x,B_force.y,B_force.z, mag(B_force) )
        part.B_force_sum += B_force

    part.B_force_sum+= background_magnetic_field
    part  .force_sum = part.e_force_sum + part.B_force_sum
    part.B_force_mag = mag(part.B_force_sum)
    part.e_force_mag = mag(part.e_force_sum)
    part.t_force_mag = mag(part.force_sum)
        # Calculate the torque.
        # magnetic fields tend to align, so the torque and consequently the rotation will be to align in the same direction.
        # The rotation will be from current location that magnetic moment vector is pointing to the magnetic moment vector of other particle
        # The axis will be the cross product from old to new direction that the vector is pointing to.
        # Our torque vectors direction is pointing in the direction of the axis of rotation.
        # Torque  http://en.wikipedia.org/wiki/Magnetic_moment#Effects_of_an_external_magnetic_field_on_a_magnetic_moment
    part.torque     = cross(part.intrinsic_magnetic_moment,part.B_force_sum)
    part.torque_mag = mag ( part.torque )
    #print "\t torque x%1.0e y%1.0e z%1.0e %g " % (part.torque.x,part.torque.y,part.torque.z, mag(part.torque))

    if part.is_electron :
        part. EF_pointer.axis = norm(part.e_force_sum) * min_ptr_len
        part. MF_pointer.axis = norm(part.B_force_sum) * min_ptr_len
        part.tot_pointer.axis = norm(part.force_sum  ) * min_ptr_len
        #print scene.scale,
        # Make length relative to strongest force
        strongest_force = part.e_force_mag
        if part.B_force_mag > strongest_force :
            strongest_force = part.B_force_mag
        if part.t_force_mag > strongest_force :
            strongest_force = part.t_force_mag
        if strongest_force != 0.0 :
            part. EF_pointer.length = min_ptr_len + (max_ptr_len * ( part.e_force_mag / strongest_force ) )
            part. MF_pointer.length = min_ptr_len + (max_ptr_len * ( part.B_force_mag / strongest_force ) )
            part.tot_pointer.length = min_ptr_len + (max_ptr_len * ( part.t_force_mag / strongest_force ) )
            #if wl.main_loop_count > 10 :
            #    print "e%2.0f b%2.0f t%2.0f" % (part.e_force_mag, part.B_force_mag, part.t_force_mag),
    if debug_sum_forces_on_one_particle > 20 :
        c_vel   = ( part.mag_velocity / c ) * 100
        #print "e %3.0f b %3.0f" % (part.e_force_mag, part.B_force_mag),
        if part.is_electron :
            print " v%3.0f" % (c_vel),
        else :
            print " v%1.0f" % (c_vel),
        if debug_sum_forces_on_one_particle > 20 :
            if part.previous_force_was_3_digits or part.e_force_mag >= 100.0:
                if part.e_force_mag >= 100.0 :
                    part.previous_force_was_3_digits = True
                    part.force_3_digits_iteration = wl.main_loop_count
                elif (part.force_3_digits_iteration + 50) < wl.main_loop_count:
                    part.previous_force_was_3_digits = False
                print "e%3.0f b%3.0f t%3.0f" % (part.e_force_mag, part.B_force_mag, part.t_force_mag),
            else :
                print "e%2.0f b%2.0f t%2.0f" % (part.e_force_mag, part.B_force_mag, part.t_force_mag),

    #if part.is_proton :
    if ( debug_sum_forces_on_one_particle > 20 ):
        percent_diff = int((part.charge_cur / part.charge_avg         )*100)
        print "%%%3d" % percent_diff,
        percent_diff = int((part.B_field_mag/ part.avg_magnetic_moment)*100)
        print "B%%%3d" % percent_diff,
            #print "d %1.1e "    % part.dist_apart,
                #print " e for %2.0f %1.0e, b for %2.0f %1.0e," % (part.e_force_mag, part.e_force_mag, part.B_force_mag, part.B_force_mag),
                #print "e %1.0e b %1.0e," % (part.e_force_mag, part.B_force_mag),
            #print " e for %3.0f, b for %3.0f " % (part.e_force_mag, part.B_force_mag), " tot ", part.force_sum
            #    print "\t part.e_force_mag %3g" % part.e_force_mag
        #print "e_force_sum %e" % mag(part.e_force_sum), part.e_force_sum
    if part.e_force_sum > 100 :
        part.was_very_close = True
    else :
        part.was_very_close = False
    # end sum_forces_on_one_particle()

def update_dt():
    # If the electrical force is strong I want to have a smaller dt so that we can model orbits better
    # Also have smaller digital energy gain issue.
    new_dt = wl.large_dt / wl.largest_ele_force
    if new_dt < wl.dt :
        if (wl.last_dt_change / new_dt ) > 1.2 :
            print "\t lowering dt from %1.2e " % wl.last_dt_change, "to %1.2e" % new_dt, "e force", int(wl.largest_ele_force)
            wl.last_dt_change = new_dt
        wl.dt = new_dt
    elif new_dt > wl.dt and wl.dt < wl.desired_dt :
        if new_dt > wl.desired_dt :
            new_dt = wl.desired_dt
            print "\t Reached desired delta time (dt)"
        if (wl.last_dt_change / new_dt ) > 0.8 :
            print "\t lowering dt from ", wl.last_dt_change, "to", new_dt, "e force", wl.largest_ele_force
            wl.last_dt_change = new_dt
        wl.dt = new_dt
    

def calc_particle_acceleration(part):
    # Find the angle of the force on the current velocity.
    # This is used for special relativity calcs.     http://en.wikipedia.org/wiki/Special_relativity#Force
    radian_angle = angle_between_two_norm_vectors(part.vel_norm,part.force_unit_v,
                                                  part.velocity,part.force_sum   ,"calc accel")
    #degrees = radian_angle * 57.2957795
    #if part.is_electron :
    #    print "%3.0f" % degrees,
    # force = gamma**3 * mass * acceleration_parallel + gamma * mass * acceleration_perendicular
    # gamma = lorentz factor                        # http://en.wikipedia.org/wiki/Lorentz_factor
    # http://en.wikipedia.org/wiki/Special_relativity#Force
    # Need to divide the force into perpendicular and parallel parts.
    # If force is completely parallel      then force = gamma**3 * mass * accel
    # If force is completely perpendicular then force = gamma    * mass * accel
    if radian_angle == 0 :
        part.acceleration = part.force_sum / ( part.lorentz_factor**3 * part.rest_mass_kg )
    else :
        r = cross( part.velocity, part.force_sum )
        perpen_to_velocity = norm(cross(r, part.velocity))              # This is perpendicular to velocity on the same place as the force and velocity.
        perpen_force_fraction = math.fabs(sin(radian_angle))**2         # If this is 90 degrees then answer is 1
        parall_force_fraction = math.fabs(cos(radian_angle))**2         # If force is parallel to velocity then result is 1
        part.perpen_force = perpen_to_velocity * perpen_force_fraction * part.t_force_mag
        #if part.is_electron :
            #print "%1.3f %1.3f" % ( parall_force_fraction, perpen_force_fraction ),
            #print "ll %1.3f per %1.3f lorentz %g" % ( parall_force_fraction, perpen_force_fraction, part.lorentz_factor ),
        part.accel_parall = ( part.force_sum * parall_force_fraction ) / ( part.lorentz_factor**3 * part.rest_mass_kg )
        part.accel_perpen =   part.perpen_force                        / ( part.lorentz_factor    * part.rest_mass_kg )
        #part.acceleration = part.force_sum / ( part.lorentz_factor * part.rest_mass_kg )
        #part.acceleration = part.force_sum / part.rest_mass_kg
        part.acceleration = part.accel_perpen
        #if part.velocity < c or degrees > 90:
        part.acceleration += part.accel_parall
    # calc_particle_acceleration()


def calc_new_particle_position(part):
    dt_changed = False
    part.force_unit_v = norm(part.force_sum)
    calc_particle_acceleration(part)

    if part.index != 0 :
        if ( debug_sum_forces_on_one_particle > 50 ):
            print "\t accelerat %e " % mag(part.acceleration), part.acceleration
        if ( debug_sum_forces_on_one_particle > 60 ):
            if ( part.velocity == vector(0,0,0) or part.velocity == 0 ):
                print "\t old direction  zero "
            else:
                print "\t old direction  ", part.velocity
                print "\t old direction  ", norm(part.velocity)             # norm = unit vector, http://vpython.org/webdoc/visual/vector.html
    part.new_velocity = part.velocity + (part.acceleration * wl.dt)         # velocity = v0 + acceleration * wl.dt
    if part.index != 0 and debug_sum_forces_on_one_particle > 60 :
        print "new v", part.new_velocity, " old", part.velocity, " accel", part.acceleration
    part.new_mag_vel  = mag (part.new_velocity)
    part.distance_traveled = part.new_velocity * wl.dt                      # used to set new position
    part.mag_dist_traveled = mag(part.distance_traveled)
    if ( debug_sum_forces_on_one_particle > 30 ):
        if part.index != 0 :
            print " v %1.1e " % mag(part.new_velocity),
        if part.index == 0 :
            if ( debug_sum_forces_on_one_particle > 60 ):
                print part.new_velocity,
    if ( debug_sum_forces_on_one_particle > 50 ):
        print "  dist trav ", part.mag_dist_traveled
        #print part.new_velocity
    part.new_pos = part.pos + part.distance_traveled
    if ( debug_sum_forces_on_one_particle > 15 and wl.main_loop_count%64 == 0):
        print "p %1.0e v%1.0e" % ( part.new_pos.mag, part.velocity.mag ),
    if part.is_electron :
        if part.new_velocity > wl.fastest_electron :
            wl.fastest_electron = part
    else :
        if part.new_velocity > wl.fastest_proton :
            wl.fastest_proton = part
    
    return dt_changed
    # end of calc_new_particle_position()

# Todo:
#   Remove rotate particle.  Want to see curve with particles
#   Have dt change while part is tossed


def rotate_particle(part):
    # Have torque.  Want to spin the particle.  How do i describe the current orientation?
    # The orientation is described by the vector: magnetic moment.
    # How do I spin it?
    # Something like angular moment + torque
    # angular acceleration = torque / moment of inertia     # http://en.wikipedia.org/wiki/Angular_acceleration
    # Similar to acceleration = force / mass  or f = m * a
    # moment_of_inertia for sphere is: 2/5 m r^2             # http://en.wikipedia.org/wiki/List_of_moments_of_inertia

    part.angular_accel = part.torque_mag / part.moment_of_inertia
        # for linear motion: velocity = v0 + at
    part.angular_velocity = part.angular_velocity + part.angular_accel * wl.dt
        # for linear motion: location = x0 + v0*t + .5at^2
        # v2 = rotate(v1, angle=theta, axis=(1,1,1))   http://vpython.org/contents/docs/visual/vector.html
    if debug_verbosity > 20 :
        print " n ", part.num,
        print " intrinsic_magnetic_moment", part.intrinsic_magnetic_moment,
        print " angular_velocity", part.angular_velocity,
        print " torque", part.torque,
       # Minor documentation on rotate: http://vpython.org/webdoc/visual/vector.html
    part.intrinsic_magnetic_moment = rotate(part.intrinsic_magnetic_moment, angle=part.angular_velocity * wl.dt, axis=part.torque)
    part.mag_moment_unit_vec       = rotate(part.mag_moment_unit_vec      , angle=part.angular_velocity * wl.dt, axis=part.torque)
    if debug_verbosity > 20 :
        print " mag_moment_unit_vec ", part.mag_moment_unit_vec,
        print " axis", part.mm_pointer.axis
    part.mm_pointer.axis     = part.mag_moment_unit_vec * part.radius*2;



def do_particle_calcs():
    """Calculate the forces acting on each particle and update positions"""
    #global wl.part_with_most_energy, wl.curr_KE_system, debug_sum_forces_on_one_particle
    wl.curr_KE_system = 0.0
    wl.part_with_most_energy = wl.pars[0]              # initialize
    wl.fastest_proton = wl.pars[0]
    wl.fastest_electron = wl.pars[1]
    #wl.part_with_highest_velocity = wl.pars[0]
    for part in wl.pars :
        sum_forces_on_one_particle(part)
    #update_dt()
    for part in wl.pars :
        dt_changed = calc_new_particle_position(part)
        if dt_changed:
            break
        rotate_particle(part)
    return dt_changed
    # end of do_particle_calcs()


def vector_copy(p1,p2):     # copy by value instead of copy by reference
    p1.x = p2.x
    p1.y = p2.y
    p1.z = p2.z

def implement_particle_calcs():
    global debug_sum_forces_on_one_particle
    highest_vel_part = wl.pars[0]
    for i in range(wl.particle_count):
        decrease_vel  = 0
        part = wl.pars[i]
        vector_copy(part.old_f_unit_v, part.force_unit_v)
        vector_copy(part.old_vel     , part.    velocity)
        part.vel_norm = norm(part.new_velocity)
            #  10. If particle is outside of the normal size of an atom then assume it is escaping due to energy gain.
            #       reset kinetic energy to zero.
        escape = False
        if part.is_electron :
            #if part.pos.mag > 350e-12 :                 # Some atoms are 260e-15 in size. http://en.wikipedia.org/wiki/Atomic_radius
            if part.pos.mag > wl.particle_count*1e-12 :                 # Some atoms are 260e-15 in size. http://en.wikipedia.org/wiki/Atomic_radius
                part.new_pos = 1e-13 * part.vel_norm
                escape = True
        elif part.pos.mag > (wl.particle_count*1e-14) :                # http://en.wikipedia.org/wiki/Atomic_nucleus
            part.new_pos = 1e-15 * part.vel_norm
            escape = True
        if escape :
            print "\t escape " , part.num, " vel", part.new_velocity.mag
            part.new_velocity = vector(0,0,0)
            #debug_sum_forces_on_one_particle += 5
        mag_velocity  = mag (part.new_velocity)
        """
        if part.getting_further_apart and part.mag_velocity > (c*(7/8)):
            #if part.closest_dist < size_of_hydrogen_nucleus and part.e_force_mag > 150 :  # If we do less than 10 then may have issue when electron orbits proton.
                #if part.
                #part.has_been_tossed = True
                #rate(20)
            print "\t slowing part", part.num, part.ptype_name, part.closest_part.num, part.closest_part.ptype_name,
            if ( part.new_velocity.mag > c/2 ) :
                part.new_velocity = part.vel_norm * c/2
                part.mag_velocity = c/2
            # slow down the proton too
            if part.closest_part.velocity.mag > c/4 and part.closest_part.is_proton:
                part.closest_part.    velocity = part.closest_part.vel_norm * c/4
                part.closest_part.new_velocity = part.closest_part.vel_norm * c/4
                part.mag_velocity = c/4
        """
        mag_velocity  = mag (part.new_velocity)
        vector_copy(part.velocity, part.new_velocity)
        if part.mag_velocity > highest_vel_part.mag_velocity:
            highest_vel_part = part
        # Handle special relativity
        part.vel2_div_c2 = part.velocity.mag2 / c2
        if part.vel2_div_c2 < 1.0 :
            divisor = math.sqrt(1.0 - part.vel2_div_c2)           # http://en.wikipedia.org/wiki/Lorentz_factor
            part.lorentz_factor = 1.0 / divisor                 # http://en.wikipedia.org/wiki/Lorentz_transformation
            #if part.is_electron :
                #print "vel %2.0f, div %g, lorentz %g" % ( (part.mag_velocity/c)*100, divisor, part.lorentz_factor),
                #print "vel %2.0f, div %g, lorentz %g" % ( (part.mag_velocity/c)*100, divisor, part.lorentz_factor),
        else :
            print "Velocity matches c. :0(", part.num, part.ptype_name,
            part.velocity = norm(part.velocity) * (c-1)
            #exit(0)
            part.lorentz_factor = 9e30                  # Some random rediculuously large number
            #part.lorentz_factor = 1                  # Some random rediculuously large number
            #print "highest vel part num %d" % highest_vel_part.num,
        #part.old_dist_trav = part.distance_traveled
        vector_copy(part.old_pos ,part.    pos)
        vector_copy(part.    pos ,part.new_pos)
        if part.velocity != vector(0,0,0) :
            #part.vel_pointer.axis   = norm(part.velocity) * min_ptr_len
            part.vel_pointer.axis   = norm(part.velocity)
            part.vel_pointer.pos    = part.pos
            part.vel_pointer.length = 2 * max_ptr_len * ( part.mag_velocity / c )
            if part.vel_pointer.length < min_ptr_len :
                part.vel_pointer.length = min_ptr_len
        if part.is_electron :
            part. EF_pointer.pos = part.pos
            part. MF_pointer.pos = part.pos
            part.tot_pointer.pos = part.pos
        part.mm_pointer.pos = part.pos              # Magnetic moment
        #part.mag_di_moment_arrow.pos = part.pos
        #print " dist apart cur ", part.dist_apart, part.old_dist_apart
        part.old_dist_apart = part.dist_apart

        # Set how long we want particle trails.
        # Could make simulation run faster by having fewer joints in the trails.
        if   part.mag_velocity > (c/2) :
            and_value = 8
        elif part.mag_velocity > (c/4) :
            and_value = 16
        elif part.mag_velocity > (c/8) :
            and_value = 32
        elif part.mag_velocity > (c/16) :
            and_value = 64
        elif part.mag_velocity > (c/32) :
            and_value = 128
        else :
            and_value = 256
        if ( wl.main_loop_count % and_value ) == 0 :
            my_retain = 64         # arbitrary, what I think looks nice
            if part.is_electron :
                my_retain = 512 / wl.particle_count        # arbitrary
            if my_retain < 8:
                my_retain = 8
            part.trail.append(pos=part.pos,retain=my_retain)
        #part.trail.append(pos=part.pos)                  # If using earlier version of vpython
    #print "wl.pars[1].pos-old_pos", (wl.pars[1].pos-wl.pars[1].old_pos),
    #print "high vel pos delta", pos_delta
    # end implement_particle_calcs()


def particle_calcs():
    wl.old_dt = wl.dt
    scene.autoscale = 0
    while True:
        dt_updated = do_particle_calcs()
        if not dt_updated:
            break           # We didn't update our delta time (wl.dt) so we are done with do_particle_calcs.
    implement_particle_calcs()
    scene.autoscale = wl.auto_zoom
    wl.prev_KE_system = wl.curr_KE_system
    #if ( debug_sum_forces_on_one_particle > 20 ):
        #if wl.main_loop_count%4 == 0 :

def usage() :
    print "particles.py <num particles>"
    print "\t a = toggle auto zoom"
    print "\t n = show numbers"
    print "\t <space bar> = toggle between pause and play mode"

def main(argv):
    try:
        opts, args = getopt.getopt(sys.argv[1:], "ho:v", ["help", "output="])   # http://docs.python.org/library/getopt.html
    except getopt.GetoptError, err:
        # print help information and exit:
        print str(err) # will print something like "option -a not recognized"
        usage()
        sys.exit(2)
    output = None
    verbose = False
    for o, a in opts:
        if o == "-v":
            verbose = True
        elif o in ("-h", "--help"):
            usage()
            sys.exit()
        elif o in ("-o", "--output"):
            output = a
        else:
            assert False, "unhandled option"
    argc = len(argv)
    if ( argc > 0 ) :
        wl.particle_count = int(argv[0])
    setup_workload()
    setup_particles()
    display_particle_labels()
    if ( wl.particle_count <= 2 ) :
        hide_particle_labels()
    while true:
        if scene.kb.keys: # did the user hit a key?
            key = scene.kb.getkey() # obtain keyboard information
            #print "\ns%s\n", key
            if len(key) == 1:
                kb_input_label.text = key # append new character
                kb_input_label.visible = 1
                if key == 'q' or key == 'e':
                    print "Exiting"
                    sys.stdout.flush()
                    break
                elif key == 's' or key == 'p' or key == ' ':   #stop, pause or play
                    wl.play_mode = not wl.play_mode
                    if wl.play_mode :
                        kb_input_label.visible = 0
                elif key == 'n' or key == 'N' :
                    if wl.display_particle_labels_bool :
                        hide_particle_labels()
                    else :
                        display_particle_labels()
                elif key == 'a' or key == 'A' :
                    wl.auto_zoom = not wl.auto_zoom
            else :
                kb_input_label.visible = 0
                kb_input_label.text = '' # erase all the text
        if wl.play_mode :
            area_of_interest = 50000000
            #area_of_interest = 50
            if   wl.main_loop_count > area_of_interest + 118 :
                rate(2)
            elif wl.main_loop_count > area_of_interest + 100:
                rate(5)
            elif wl.main_loop_count > area_of_interest + 50:
                rate(10)
            elif wl.main_loop_count > area_of_interest + 25 :
                rate(25)
            elif wl.main_loop_count > area_of_interest :
                rate(50)            # 10 = 1/10 of a second wait.  http://vpython.org/webdoc/visual/rate.html
            #rate(2)
            #rate(6500/(wl.main_loop_count+1))            # http://vpython.org/webdoc/visual/rate.html
            if ( debug_sum_forces_on_one_particle > 20 or 
               ( debug_sum_forces_on_one_particle > 15 and wl.main_loop_count%64==0 )):
                print
                sys.stdout.flush()
                print "it %4d" % wl.main_loop_count,
            for i in range(wl.particle_count):
                part = wl.pars[i]
                update_particle_em_wave_heights(part)        # Update the current magnitude of the e and m wave.

            particle_calcs()
            wl.total_time += wl.dt
            wl.main_loop_count += 1
            #rate(1)                 # Remove this after debug
            if wl.display_particle_labels_bool and wl.main_loop_count > wl.hide_number_label_count :
                hide_particle_labels()
        else :
            rate(100)

if __name__ == "__main__":
    main(sys.argv[1:])

# This theory explains the following:
# 1) All mass as standing EM waves
#     Electric charge as an oscillation of electric wave between 0 and 2e.  Average is e.
#     There are no particles, just standing EM waves or from a different perspective: particles are standing EM waves.
#     
# 2) Physics spin as the amplitude of the magnetic wave
#   If you look close at the stern-gerlach experiment.  If you get exact measurements for the particle distribution you will
#   see the result is sinusoidal distribution not a split.
#   This also explains other results from Japan that should the electrons flowing sinusoidally.
# 3) Weak nuclear force
#   As electrical and magnetic forces and particles constantly in motion in the nucleus
#   Nucleus not a rigid body.  Motion creates strong magnetic forces.
#   Explains beta decay.
# 4) Why a neutron is not stable
#   Electron is hurling around the proton at relativistic speeds
#   Both particles oscillate between 0 and 2e charge.
#   When electron and proton are close to zero charge the two particles separate.
# 5) Explain spin as intrinsic magnetic field
#   A changing electric field will give rise to an apposing magnetic field.
#
# May help explain strong nuclear force:
#   No neutrons.  Just electrons and protons.
#   Neutron is a proton and electron superimposed on each other.
#   a. Electron helps mask charge of nearby proton
#       When electron at 2e charge will cause protons to come closer together
#       but when electron at zero charge then protons will fly away
#   b. Magnetic fields forces particles into circular path
#       Stops them from flying directly away from each other.
#   c. Orbiting electrons spend time near the nucleus helping nucleus be stable
#       Can fly right thru center and possibly bump other electron out
#   d. Magnetic force from all particles can combine and attract.
#       Doesn't have to oppose like protons positive repelling each other.
#       Magnetic force from qv and intrinsic magnetic field
#       Magnetic force from qv high because electron buzzing around protons at nearly the speed of light.
#       Two protons traveling in the same direction will have a magnetic force that draws them together.
#           If speed is high enough then perhaps strong compared to electric repulsive force.
#   e. When measure size of the nucleus, perhaps what is being measured is a strong magnetic field.
#      Perhaps a protons are further apart then believed but still contribute to magnetic field of
#      the nucleus.
#   
#   
#
# Questions:
#   Is split in Stern-Gerlach experiment due mostly to qv, intrinsic magnetic field, or both?
#
#
"""
to do:
    1. Fix Only reduce velocity of electrons that have been tossed by a proton and vice versa.
        Bring it back much sooner?
       Quit the experiment if it has escape velocity?
       Calculate escape velocity?

Enhance the program to?
   1. Specify opacity on command line and thru keyboard events.
   2. show number for particles
   3. Display info on a particular particle
   4. Specify what to log or not log anything
   5. Only show total force.  Force is colored for electrical or magnetic force.
   6. Specify trail length or no trails
   7. Specify wl.dt as a multiple of proton frequency of 1/16.
   8. Other things that can be done to improve speed.
   9. Initial spacing
   0. Specify trail length. retain=
   a. Different algorithms to compensate for energy gain or no algo.
   b. Different algorithms for how the magnetic field lines up.
"""
"""
Simulation to do:
    1. http://en.wikipedia.org/wiki/Helium-6#Table
        Does it predect helium isotope decay?
    2.
"""