bmarechal / sicav_sim

Download as
Browse Code ยป

Commit 01f78001077ca38a49c870bf40a7920037453f44

Authored by bmarechal
1 parent e5671c2970
Exists in master

rename ram to phase_modulation (PM) sim

Showing 4 changed files with 92 additions and 0 deletions Side-by-side Diff

  1 +from sympy import *
  2 +from sympy.simplify.fu import *
  3 +
  4 +init_printing()
  5 +
  6 +#constants
  7 +E0, W, t, m, b= symbols('E0, Omega, t, m, beta', real=True, imaginary=False)
  8 +
  9 +#AM+PM
  10 +E_ampm = E0*(1+m*cos(W*t))*exp(I*(W*t+b*sin(W*t)))
  11 +I_ampm = abs(E_ampm)**2
  12 +I_ampm = expand(TR8(expand(expand_complex(I_ampm))))
  13 +
  14 +#AM
  15 +E_am = E0*(1+m*cos(W*t))*exp(I*(W*t))
  16 +I_am = abs(E_am)**2
  17 +I_am = expand(TR8(expand(expand_complex(I_am))))
  18 +
  19 +#PM
  20 +E_pm = E0*exp(I*(W*t+b*sin(W*t)))
  21 +I_pm = abs(E_pm)**2
  22 +I_pm = expand(TR8(expand(expand_complex(I_pm))))
  1 +% Phasor demo
  2 +%
  3 +% Tobin Fricke <tobin.fricke@ligo.org> 2011-06-08
  4 +
  5 +% modulation depth
  6 +%m = pi/4;
  7 +m = 1;
  8 +
  9 +% vector of (normalized) frequencies to include
  10 +%freqs = [0 1 -1 2 -2 3 -3 4 -4];
  11 +freqs = [0 1 -1 2 -2];
  12 +
  13 +% initial amplitudes of the phasors
  14 +amp = besselj(freqs, m);
  15 +
  16 +set(0, 'defaultlinelinewidth', 5);
  17 +
  18 +% Set up the axes
  19 +cla reset;
  20 +axis equal;
  21 +xlim([-1 1.5]);
  22 +ylim([-1 1]);
  23 +
  24 +% Draw the unit circle
  25 +rectangle('Position', [-1 -1 2 2], 'Curvature', [1 1])
  26 +
  27 +% make some lines
  28 +L = [];
  29 +L(1) = line([0], [0]);
  30 +L(2) = line([0], [0], 'color', 'red');
  31 +
  32 +save_to_disk = 0;
  33 +
  34 +% how many frames to render?
  35 +if save_to_disk,
  36 + n_frames = 200;
  37 +else
  38 + n_frames = inf;
  39 +end
  40 +
  41 +% how much the 1-Omega should rotate in one frame of the animation?
  42 +if isinf(n_frames)
  43 + dphi = m/20;
  44 +else
  45 + dphi = (2*pi)/(n_frames + 1);
  46 +end
  47 +
  48 +frameno = 0;
  49 +ret_amp = [];
  50 +while (frameno < n_frames),
  51 + frameno = frameno + 1;
  52 +
  53 + % Add up the phasors cumulatively
  54 + seg = cumsum([0 amp]);
  55 +
  56 + % Draw both the collection of phasors and the resultant
  57 + set(L(1), 'xdata', real(seg), 'ydata', imag(seg));
  58 + set(L(2), 'xdata', real(seg([1,end])), 'ydata', imag(seg([1,end])));
  59 +
  60 + ret_amp = [ret_amp, abs(seg([end]))];
  61 +
  62 + % Update the display
  63 + drawnow;
  64 +
  65 + if (save_to_disk),
  66 + print('-dpng', '-r72', sprintf('pm%03d.png', frameno));
  67 + else
  68 + pause(0.02);
  69 + end
  70 +
  71 + % Advance the phasors for the next frame
  72 + amp = amp .* exp(1i * dphi * freqs);
  73 +end
  1 +clc
  2 +clear all
  3 +close all
  4 +
  5 +N = 10;
  6 +
  7 +set(0, 'defaultlinelinewidth', 5);
  8 +
  9 +phi = [0:0.01:1]*2*pi;
  10 +m = 1;
  11 +
  12 +for n = 2:N
  13 + freqs = -n:n;
  14 + amp = besselj(freqs, m);
  15 + plot(phi, abs(sum(amp'.*exp(1i.*phi.*freqs'))))
  16 + hold on
  17 +end
  18 +
  19 +legend(strsplit(num2str(2:N)))
1   -from sympy import *
2   -from sympy.simplify.fu import *
3   -
4   -init_printing()
5   -
6   -#constants
7   -E0, W, t, m, b= symbols('E0, Omega, t, m, beta', real=True, imaginary=False)
8   -
9   -#AM+PM
10   -E_ampm = E0*(1+m*cos(W*t))*exp(I*(W*t+b*sin(W*t)))
11   -I_ampm = abs(E_ampm)**2
12   -I_ampm = expand(TR8(expand(expand_complex(I_ampm))))
13   -
14   -#AM
15   -E_am = E0*(1+m*cos(W*t))*exp(I*(W*t))
16   -I_am = abs(E_am)**2
17   -I_am = expand(TR8(expand(expand_complex(I_am))))
18   -
19   -#PM
20   -E_pm = E0*exp(I*(W*t+b*sin(W*t)))
21   -I_pm = abs(E_pm)**2
22   -I_pm = expand(TR8(expand(expand_complex(I_pm))))