GSoC Week9 : Beginning of electrostatics

Week 8 of GSoC is now over, and just 5 more weeks remain till the hard pencils-down date.
However, my work on the project will most likely exceed that time, since I will be having my first tests around September 17. I am sure I will be able to get all the code in by that time, but to make the entire work shippable, things like Sphinx documentation, doctests etc will also have to be finished, which may take a little more time beyond the actual deadline.

Anyways, coming to this week's work, I wrote some code for electrostatics. I added some functions and the ParticleCharge class (inheriting from Particle of mechanics module).
However, the E-M module's handling of charged particles required me to modify my earlier implementation of Particle- most importantly, adding a 'set_motion' method. Like I mentioned earlier, Particle earlier had a frame attached to it. Gilbert and I discussed at length on this, and we agreed that making the user initialize a new frame everytime, just to attach it to a Particle, was quite cumbersome.
Hence, I modified the class in the mechanics core- It still has a frame as an attribute, but its now 'hidden' from the user (nothing in Python in hidden from the user, but what I mean here is that the user can do everything with a Particle without caring about the frame). The set_motion method I mentioned earlier just initializes a new frame according to the motion parameters set by the user (similar to MovingRefFrame).

The main reason I want a frame attached to a Particle is the ease of inner working. Since a Particle will not have an 'angular acceleration', those parts of the MovingRefFrame code will never really get executed in Particle. However, Gilbert has a point in pointing out that its unnecessarily adding a complex attribute to the class. We havent finalized anything yet, but re-writing the motion-setting methods for Particle would be tedious, and more importantly, will lead to quite a lot of code duplication. Lets see what we end up doing.

Anyways, I also finished writing tests for the code I wrote this week, and I am waiting for a review from Gilbert. Jason and he also encouraged me to add example-tests to ensure everything works smoothly.
Here are the example tests I wrote (they are really basic, but its better that way, since strenuous tests only make things worse)-

 def test_example_1():  
   """  
   This is a sample example meant to test the basic functionality of  
   electrostatics.py  
   
   Consider a frame R, and a particle charge P(mass m, charge q) situated  
   at its origin initially.  
   Now, at time t=0, a field R.x+R.y+R.z is 'switched on' and kept that way  
   for time 't0'.  
   The user needs the x-coordinate of the particle charge as a function of  
   time, after the field is switched off at time=t0.  
   """  
     
   #Basic initialization  
   m, q = symbols('m q')  
   P = ParticleCharge('P', m, q)  
   P.set_motion(R, pos_vector = 0)  
   field = R.x + R.y + R.z  
   time = dynamicsymbols._t  
   t0 = Symbol('t0')  
   #The acceleration is equal to the electrostatic force experience by the  
   #particle charge, divided by its mass  
   acceleration = P.electrostatic_force(field)/P.mass  
   #Use get_motion_acc from the mechanics core to find velocity and position  
   #parameters using acceleration function and boundary conditions  
   translation = get_motion_acc(acceleration, \  
                  0, P.pos_vector_wrt(R), frame = R)  
   #Calculate the motion parameters of the particle charge at time t0  
   acceleration = translation[0]  
   velocity_at_t0 = translation[1].subs({time:t0})  
   position_at_t0 = translation[2].subs({time:t0})  
   #Set motion of P accordingly  
   P.set_motion(R, trans_acc = 0, trans_vel_b = velocity_at_t0,  
          pos_vector_b = position_at_t0)  
   #assert that the dot product of P's pos_vector wrt R, with R.x  
   #matches the actual result  
   assert P.pos_vector_wrt(R).dot(R.x) == \  
       q*t*t0/m + q*t0**2/(2*m)  
   
   
 def test_example_2():  
   """  
   This is a sample example meant to test the basic functionality of  
   electrostatics.py  
   
   Consider a set of 3 particle charges of equal mass and charge placed  
   symmetrically on a circle of unit radius around the origin, with one  
   point lying at position vector of R.y  
   We calculate the energy required to assemble this system of charges from  
   infinity(with no kinetic energy in any particle) and find out values of  
   the electrostatic fields and potentials at the origin  
   """  
   
   #Basic initialization  
   from sympy import sin, cos, pi  
   m, q = symbols('m q')  
   P = ParticleCharge('P', m, q)  
   Q = ParticleCharge('Q', m, q)  
   S = ParticleCharge('S', m, q)  
   #Place particle charges as per required configuration  
   P.set_motion(R, pos_vector = R.y)  
   Q.set_motion(R, pos_vector = -cos(pi/6) * R.x + -sin(pi/6) * R.y)  
   S.set_motion(R, pos_vector = cos(pi/6) * R.x + -sin(pi/6) * R.y)  
   #Check outputs of required values  
   assert charge_assembly_energy(P, Q, S) == 3*3**(-0.5)*k*q**2  
   assert S.electrostatic_field(R, 0) + \  
       Q.electrostatic_field(R, 0) + \  
       P.electrostatic_field(R, 0) == 0  
   assert P.electrostatic_potential(R, 0) == \  
       Q.electrostatic_potential(R, 0) == \  
       S.electrostatic_potential(R, 0)  

If you think the API is a little complex, I think so too, but those modifications are for later. It would be good, according to me, to have something like an 'environment' infrastructure. Basically, the user will keep adding components (fields, particles etc etc) to it, and the constituents will keep getting modified accordingly. However, it will be good to have the module at a certain level of working before thinking on those lines.

That's all there is to report for this week. The coming week, I will be reading up on magnetostatics from Griffith's and start coding for the same. However, I am afraid things are only going to get more and more complex hereon. Lets hope for the best.

Have a great week :-)