#!/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 __future__ import division from visual import * # See vpython.org for more info. #from collections import deque import math import getopt import pprint #import Queue import random # http://effbot.org/pyfaq/how-do-i-generate-random-numbers-in-python.htm import sys # for argv, flush #import threading # Have computations on one thread and drawing on another import types #import copy print "version 2010.03.22" 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.num_particles = 2 # hydrogen #wl.num_particles = 3 # Neutron and proton wl.num_particles = 4 # heavy hydrogen, Deuterium #wl.num_particles = 6 # http://en.wikipedia.org/wiki/Tritium #wl.num_particles = 8 # Helium, 2 neutrons #wl.num_particles = 12 # lithium 7.5% #wl.num_particles = 14 # lithium 92.5% #wl.num_particles = 18 # Beryllium #wl.num_particles = 20 # Boron 20% #wl.num_particles = 22 # Boron 80% wl.num_particles = 24 # Carbon num_particles_to_atom_type = [ "1", "Hydrogen", "Neutron and proton", "heavy hydrogen, Deuterium", "5", "Tritium, hydrogen and two neutrons", "7", "Helium, 2 neutrons", "9", "10", "11", "Lithium 7.5%", "13", "lithium 92.5%", "15", "16", "17", "Beryllium", "19", "Boron 20%", "21", "Boron 80%", "23", "Carbon" ] atom_type_to_num_particles = {'hydrogen':2, 'heavy':4, 'deut':4, 'deuterium':4, 'trit':6, 'tritium':6, 'helium': 8, 'lithium': 14, 'beryllium':18, 'boron':22, 'carbon': 24 } debug_verbosity = 15 debug_sum_forces_on_one_particle = 15 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 # http://en.wikipedia.org/wiki/Proton 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. There is something wrong with this. The way I'm calculating this is wrong, but the frequency I think is correct. proton .frequency = ( 2 * proton.mass_eV ) / h_eVs # e = h*f , f = e /h, http://en.wikipedia.org/wiki/Pair_production 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 = 1.410607e-26 # http://en.wikipedia.org/wiki/Magnetic_moment#Elementary_particles # Doesn't make sense that a proton would have a small magnetic moment. # http://physics.nist.gov/cgi-bin/cuu/Value?mup|search_for=proton+magnetic+moment #electron.magnetic_moment *= 100 # 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. proton .magnetic_moment = electron.magnetic_moment * el_pr_weight_diff # What it should be #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. initial_velocity_max = 3e4 # arbitrary max_ptr_len = proton.radius * 7 min_ptr_len = proton.radius / 2 #vel_max_ptr_len = max_ptr_len / 2 proton .vel_max_ptr_len = proton.radius * 16 electron.vel_max_ptr_len = proton.radius * 8 min_part_opacity = 0.001 # Arbitrary. Just chose what I think looks good. len_part_opacity = 0.99 kb_input_label = label(visible=0) def scene_setup() : wl.size_adjust = wl.num_particles / 2 # Larger particles when there are more of them because they are further apart. # See setup_globals(). Just eye candy improvement. wl.size_adjust = 1 if wl.num_particles < 25 : scene.title = num_particles_to_atom_type[wl.num_particles-1] else : 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 # offset from top of your monitor scene.y = 0 # offset fromm left side of your monitor axis_length = size_of_hydrogen_nucleus * 4 * wl.size_adjust y_axis_label = label(pos=(0,axis_length,0), text="Y", opacity=.2, box=0, line=0) y_axis = arrow(pos=(0,0,0), axis=(0,size_of_hydrogen_nucleus,0), length=axis_length, opacity=0.3, 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=.2, box=0, line=0) x_axis = y_axis.__copy__() x_axis.axis = (size_of_hydrogen_nucleus,0,0) x_axis.color = (0,0.5,0.5) x_axis.length= axis_length #scene.show_rendertime = True def set_dt() : wl.desired_dt = proton.period / wl.num_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 num_dt_per_proton_period %d, wl.dt %1.1e" % ( wl.num_dt_per_proton_period, wl.dt ) sys.stdout.flush() wl.last_dt_change = wl.dt wl.max_force = 8 * wl.num_dt_per_proton_period / 32 # In newtons. Larger force means particles get thrown. Smaller force means less acurate. def workload_setup() : if wl.num_dt_per_proton_period == 0 : wl.num_dt_per_proton_period = 8 wl.num_dt_per_proton_period = 16 #wl.num_dt_per_proton_period = 32 wl.num_dt_per_proton_period = 64 #wl.num_dt_per_proton_period = 128 #wl.num_dt_per_proton_period = 256 #wl.num_dt_per_proton_period = 512 #wl.num_dt_per_proton_period = 1024 #wl.num_dt_per_proton_period = 2048 #wl.num_dt_per_proton_period = 4096 #wl.num_dt_per_proton_period = 8192 #wl.num_dt_per_proton_period = 8192*2 #wl.num_dt_per_proton_period = 8192*8 set_dt() wl.parts = [] # wl.parts = 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 wl.hide_objects_far_from_center = False wl.nucleus_radius = 16 * 1.25e-15 * ((wl.num_particles * 3.0/4.0) ** (1.0/3.0)) #http://en.wikipedia.org/wiki/Atomic_nucleus#Nuclear_models print "nucleus radius ", wl.nucleus_radius wl.zoom_radius = wl.nucleus_radius wl.num_linear_steps = 5 wl.proton_flash = wl.num_particles > ( 1024 / wl.num_dt_per_proton_period ) """ wl.proton_flash = True if wl.num_dt_per_proton_period >= 1024 : wl.proton_flash = wl.num_particles > 1 elif wl.num_dt_per_proton_period >= 512 : wl.proton_flash = wl.num_particles > 2 elif wl.num_dt_per_proton_period >= 256 : wl.proton_flash = wl.num_particles > 3 elif wl.num_dt_per_proton_period >= 128 : wl.proton_flash = wl.num_particles > 4 elif wl.num_dt_per_proton_period >= 64 : wl.proton_flash = wl.num_particles > 6 else : wl.proton_flash = wl.num_particles > 8 #if wl.num_particles == 2 """ wl.rotation = False wl.furthest_electron = simple_class() wl.furthest_electron.pos = vector(0,0,0) wl.furthest_electron.label = simple_class() 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 my_arrow(pos,color,opacity=0.3) : # part.EF_pointer = my_arrow(pos=part.pos, color=(1,0,0), opacity=0.2 ) return arrow(pos=pos, length=min_ptr_len, color=color, opacity=opacity, fixedwidth=1, shaftwidth=min_ptr_len*0.2) def put_equidistant_on_a_sphere(electrons_flag): # Thanks to Paul Bourke for the algorithm # http://local.wasp.uwa.edu.au/~pbourke/geometry/spherepoints/source1.c num_protons = int(round(wl.num_particles / 2.0 )) num_iterations = wl.num_particles * 32 if electrons_flag : i_start = num_protons + 1 i_end = wl.num_particles else : i_start = 0 i_end = num_protons num_iterations = num_iterations * 2 #print "start", i_start, "end", i_end and_value = num_protons * 32 radius = wl.nucleus_radius * 0.5 p1 = vector(0,0,0) p2 = vector(0,0,0) # How far to move two points that are two close? # Would be nice to know what is the ideal distance, or average distance move # two close parts that far apart. accuracy = 0.9 # Put the protons on the edge of a sphere #for i in range(0,num_protons) : for i in range(i_start,i_end) : wl.parts[i].pos.mag = radius for iteration_count in range(0,num_iterations) : # Find the closest two points and separate them some closest_1 = wl.parts[i_start] closest_2 = wl.parts[i_start+1] min_dist = mag(closest_1.pos - closest_2.pos) tot_dist = min_dist max_dist = min_dist num_comparisons = 1 for i in range(i_start,i_end-1) : part = wl.parts[i] for j in range(i+1,i_end) : distance = mag(part.pos - wl.parts[j].pos) #print i, j, distance tot_dist = tot_dist + distance num_comparisons += 1 if ( distance < min_dist ) : closest_1 = part closest_2 = wl.parts[j]; min_dist = distance if ( distance > max_dist ) : max_dist = distance error_dist = max_dist - min_dist avg_distance = tot_dist / num_comparisons distance_vector = closest_2.pos - closest_1.pos distance_vector.x = distance_vector.x + ((random.random()-0.5) * accuracy * error_dist) distance_vector.y = distance_vector.y + ((random.random()-0.5) * accuracy * error_dist) distance_vector.z = distance_vector.z + ((random.random()-0.5) * accuracy * error_dist) closest_2.pos = closest_1.pos + ( distance_vector * (1 + accuracy )); closest_1.pos = closest_1.pos - ( distance_vector * accuracy ); if accuracy > 0.15 : accuracy = accuracy - 0.05 elif accuracy > 0.01 : accuracy = accuracy - 0.005 #print "clos_1 ", closest_1.index, "clos_2", closest_2.index, #print " accuracy ", int(accuracy*100), #print " avg", avg_distance, #print " amount_to_move", amount_to_move if iteration_count % and_value == 0 : print "min_dist", min_dist, "max_dist", max_dist, "avg_distance", avg_distance sys.stdout.flush() closest_1.pos.mag = radius closest_2.pos.mag = radius print "min_dist", min_dist, "max_dist", max_dist, "avg_distance", avg_distance print def particles_setup(): initial_spacing = wl.nucleus_radius electrons_in_center = 0 electron_count = 0 proton_count = 0 num_particles_of_same_type = int ( (wl.num_particles / 2) + 0.5 ) want_electrons_between_protons = wl.num_particles > 4 # Sometimes electrons randomly appear next to each other then fly away fast causing issues. max_electrons_in_center = wl.num_particles / 4 for i in range(0,wl.num_particles): # For describing sub atomic particles particles we start with visual python sphere object then add to it. wl.parts.append(sphere(radius=electron.visual_radius)) # Specify a radius that gets reset so that we don't get flicker with autoscale. part = wl.parts[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.num_particles / 2.0 ) ; part.is_electron = not part.is_proton part.velocity = vector(0,0,0) # For info on vector see: http://vpython.org/webdoc/visual/vector.html spacing = initial_spacing #if i <= 2 : # part.pos = vector (1e-15,-1e-15+(3*i*1e-15),0) # part.velocity = vector(c*.99999999,0,0) # need_random_loc_assigned = False #else : need_random_loc_assigned = True if part.is_electron : # Put electron between protons. if i >= 2 and want_electrons_between_protons and electrons_in_center < max_electrons_in_center : electrons_in_center = electrons_in_center + 1 need_random_loc_assigned = False #part.pos = mid_point(previous_proton, prev_prv_proton) proton_index = int(i-(wl.num_particles/2.0)) print "proton index ", proton_index, p = wl.parts[proton_index] np = wl.parts[proton_index+1] # np = next proton if np.is_electron : if i > 4 : np = wl.parts[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 " electron in the middle of two protons", if need_random_loc_assigned : 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 #if i == 1 : # part.velocity = vector(c*0.9,0,0) # Make experiment go faster #else : part.vel_pointer = my_arrow(pos=part.pos, color=(0,1,0) ) # direction of velocity part.mag_moment_unit_vec = vector(0,1,0) # initial z direction is arbitrary if part.is_electron : part.p_type = "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(),random.random()) # Redish color, ^2 the other colors makes them closer to zero part.vel_max_ptr_len = electron.vel_max_ptr_len part.EF_pointer = my_arrow(pos=part.pos, color=(1,0,0) ) # Where is the electric force pushing us? part.MF_pointer = my_arrow(pos=part.pos, color=(0,0,1) ) # Where is the magnetic force pushing us? part.tot_pointer = my_arrow(pos=part.pos, color=(1,1,1) ) # Where is the total force pushing us? #part.mm_pointer = arrow(pos=part.pos, length=part.radius*2, color=(0.5,0.5,0.0), opacity=.1, axis=part.mag_moment_unit_vec) # What direction is the particle oriented to. Where is the magnetic moment vector pointing. #part.tot_pointer = simple_class() part.avg_magnetic_moment = electron.magnetic_moment; part.this_type_count = electron_count electron_count += 1; part.label = label(pos=part.pos,text='e%d' % part.index) else: part.p_type = "proton " part.charge_avg = proton.charge 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(),random.random(),1.0) # Blue-ish color, ^2 the other colors makes them closer to zero part.vel_max_ptr_len = proton.vel_max_ptr_len part. avg_magnetic_moment = proton.magnetic_moment; part.this_type_count = proton_count proton_count += 1; part.EF_pointer = simple_class() # Make logic for turning on and off visibility easier. part.MF_pointer = simple_class() #part.mm_pointer = simple_class() #part.mm_pointer.axis = " " part.tot_pointer = simple_class() part.label = label(pos=part.pos,text='p%d' % part.index) part.intrinsic_magnetic_moment = part.mag_moment_unit_vec * part.avg_magnetic_moment # http://en.wikipedia.org/wiki/Spin_%28physics%29#Magnetic_moments #part.tot_pointer = simple_class() part.mm_pointer = simple_class() part.mm_pointer.axis = " " #part.EF_pointer = my_arrow(pos=part.pos, color=(1,0,0) ) # Where is the electric force pushing us? #part.EF_pointer = simple_class() # Make logic for turning on and off visibility easier. #part.MF_pointer = my_arrow(pos=part.pos, color=(0,0,1) ) # Where is the magnetic force pushing us? #part.type = part.p_type #print "part.radius", part.radius # magnetic Dipole moment pointer. http://en.wikipedia.org/wiki/Magnetic_dipole_moment #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 phases in their cycles. previous_part_of_same_type = wl.parts[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.old_dist_apart = 1.0 # Initialize to huge value. 1.0 = 1 meter print " ", part.p_type, #print "position ", part.pos, #print "\t wl.dt offset random ", part.initial_wave_height_in_terms_of_dt #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.moving_further_apart = False prev_part = part part.getting_closer = False part.old_getting_closer = False part.previous_force_was_3_digits = False part.label.visible = False 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 part.accel_rq = [] # create results queue (rq). Used for Adams-bashforth numerical integration part.veloc_rq = [] for i in range(wl.num_linear_steps) : part.accel_rq.append(vector(0,0,0)) # prime our results queue (rq) part.veloc_rq.append(vector(0,0,0)) # prime our results queue (rq) print if wl.num_particles > 4 : put_equidistant_on_a_sphere(False) put_equidistant_on_a_sphere(True) for part in wl.parts: if part.is_electron : #print "\n \t period ", part.period, " initial_wave_height_in_terms_of_dt ", part.initial_wave_height_in_terms_of_dt, wave_segment = int(part.initial_wave_height_in_terms_of_dt / ( part.period / 16 )) print " %d wave_segment %d " % (part.index, wave_segment), # sine wave # electrons at 2e charge should be in center # 0 e charge out side # 0-2 = center # 3 - getting ready to leave # 4 - leaving # 5-6 - out # 7 approaching if wave_segment == 0 : part.pos.mag = part.pos.mag * 1.5 part.velocity = - norm(part.pos) * (c/3) # Should be heading away if electric force is on the down side. #if wl.num_particles == 2 : #part.pos.mag = part.pos.mag * 4 #else : # part.velocity = -norm(part.pos) * (c/4) # Heading towards center # part.pos.mag = part.pos.mag * 2 elif wave_segment < 6 : part.pos.mag = part.pos.mag * (1+(random.random()-.5)) print "center electron, ratio to nucleus radius %4.2f" % (part.pos.mag/wl.nucleus_radius) elif wave_segment < 8 : part.velocity = norm(part.pos) * (c/2.0) # Should be heading away if electric force is on the down side. elif wave_segment == 8 : part.pos.mag = part.pos.mag * 1.5 # Calibrated for helium part.velocity = norm(part.pos) * (c/4.0) # Should be heading away if electric force is on the down side. elif wave_segment == 9 : part.pos.mag = part.pos.mag * 1.5 # Lithium part.velocity = norm(part.pos) * (c/6) # Should be heading away if electric force is on the down side. elif wave_segment == 10 : part.pos = norm(part.pos) * wl.nucleus_radius * 1.1 # Tritium part.velocity = norm(part.pos) * (c/32) # Should be heading away if electric force is on the down side. elif wave_segment == 11 : part.pos = norm(part.pos) * wl.nucleus_radius * 1.3 # Lithium elif wave_segment == 12 : part.pos = norm(part.pos) * wl.nucleus_radius * 1.4 # Helium elif wave_segment == 13 : part.pos = norm(part.pos) * wl.nucleus_radius * 1.25 # Adjusted for lithium part.velocity = - norm(part.pos) * (c/8) # Should be heading away if electric force is on the down side. elif wave_segment == 14 : part.pos = norm(part.pos) * wl.nucleus_radius * 1.25 part.velocity = - norm(part.pos) * (c/3.5) # Should be heading away if electric force is on the down side. elif wave_segment < 16 : part.pos = norm(part.pos) * wl.nucleus_radius * 0.9 part.velocity = - norm(part.pos) * (c/3) # Should be heading away if electric force is on the down side. else : print "Something wrong with wave segment calc" sys.exit(1) print " velocity %1.1e , pos %5.1e" % (mag(part.velocity), part.pos.mag ) part.trail = curve(color=part.color, radius = part.radius / 8, pos=part.pos ) #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.parts : part.label.visible = True part.label.pos = part.pos wl.display_particle_labels_bool = True num_loops_to_display_part_index = 512 wl.hide_number_label_count = wl.main_loop_count + num_loops_to_display_part_index def hide_particle_labels() : for part in wl.parts : 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.parts : 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_EM_wave_height_of_particle(part): # Define where in the sine wave cycle we are in. cycle = ( wl.total_time + part.initial_wave_height_in_terms_of_dt ) * part.frequency * 2 * pi # The height of the sine wave defines are current charge part.charge_cur = ( part.charge_avg * sin(cycle) ) + part.charge_avg # Theory is that it oscillates around the average # The height of the wave defines the magnitude of our spin or intrinsic B field part.B_field_mag = ( part.avg_magnetic_moment * cos(cycle) ) + part.avg_magnetic_moment part.qv = part.charge_cur * part.velocity # http://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law#Point_charge_at_constant_velocity part.B_from_qv = K_m * part.qv # http://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law#Point_charge_at_constant_velocity 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.p_type, #print "b", part.B_field_mag, " uv", part.mag_moment_unit_vec, print " qv ", part.B_from_qv, print " spin ", part.intrinsic_magnetic_moment print " intrinsic_B_field_plus_qv", part.intrinsic_B_field_plus_qv, part.charge_ratio = part.charge_cur / part.charge_avg if (part.is_proton and wl.proton_flash) or part.is_electron : part.opacity = min_part_opacity + (part.charge_ratio*len_part_opacity/2) if ( debug_sum_forces_on_one_particle >= 20 or ( debug_sum_forces_on_one_particle >= 15 and wl.main_loop_count%256==0 )): if debug_sum_forces_on_one_particle > 25: print "charge", percent_charge_ratio = int(part.charge_ratio*100) #if debug_sum_forces_on_one_particle >= 25 or part.is_electron : # print "%%%3d" % percent_charge_ratio, if debug_verbosity > 25 : percent_diff = int((part.B_field_mag / part.avg_magnetic_moment)*100) print "spin %%%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 # end update_EM_wave_height_of_particle() ################################################################################# # 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) : if v1.mag == 0 or v2.mag == 0 : return 0 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 arrow_length_set(arrow_obj,newtons) : compared_to_max = newtons / wl.max_force # If it is 16 newtons or greater force than arrow is max length if compared_to_max > 1.0 : arrow_obj.length = max_ptr_len #arrow_obj.shaftwidth = 0.05 else : arrow_obj.length = max_ptr_len * compared_to_max if arrow_obj.length < min_ptr_len : arrow_obj.length = min_ptr_len 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 > 25 ): print " #%d" % part.index, for j in range(wl.num_particles): # 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.parts[j] part.dist_vec = part.pos - other_part.pos # dist is a vector mag2_dist = mag2(part.dist_vec) # mag2 = magnitude squared part.dist_apart = mag (part.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 # 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 : if j == 0 : if debug_sum_forces_on_one_particle > 25 and part.is_electron: print "D%1.2e" % part.dist_apart, part.current_dist_delta = part.dist_apart - part.old_dist_apart if debug_sum_forces_on_one_particle >= 20 and part.is_electron: print "d%6.0e" % (part.current_dist_delta), #if ( old_delta < 0 and current_delta > 0 ) : # part.velocity = part.velocity / 2 part.moving_further_apart = (part.dist_apart > part.old_dist_apart) # Don't know if these attributes are used. if debug_sum_forces_on_one_particle >= 20 and part.is_electron: print part.moving_further_apart, part.dist_apart, part.old_dist_apart, part.old_getting_closer = part.getting_closer part. getting_closer = not part.moving_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(part.dist_vec) e_force = dist_unit_vector * (Ke * part.charge_cur * other_part.charge_cur / mag2_dist) if ( debug_sum_forces_on_one_particle > 25 ): print "ele_force %3.0f " % mag(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: # http://en.wikipedia.org/wiki/Magnetic_moment#Two_kinds_of_magnetic_sources # # http://en.wikipedia.org/wiki/Magnetism#Sources_of_magnetism # 1. Magnetic field created from moving charge # 2. Magnetic field from standing EM wave. In this simulation also called intrinsic B field. # Commonly known as 'spin', such as electron spin # # 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 # Scroll down to: "SOURCES OF MAGNETIC FIELD (Chapter 28)" # http://faculty.sdmiramar.edu/gbochicchio/magnetic.ppt # http://physicsx.pr.erau.edu/Courses/CoursesSummerA2006/Lectures/chapter_28/chapter_28.pdf # 2. Magnetic field due to intinsic magnetic field ie spin. # Similar to that of a current loop. # Also known as a magnetic dipole. http://en.wikipedia.org/wiki/Spin_%28physics%29 # http://en.wikipedia.org/wiki/File:Magnetic_ring_dipole_field_lines.svg # http://en.wikipedia.org/wiki/Magnetic_moment#Magnetic_field_produced_by_a_magnetic_moment # # 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. # # http://en.wikipedia.org/wiki/Magnetic_moment#Magnetic_field_produced_by_a_magnetic_moment # On the wiki page look down towards the section that describes vector notation # U0/4pi=10e-7 r=dist_vec # u can be either: # 1. other_part.intrinsic_magnetic_moment or # 2. other_part.intrinsic_B_field_plus_qv r = part.dist_vec u = other_part.intrinsic_magnetic_moment #print "r", r, ", u", u, ", dot(u,r)", dot(u,r) b_field_from_other_part_at_this_point_in_space_spin = ( 10e-7 * (( 3 * r * dot(u,r) ) - ( u * mag2_dist ))) / pow(part.dist_apart,5) B_force_from_spin = cross(b_field_from_other_part_at_this_point_in_space_spin, part.intrinsic_B_field_plus_qv) b_field_from_other_part_at_this_point_in_space_qv = cross(other_part.B_from_qv,-dist_unit_vector/mag2_dist) B_force_from_qv = cross(b_field_from_other_part_at_this_point_in_space_qv, part.intrinsic_B_field_plus_qv) B_force = B_force_from_spin + B_force_from_qv if ( debug_sum_forces_on_one_particle >= 30 ): #if part.index >= 0 : #print "\t type", part.p_type, " dist apart", part.dist_apart print " i+qv %1.0e" % part.intrinsic_B_field_plus_qv.mag, print " other ", other_part.intrinsic_B_field_plus_qv print "\t b_field_from_other_part_at_this_point_in_space spin %1.0e %1.0e %1.0e" % ( b_field_from_other_part_at_this_point_in_space_spin.x, b_field_from_other_part_at_this_point_in_space_spin.y, b_field_from_other_part_at_this_point_in_space_spin.z ), print "\t qv %1.0e %1.0e %1.0e" % ( b_field_from_other_part_at_this_point_in_space_qv.x, b_field_from_other_part_at_this_point_in_space_qv.y, b_field_from_other_part_at_this_point_in_space_qv.z ) b_field_from_other_part_at_this_point_in_space = b_field_from_other_part_at_this_point_in_space_spin + \ b_field_from_other_part_at_this_point_in_space_qv print "\t total %1.0e %1.0e %1.0e" % ( 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 ), print "B %1.0e" % 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) ) ################################################################# # What happens when an electron and proton are on top of each other? # The force of attraction is infinite. This gives the simulation a big problem. # Causes particles to be thrown out. In other words electrons are ejected from the atom. # To fix digital energy gain and throw problem, limit the max force. # Throw happens when electron and proton are almost on top of each other. Force approaches infinity. ################################################################# if part.closest_part == other_part : max_force = wl.max_force # Chose a large force if e_force.mag > wl.max_force or B_force.mag > wl.max_force : if part.is_electron and debug_sum_forces_on_one_particle >= 20 : print "#%d e%3.0f b%3.0f" % ( part.index, e_force.mag, B_force.mag ), if e_force.mag > B_force.mag : ratio = wl.max_force / e_force.mag else : ratio = wl.max_force / B_force.mag e_force.mag = e_force.mag * ratio B_force.mag = B_force.mag * ratio #else : # if part.index != 0 and debug_sum_forces_on_one_particle > 20 : # print " \t \t", part.e_force_sum += e_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 if debug_sum_forces_on_one_particle >= 25 and part.index != 0 : print "e%7.0e%7.0e%7.0e%7.0e" % ( part.e_force_sum.mag, part.e_force_sum.x, part.e_force_sum.y, part.e_force_sum.z), print "b%7.0e%7.0e%7.0e%7.0e" % ( part.B_force_sum.mag, part.B_force_sum.x, part.B_force_sum.y, part.B_force_sum.z) 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)) 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 ) ) arrow_length_set(part. EF_pointer,part.e_force_mag) arrow_length_set(part. MF_pointer,part.B_force_mag) arrow_length_set(part.tot_pointer,part.t_force_mag) #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 > 50 : c_vel = ( part.velocity.mag / 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 > 25 ): 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 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 velocity_norm = norm(part.velocity) radian_angle = angle_between_two_norm_vectors(velocity_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.accel_perpen 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 >= 30 ): 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 # Below is the simple formula. Not used. Shown for reference only. part.new_velocity = part.velocity + (part.acceleration * wl.dt) # velocity = v0 + acceleration * wl.dt ac = part.accel_rq # Create a short cut because using frequently below, # rq = results queue. del ac[0] # Don't need that one anymore ac.append(part.acceleration) # put are latest result into the results queue # numerical integration using the Adams-Bashforth method # http://en.wikipedia.org/wiki/Linear_multistep_method#Adams.E2.80.93Bashforth_methods # http://kr.cs.ait.ac.th/~radok/math/mat7/step34.htm new_v_corrected = part.velocity + (wl.dt * ( \ ((1901/720)*ac[4]) - ((1387/360)*ac[3]) + ((109/30)*ac[2]) - ((637/360)*ac[1]) + ((251/720)*ac[0]) )) #part.new_velocity = new_v_corrected #if ( debug_sum_forces_on_one_particle > 15 and wl.main_loop_count%256 == 0 ) or \ if ( debug_sum_forces_on_one_particle >= 20 or ( debug_sum_forces_on_one_particle >= 15 and wl.main_loop_count%256==0 )): #if part.is_electron or wl.num_particles < 5 : #if part.is_electron or wl.num_particles < 8 : if part.is_electron : percent_charge_ratio = int(round(part.charge_ratio*100)) percent = int( 1000 * part.new_velocity.mag / c ) print "%d:%3d%%%3d" % ( part.index, percent_charge_ratio, percent ), if part.is_electron and wl.num_particles < 12 : if wl.num_particles < 10 and debug_sum_forces_on_one_particle > 35 : percent_diff = int(round((part.charge_cur / part.charge_avg)*100)) print "%3d%%" % percent_diff, v_x = int( 1000 * part.new_velocity.x / c ) v_y = int( 1000 * part.new_velocity.y / c ) v_z = int( 1000 * part.new_velocity.z / c ) print "%4d%5d%5d" % (v_x, v_y, v_z), ############################################################################################### # Below is to account for digital energy gain problem. # http://en.wikipedia.org/wiki/Energy_drift # Digital energy gain problem occurs because our DT is not ifinitely small. # Electron gains energy from proton because dt space apart greater leaving then when coming, # because of velocity diff. # Take away energy when leaving #moving_away = part.old_pos < part.pos #if moving_away or part.velocity.mag > c/16 : part.new_velocity = part.new_velocity * ( 1.0 - (1.0/(wl.energy_damp_factor*wl.num_dt_per_proton_period)) ) # Account for digital energy increase. I.E. error in simulation 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 > 50 ): print " dist trav ", part.mag_dist_traveled #print part.new_velocity part.new_pos = part.pos + part.distance_traveled ve = part.veloc_rq # Create a short cut because using frequently below, # rq = results queue. del ve[0] # Don't need that one anymore ve.append(part.new_velocity) # put are latest result into the results queue # numerical integration using the Adams-Bashforth method # http://en.wikipedia.org/wiki/Linear_multistep_method#Adams.E2.80.93Bashforth_methods # http://kr.cs.ait.ac.th/~radok/math/mat7/step34.htm # Also example in xstar code new_pos_corrected = part.pos + (wl.dt * ( \ ((1901/720)*ve[4]) - ((1387/360)*ve[3]) + ((109/30)*ve[2]) - ((637/360)*ve[1]) + ((251/720)*ve[0]) )) #part.new_pos = new_pos_corrected if debug_sum_forces_on_one_particle > 30 : print "p%7.0e%7.0e%7.0e" % ( part.new_pos.x, part.new_pos.y, part.new_pos.z ), 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() # Below is not currently used 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 %2.0e" % 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.parts[0] # initialize wl.fastest_proton = wl.parts[0] wl.fastest_electron = wl.parts[1] #wl.part_with_highest_velocity = wl.parts[0] for part in wl.parts : sum_forces_on_one_particle(part) #update_dt() for part in wl.parts : 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 record_statistics(): # What was the furthest electron? At what time did that occur? for part in wl.parts : if part.is_electron : if part.pos.mag > wl.furthest_electron.pos.mag : #wl.furthest_electron = copy.copy(part) vector_copy(wl.furthest_electron.pos, part.pos) wl.furthest_electron.time = wl.total_time wl.furthest_electron.main_loop_count = wl.main_loop_count wl.furthest_electron.label.text = part.label.text def print_details(key): furthest_proton_distance = 0.0 tot_dist = 0.0 closest_proton_distance = 1e30 num_comparisons = 0 print " iteration %d, total time %9.2e " % ( wl.main_loop_count, wl.total_time ) for part in wl.parts : if part.is_proton : # What is the average proton spread. How does it compare to start of simulation? # What is the min and the max? tot_dist = tot_dist + part.pos.mag if part.pos.mag > furthest_proton_distance : furthest_proton_distance = part.pos.mag elif part.pos.mag < closest_proton_distance : closest_proton_distance = part.pos.mag num_comparisons += 1 if key=='d' : continue; print " %3s:" % part.label.text, print "p%8.2e %8.1e %8.1e %8.1e" % ( part.pos.mag, part.pos.x, part.pos.y, part.pos.z ), percent = int( 1000 * part.velocity.mag / c ) v_x = int( 1000 * part.velocity.x / c ) v_y = int( 1000 * part.velocity.y / c ) v_z = int( 1000 * part.velocity.z / c ) print " v%3d%5d%5d%5d" % ( percent, v_x, v_y, v_z), print " a%8.1e" % ( part.acceleration.mag ) if part.index % 4 == 0 : print #print " %3d: %s p %8.1e %8.1e %8.1e %8.1e, v %8.1e %8.1e %8.1e %8.1e, a %8.1e %8.1e %8.1e %8.1e" % \ # ( part.index, part.p_type, part.pos.mag, part.pos.x, part.pos.y, part.pos.z, # part.velocity.mag, part.velocity.x, part.velocity.y, part.velocity.z, # part.acceleration.mag, part.acceleration.x, part.acceleration.y, part.acceleration.z ) print "\t furthest electron %s, distance %10.4e, time %8.2e, iter %d" % ( wl.furthest_electron.label.text, wl.furthest_electron.pos.mag, wl.furthest_electron.time, wl.furthest_electron.main_loop_count ) avg_distance = tot_dist / num_comparisons print "\t protons: closest %10.3e, furthest %10.3e, avg %10.3e, ratio to nucleus radius %7.4f" % ( closest_proton_distance, furthest_proton_distance, avg_distance, (avg_distance/wl.nucleus_radius) ) sys.stdout.flush() def implement_particle_calcs(): global debug_sum_forces_on_one_particle highest_vel_part = wl.parts[0] for i in range(wl.num_particles): decrease_vel = 0 part = wl.parts[i] vector_copy(part.old_f_unit_v, part.force_unit_v) vector_copy(part.old_vel , part. 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.num_particles*1e-12 : # Some atoms are 260e-15 in size. http://en.wikipedia.org/wiki/Atomic_radius part.new_pos = 1e-13 * norm(part.velocity) escape = True elif part.pos.mag > (wl.nucleus_radius * 2 ): # http://en.wikipedia.org/wiki/Atomic_nucleus part.new_pos = 1e-15 * norm(part.velocity) escape = True if escape : print "\t escape " , part.num, " vel", part.new_velocity.mag part.new_velocity = vector(0,0,0) part. velocity = vector(0,0,0) #debug_sum_forces_on_one_particle += 5 vector_copy(part.velocity, part.new_velocity) if part.velocity.mag > highest_vel_part.velocity.mag: 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 else : print "Velocity matches c. :0(", part.num, part.p_type, 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 ( debug_sum_forces_on_one_particle > 30 and wl.main_loop_count%128 == 0 ): # print " p%d%8.1e%6.0e%6.0e%6.0e" % (part.index, part.pos.mag, # part.pos.x, part.pos.y, part.pos.z), # # Give the user the ability to view/focus on nucleus using 'h' for hide. if part.pos.mag > wl.zoom_radius : part. visible = not wl.hide_objects_far_from_center part.vel_pointer.visible = not wl.hide_objects_far_from_center part.EF_pointer. visible = not wl.hide_objects_far_from_center part.MF_pointer. visible = not wl.hide_objects_far_from_center part.tot_pointer.visible = not wl.hide_objects_far_from_center part.mm_pointer. visible = not wl.hide_objects_far_from_center part.trail. visible = not wl.hide_objects_far_from_center else : part. visible = True part.vel_pointer.visible = True part.EF_pointer. visible = True part.MF_pointer. visible = True part.tot_pointer.visible = True part.mm_pointer. visible = True part.trail. visible = True part.vel_pointer.pos = part.pos if part.velocity != vector(0,0,0) : #part.vel_pointer.axis = norm(part.velocity) * min_ptr_len arrow_length = part.vel_max_ptr_len * ( part.velocity.mag / c ) if arrow_length < min_ptr_len : arrow_length = min_ptr_len elif arrow_length > part.vel_max_ptr_len : arrow_length = part.vel_max_ptr_len part.vel_pointer.axis = norm(part.velocity) * arrow_length 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.is_electron : if part.pos.mag > (wl.nucleus_radius*4) : and_value = wl.num_dt_per_proton_period * 16 elif part.pos.mag > (wl.nucleus_radius*3) : and_value = wl.num_dt_per_proton_period * 8 elif part.pos.mag > (wl.nucleus_radius*2) : and_value = wl.num_dt_per_proton_period * 4 elif part.pos.mag > wl.nucleus_radius : and_value = wl.num_dt_per_proton_period * 2 else : and_value = wl.num_dt_per_proton_period else : and_value = 512 * wl.num_dt_per_proton_period if ( wl.main_loop_count % and_value ) == 0 : my_retain = 16 # arbitrary, what I think looks nice if part.is_electron : my_retain = 64 part.trail.append(pos=part.pos,retain=int(my_retain)) record_statistics() # 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 " print "\t a = toggle auto zoom" print "\t h = hide objects far from center" print "\t n = show numbers" print "\t = 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 ) : if argv[0] in atom_type_to_num_particles : wl.num_particles = atom_type_to_num_particles[argv[0]] else : wl.num_particles = int(argv[0]) wl.num_dt_per_proton_period = 0 if ( argc > 1 ) : wl.num_dt_per_proton_period = int(argv[1]) wl.energy_damp_factor = 512.0 + 256.0 + 128.0 if ( argc > 2 ) : wl.energy_damp_factor = int(argv[2]) print "\t wl.energy_damp_factor", wl.energy_damp_factor #stdscr = curses.initscr() #curses.cbreak() #stdscr.nodelay(1) scene_setup() workload_setup() particles_setup() #display_particle_labels() #if ( wl.num_particles <= 2 ) : hide_particle_labels() while true: #key = ' ' #key_pressed = False #if scene.kb.keys: # did the user hit a key? # key_pressed = True #else : # key = stdscr.getch() # if key != -1 : # key_pressed = True #if key_pressed : if scene.kb.keys: # did the user hit a key? 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 == ' ': #stop, pause or play wl.play_mode = not wl.play_mode if wl.play_mode : kb_input_label.visible = 0 elif key == '-' or key == '_' : wl.num_dt_per_proton_period *= 2 # Decrease speed of simulation. Increase accuracy set_dt() elif key == '+' or key == '=': wl.num_dt_per_proton_period /= 2 # Faster simulation, less accurate. set_dt() elif key == '1' : wl.energy_damp_factor = wl.energy_damp_factor - 64.0 print "\t wl.energy_damp_factor", wl.energy_damp_factor sys.stdout.flush() elif key == '2' : wl.energy_damp_factor = wl.energy_damp_factor + 64.0 print "\t wl.energy_damp_factor", wl.energy_damp_factor sys.stdout.flush() elif key == 'a' or key == 'A' : wl.auto_zoom = not wl.auto_zoom elif key == 'd' or key == 'D' : print_details(key) elif key == 'E' or key == 'e' : # Toggle display of electric force arrow for part in wl.parts: if isinstance(part.EF_pointer,simple_class) : part.EF_pointer = my_arrow(pos=part.pos, color=(1,0,0)) else : part.EF_pointer.visible = 0 part.EF_pointer = simple_class() elif key == 'h' or key == 'H' : wl.hide_objects_far_from_center = not wl.hide_objects_far_from_center elif key == 'm' or key == 'M' : # Toggle display of magnetic force. for part in wl.parts: if isinstance(part.MF_pointer,simple_class) : part.MF_pointer = my_arrow(pos=part.pos, color=(0,0,1) ) # Where is the magnetic force pushing us? else : part.MF_pointer.visible = 0 part.MF_pointer = simple_class() elif key == 'n' or key == 'N' : if wl.display_particle_labels_bool : hide_particle_labels() else : display_particle_labels() elif key == 'p' or key == 'P': #stop, pause or play wl.proton_flash = not wl.proton_flash if not wl.proton_flash : for part in wl.parts : part.opacity = 0.3 # Arbitrary. Just chose what I think looks good. elif key == 'q' or key == 'Q': print "Exiting" sys.stdout.flush() break elif key == 'r' or key == 'R' : #print "scene.forward", scene.forward wl.rotation = not wl.rotation if not wl.rotation : scene.forward = vector(0,0,-1) rotate_count = 0 if key == 'R' : rotate_angle = 0.0001 * wl.num_particles else: rotate_angle = 0.0005 * wl.num_particles elif key == 's' or key == 'S' : # Toggle spin display for electrons for part in wl.parts: if part.is_electron : if isinstance(part.mm_pointer,simple_class) : part.mm_pointer = my_arrow(pos=part.pos, color=(0.5,0.5,0.0)) # What direction is the particle oriented to. Where is the magnetic moment vector pointing. else : part.mm_pointer.visible = 0 part.mm_pointer = simple_class() part.mm_pointer.axis = " " elif key == 't' or key == 'T' : for part in wl.parts: if isinstance(part.tot_pointer,simple_class) : part.tot_pointer= my_arrow(pos=part.pos, color=(1,1,1)) # Where is the total force pushing us? else : part.tot_pointer.visible = 0 part.tot_pointer = simple_class() #elif key == 'v' or key == 'V' : # scene.visible = not scene.visible # print "3d display is ", ("on " if scene.visible else "off") else : kb_input_label.visible = 0 kb_input_label.text = '' # erase all the text if wl.play_mode : area_of_interest = 500000000 #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%256==0 )): sys.stdout.flush() print for i in range(wl.num_particles): part = wl.parts[i] update_EM_wave_height_of_particle(part) # Update the current magnitude of the e and m wave. particle_calcs() if debug_verbosity >= 15 and (wl.main_loop_count%256) == 0 : print "%d" % (wl.main_loop_count/256), wl.total_time += wl.dt wl.main_loop_count += 1 #if wl.main_loop_count > 25 : # sys.exit(0) #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() if wl.rotation : scene.forward = scene.forward.rotate(angle=rotate_angle, axis=(0,1,0)) #print "scene.forward", scene.forward, "rotate_count", rotate_count if scene.forward.z < -0.999 and rotate_count > 200: wl.rotation = False rotate_count = rotate_count + 1 else : rate(1) 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. """