Well, I did something productive: I wrote Python code that simulates the paths of ions in a quadrupole mass spectrometer.
I still need to write something to actually put this to use.
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 21 23:38:33 2018
@author: ikrase
Code to simulate the paths and stability of ions in quadrupole mass
spectrometers.
This is heavily based on the reference of
http://www.chem.uky.edu/research/Lynn/pdfs/JChemEd_quads.pdf
"""
import numpy as np
import scipy as sp
import scipy.integrate
from math import *
from numpy import array
def fMathieu_qp(t, u, au, qu):
''' RHS differential equation of Matthieu for quadrupoles
u = X/Y position, combined (ODE vector)
t = quasi-time -- normalized to angular frequency.
au, qu = dimensionless parameters, functions of all quadrupole properties.
(See referenced paper)
'''
return np.array([u[1], -(au + 2*qu * cos(2*t))*u[0]])
#return np.array([u[2], - (w0 * w0 + eps * math.cos(w1*t))*u[1]])
def alt_fMathieu_qp(t,u,au,qu,freq):
'''Alternative non-combined equation of Matthieu, in strict time domain.
Freq is the actual frequency in hertz.
u is [x, x', y, y']. '''
w = 2 * pi * freq
xpp = - ((w/2)**2) * (au + 2*qu*cos(w*t))*u[0]
# at some point I'll figure out which is the squshed axis and which is stretched.
ypp = ((w/2)**2) * (au + 2*qu*cos(w*t))*u[2]
return array([u[1], xpp, u[3], ypp])
def quadrupol_test(r, m, Vrf, Udc, freq, duration=False, tstep = 0.1E-7):
'''Test whether a given quadrupole system is stable or unstable.
r: quadrupole center-to-electrode radius / inscribed radius in mm.
m: ion M/Z ratio in amu to electrons
Vrf: RF zero-to-peak voltage in volts
Udc: DC spread voltage
freq: Actual frequency (not angular) of the quadruopole in hertz
'''
rspan = r * 0.125 # Radius of input particle beam
aqbot = ((2 * pi * freq * r)**2 * 1.036427E-8 * m ) # bottom of A and Q
a = 4 * Udc / aqbot
q = 2 * Vrf / aqbot
tb = min(0.002, 2, sqrt(a),abs(2-sqrt(a)))
if tb < 1E-3: tb = 5E-2 #divide by zero and over-value guard.
tlong = 2 * pi / tb # Estimate of "natural" frequency
# (if there is one.)
if duration == False:
dur = 10 * tlong
else:
dur = duration
#oderes = sp.integrate.solve_ivp(lambda tt, yy: fMathieu_qp(tt, yy, a,q),
# (0,dur), np.array([rspan, 0]), max_step = 0.05)
oderes = sp.integrate.solve_ivp(lambda tt, yy: alt_fMathieu_qp(tt, yy, a,q, freq),
(0,dur), np.array([rspan, 0, rspan, 0]),max_step = tstep)
if not oderes.success:
raise RuntimeError("Failed to integrate differential equation")
return oderes