C++ Model of a charged particle in a magnetic field

[Dissident]

Homework Statement

Im taking my first physics class and We have to do a research project. The project can be anything we want. I decided to make a c++ program to model a charged particle in a magnetic field. When I first researched the problem I't seemed straight forward enough. The force on a particle in a magnetic field is givin by F = q(v X B). I thought I could solve for acceleration and then integrate via RK4 to produce a velocity for time t+delta t. Then using the previous velocity and the acceleration form the Lorentz force I would get the next position of the particle using x = v*t + 1/2 * a * t^2. I thought I could do this in a loop to generate a data file witch I could then graph using gnuplot.

The problem is that the force caused my the B field will not change the speed of the particle only the direction. When I integrate numerical the speed changes. If I try to normalize the new velocity and scale by the previous velocity the radius of the path does not line up with a the radius I calculate by hand. I have tried RK4 and the trapezoid method. I am now thinking about modeling the particle and a constant circular motion around a drift point, but I cannot for the life of me figure out how to model constant circular motion in three dimentions.

I have hit a wall and cannot get this to work.

Homework Equations

F = ma
F = qv X B
F = mv^2 / r
r = v^2 / a
r = mv / |q|B

The Attempt at a Solution

I don't think I should post 1000+ lines of code. Here is my latest attempt at integration. I am using the mpf_t form the GNU MP library which I have wrapped in a C++ class. hpv3d_t is a vector class that uses the wrapped mpf_t.
hpv3d_t k1, k2, k3, k4;
k1 = ((eField + velocity % bField) * chargeDivMass) * deltaTime;
k2 = ((eField + (velocity + k1 * 0.5) % bField) * chargeDivMass) * deltaTime;
k3 = ((eField + (velocity + k2 * 0.5) % bField) * chargeDivMass) * deltaTime;
k4 = ((eField + (velocity + k3) % bField) * chargeDivMass) * deltaTime;
velocity = velocity + (k1 + k2 * 2.0 + k3 * 2.0 + k4) * (deltaTime / 6);

the eField is <0,0,0> for now

 

[jbsteele]

and the bField is <0, 0, B>. I have tried using a constant velocity but that doesn't seem to work either.