Moving over the paultrap repository to be a gist since it was only one file anyway.

paultrap.py

"""paultrap.py - Provides a simple object for calculating stability
boundaries for a linear Paul trap.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""

import numpy
import scipy.constants

q_e = scipy.constants.e
amu = scipy.constants.u
pi = numpy.pi

A = [-6.43e-2, 0.5256, -3.23e-3, 0]
B = [-1.173, 1.0657]
q1, q2 = numpy.arange(0, 0.706, 0.005), numpy.arange(0.706, 0.915, 0.005)
q = numpy.append(q1, q2)

class PaulTrap:
    def __init__(self, radius, frequency):
        """Initialize PaulTrap object with radius r0 and frequency f
        in SI units.

        """
        self.r0 = radius
        self.Omega = 2*pi*frequency

    def stability_qa(self):
        """Returns the q, a boundary of stability for the Mathieu
        equations.

        """
        return q, numpy.append(numpy.polyval(A, q1), numpy.polyval(B, q2))
        
    def stability_VU(self, m, Z=1.):
        """Returns the Mathieu stability boundary in voltages for mass
        given in amu and charge in units of e.

        """
        q, a = self.stability_qa()
        V = q*m*amu*self.r0**2*self.Omega**2/(2.*Z*q_e)
        U = a*m*amu*self.r0**2*self.Omega**2/(4.*Z*q_e)
        return V, U

    def qa_to_voltage(self, q, a, m, Z=1.):
        """Return the voltages for given stability parameters. Mass
        given in amu.

        """
        return [q*m*amu*self.r0**2*self.Omega**2/(2.*Z*q_e),
                a*m*amu*self.r0**2*self.Omega**2/(4.*Z*q_e)]

    def voltage_to_qa(self, V, U, m, Z=1.):
        """Returns the stability parameters for given voltages. Mass
        given in amu.

        """
        return [2*Z*q_e*V/(m*amu*self.r0**2*self.Omega**2),
                4*Z*q_e*U/(m*amu*self.r0**2*self.Omega**2)]

if __name__ == "__main__":
    trap = PaulTrap(3.18e-3, 2.7e6)
    print trap.stability_qa()
    print trap.stability_VU(138.)
    print trap.qa_to_voltage(0.6, 0., 138.)
    print trap.voltage_to_qa(230., 0., 138.)