Check math in code for em-force on moving particle due to finite straight wire

[happyjack27]

please check math in code for em-force on moving particle due to finite straight wire

could someone pretty please make sure I'm doing the calculations correctly?

this is for a computer simulation of charged particles in an electromagnetic field.

it's to calculate the electro-static and magneto-static force on a moving charged particle due to a voltage and current in a finite straight wire segment.

the point is at position "p" and velocity "v", and the wire is from point "a" to point "b".
the wire has current density a.w*current_density and voltage density b.w*voltage_density.

i've annotated the code for you.

and there is sugar on top of this please.

/*========= declaration of basic data types (T = scalar, T4 = vector) =========*/

typedef float T;
struct T3 { T x,y,z; };
struct T4 { T x,y,z,w; };/*========= basic vector geometry functions =========*/

#define translate(a,b) { a.x-b.x, a.y-b.y, a.z-b.z } //translate vector a by b
#define dot(a,b) (a.x*b.x + a.y*b.y + a.z*b.z) //dot product
#define SCALE(a) rsqrt(dot(a,a))  //reciprocal of the length of vector a.

//cross product stuff
#define crossx(a,b) (a.y*b.z - a.z*b.y)
#define crossy(a,b) (a.z*b.x - a.x*b.z)
#define crossz(a,b) (a.x*b.y - a.y*b.x)
#define cross(a,b) { crossx(a,b), crossy(a,b), crossz(a,b), a.w*b.w }

#define norm(a) { a.x*a.w, a.y*a.w, a.z*a.w, 1.0f }  //norm of vector a

#define project(a,b) (dot(a,b)*b.w) //length of projection of a onto b
#define cosangle(a,b) (project(a,b)*a.w) //cosine of the angle between vector a and b./*========= given a point charge p, traveling at velocity v, and a finite straight wire from point a to b, find the acceleration of p due to the magnetic and electric field induced from the current and voltage in the wire onto the point p: =========*/

T3 apply_static_fields_of_seg_on_point(T3 accel, T4 p, T4 v, T4 a, T4 b, T current_density, T voltage_density) {

   /*========= convert point pairs to vectors =========*/

   T4 lvector = translate(b,a);
   lvector.w = SCALE(lvector);

   T4 v1 = translate(p,a);
   v1.w = SCALE(v1);

   T4 v2 = translate(p,b);
   v2.w = SCALE(v2);   /*========= find the line perdenticular to the wire, through the point p ("r") ========*/
   // used the projection of a->p onto a->b to find the point where the perpendicular line intersects a->b (call it "c").  the perpendicular line is then c->p.

   //proj = length of projection of v1 onto lvector, times the reciprocal magnitude of lvector.
   T proj = project(v1,lvector);
   T adjproj = proj*lvector.w;

   //r = perpendicular vector from line to point.
   T4 r;
   r.x = v1.x-lvector.x*adjproj; //=p.x-(a.x+lvector.x*adjproj);
   r.y = v1.y-lvector.y*adjproj;
   r.z = v1.z-lvector.z*adjproj;
   T h2 = dot(r,r);
   r.w = rsqrt(h2);   /*========= now calculate the magnetic field ("bfield") =========*/
   // formulas taken from http://cnx.org/content/m31103/latest/  (thou that site doesn't calculate direction of field, i added that per the basic definition)

   T4 cr = cross(lvector,r);
   T4 bfield = norm(cr); //this is the direction of the magnetic field.  <-correct??
   
   T bstrength = cosangle(lvector,v1);
   lvector.w = -lvector.w;   //this effectively switches the direction of the vector.
   bstrength += cosangle(lvector,v2);
   lvector.w = -lvector.w;   //dont forget to switch it back

   bfield.w = bstrength*r.w*current_density*a.w*m0/(4*pi);  //a.w is current density. m0 is mag constant.  this is the relative strength of the magnetic field.  <-correct??   /*========= now calculate the electric field ("efield") =========*/
   // formulas taken from http://planetphysics.org/encyclopedia/ElectricFieldOfALineOfCharge.html

   T nproj = 1.0/lvector.w - proj;
   T comp1 = rqsrt(nproj*nproj+h2);
   T comp2 = rqsrt(proj*proj+h2);

   T xcomp = comp1-comp2;  //component in direction of a->b (lvector)
   T ycomp = (comp1*nproj-comp2*proj)*r.w; //component in direction of r

   T eforce_multiplier = voltage_density * b.w / (4*pi*e0); //b.w is charge density e0 is the electric constant
   xcomp *= lvector.w * eforce_multiplier; 
   ycomp *= r.w * eforce_multiplier;

   T4 efield;
   efield.x = lvector.x*xcomp + r.x*ycomp; //strength is included in the vector components now.
   efield.y = lvector.y*xcomp + r.y*ycomp;
   efield.z = lvector.z*xcomp + r.z*ycomp;   /*========= finally, add their forces to the acceleration vector ("accel") =========*/
   // basic "lorentz force" calculation.  (the charge-to-mass ratio of particle p will be factored in elsewhere.)

   accel.x += (efield.x + crossx(v,bfield)*bfield.w);
   accel.y += (efield.y + crossy(v,bfield)*bfield.w);
   accel.z += (efield.z + crossz(v,bfield)*bfield.w);   /*========= and we're done. =========*/

   return accel;
} /*========= loop thru all the wire segments and apply their forces to the acceleration vector for the point charge =========*/

T3 apply_static_field(T3 accel, T4 pos, T4 vel, T4* as, T4* bs, int num, T current, T voltage) {
	for( int i = 0; i < num; i++) {
		accel = apply_static_fields_of_seg_on_point(accel,pos,vel,as[i],bs[i],current,voltage);
	}
	return accel;
}

thank you!

 

[AndyW]

Here are some potential issues I noticed in the code:

1. In the function "apply_static_fields_of_seg_on_point", there seems to be a typo in the calculation of "bfield.w". It should be multiplied by "a.w" instead of "r.w".

2. In the calculation of "eforce_multiplier", there is a missing parenthesis. It should be "voltage_density * b.w / (4*pi*e0)".

3. The calculation of "xcomp" and "ycomp" in the electric field calculation seems to be incorrect. It should be "xcomp = comp1 * lvector.w - comp2 * proj" and "ycomp = (comp1 * nproj - comp2 * proj) * r.w".

4. The calculation of "efield" in the electric field calculation seems to be incorrect. It should be "efield.x = lvector.x * xcomp + r.x * ycomp" and so on.

Overall, it is recommended to double-check the code and make sure all the calculations are correct before using it for simulations.