diff --git a/ram/ram.py b/PM/PM.py
similarity index 100%
rename from ram/ram.py
rename to PM/PM.py
diff --git a/PM/phasordemo.m b/PM/phasordemo.m
new file mode 100644
index 0000000..0a5f45b
--- /dev/null
+++ b/PM/phasordemo.m
@@ -0,0 +1,73 @@
+% Phasor demo
+%
+% Tobin Fricke <tobin.fricke@ligo.org> 2011-06-08
+
+% modulation depth
+%m = pi/4;
+m = 1;
+
+% vector of (normalized) frequencies to include
+%freqs = [0 1 -1 2 -2 3 -3 4 -4];
+freqs = [0 1 -1 2 -2];
+
+% initial amplitudes of the phasors
+amp = besselj(freqs, m);
+
+set(0, 'defaultlinelinewidth', 5);
+
+% Set up the axes
+cla reset;
+axis equal;
+xlim([-1 1.5]);
+ylim([-1 1]);
+
+% Draw the unit circle
+rectangle('Position', [-1 -1 2 2], 'Curvature', [1 1])
+
+% make some lines
+L = [];
+L(1) = line([0], [0]);
+L(2) = line([0], [0], 'color', 'red');
+
+save_to_disk = 0;
+
+% how many frames to render?
+if save_to_disk,
+    n_frames = 200;
+else
+    n_frames = inf;
+end
+
+% how much the 1-Omega should rotate in one frame of the animation?
+if isinf(n_frames)
+    dphi = m/20;
+else
+    dphi = (2*pi)/(n_frames + 1);
+end
+
+frameno = 0;
+ret_amp = [];
+while (frameno < n_frames),
+    frameno = frameno + 1;
+
+    % Add up the phasors cumulatively
+    seg = cumsum([0 amp]);
+
+    % Draw both the collection of phasors and the resultant
+    set(L(1), 'xdata', real(seg), 'ydata', imag(seg));
+    set(L(2), 'xdata', real(seg([1,end])), 'ydata', imag(seg([1,end])));
+
+	ret_amp = [ret_amp, abs(seg([end]))];
+
+    % Update the display
+    drawnow;
+
+    if (save_to_disk),
+        print('-dpng', '-r72', sprintf('pm%03d.png', frameno));
+    else
+        pause(0.02);
+    end
+
+    % Advance the phasors for the next frame
+    amp = amp .* exp(1i * dphi * freqs);
+end
diff --git a/PM/pm_vs_besselj.m b/PM/pm_vs_besselj.m
new file mode 100644
index 0000000..eb4718d
--- /dev/null
+++ b/PM/pm_vs_besselj.m
@@ -0,0 +1,19 @@
+clc
+clear all
+close all
+
+N = 10;
+
+set(0, 'defaultlinelinewidth', 5);
+
+phi = [0:0.01:1]*2*pi;
+m = 1;
+
+for n = 2:N
+	freqs = -n:n;
+	amp = besselj(freqs, m);
+	plot(phi, abs(sum(amp'.*exp(1i.*phi.*freqs'))))
+	hold on
+end
+
+legend(strsplit(num2str(2:N)))