% compile and create simulation class
clear classes
MIST('fortymeters_pop_qpd.mist');
s = FortyMetersPOP_QPD(4);

% set angular motion of ITMs, ETMs and PRM
s.ETMX.setMotionShape('pitch');
s.ETMY.setMotionShape('pitch');
s.ITMX.setMotionShape('pitch');
s.ITMY.setMotionShape('pitch');
s.PRM.setMotionShape('pitch');

% list of mirrors
mirrors = {'ETMX', 'ETMY', 'ITMX', 'ITMY', 'PRM'};

% Gouy phases
npt = 100;
phi = linspace(0,180,npt);

%% Loop over all listed mirrors
results = zeros(numel(mirrors), npt);
results22 = zeros(numel(mirrors), npt);
for i=1:numel(mirrors)
    fprintf(1, '%s \n', mirrors{i});
    % reset all amplitudes to zero
    for j=1:numel(mirrors)
        s.(['z_mir_', mirrors{j}]) = 0;
    end
    % set the mirror amplitude
    s.(['z_mir_', mirrors{i}]) = 1;

    % scan Gouy phase of space before the QPD
    for k=1:npt
        s.sPOPphase.gouy = phi(k)/180*pi;
        results(i,k) = s.POP_QPD_y_TF_mir();
        results22(i,k) = s.POP_QPD_22_y_TF_mir();
    end
end

figure()
subplot(211)
plot(phi, abs(results), 'LineWidth', 2)
xlabel('Gouy phase from PR2 [deg.]')
ylabel('Signal [W/rad]')
xlim([min(phi), max(phi)])
title('POP QPD')
legend(mirrors)

subplot(212)
plot(phi, abs(results22), 'LineWidth', 2)
xlabel('Gouy phase from PR2 [deg.]')
ylabel('Signal [W/rad]')
xlim([min(phi), max(phi)])
title('POP QPD (22 MHz)')
legend(mirrors)