%% Simulation of quantum noise with MIST - part 4
%
% This example is the fourth part of a series demonstrating how to simulate
% quantum noise with MIST. In this example we simulate a system that
% exhibits ponderomotive squeezing.

clear all
close all
addpath(genpath('/Users/keiko/!Work/MIST'))

MIST('KAGRA_homodyne.mist');
s = PonderomotiveSqueezer(0);
%%
% All mirrors are suspended to pendulums with a resonant frequency of 1 Hz.
% The two ITMs have a mass of 250g, while the ETM has a mass of only 1g:
massITM = 30;%0.250;
massETM = 30;%0.001;
w = 2 * pi * 0.8;
dampRes = [0.001 + 1i, 0.001 - 1i];

s.ITMX.setMechanics(zpk([], -w*dampRes, 1/massITM));
s.ITMY.setMechanics(zpk([], -w*dampRes, 1/massITM));
s.ETMX.setMechanics(zpk([], -w*dampRes, 1/massETM));
s.ETMY.setMechanics(zpk([], -w*dampRes, 1/massETM));


%% Transfer functions
%
% As a first step we simply simulate the opto-mechanical transfer
% functions from common and differential dephasing in the arms to the
% anti symmetric port detector. Let's start with the common mode
fr = logspace(-1,4,100);

s.z_drv_sXCAV = 1;
s.z_drv_sYCAV = 1;
Tcomm = s.TransferFunction('drv', 'AP', fr);

%%
% and then continue with the differential mode
s.z_drv_sXCAV = 0.5;
s.z_drv_sYCAV = -0.5;
Tdiff = s.TransferFunction('drv', 'AP', fr);

%%
% and we plot both transfer functions 
figure()
loglog(fr, abs(Tcomm), 'r-', fr, abs(Tdiff), 'b-', 'LineWidth', 2)
legend('Comm. mode', 'Diff. mode', 'location', 'best')
xlabel('Frequency [Hz]')
ylabel('Opto-mechanical response [W/m]')

%% Quantum noise
%
% In the configuration file we are injecting quantum noise only from the
% anti symmetric port. This is the most relevant contribution. We then
% measure the quantum noise on the anti symmetric port detector. 

noise = s.QuantumNoise('AP', fr);
shot = sqrt(2*s.h_ * s.c_/s.lambda * s.AP);

figure(10)
semilogx(fr, 20*log10(noise/shot), 'LineWidth', 2)
xlabel('Frequency [Hz]')
ylabel('Squeezing [db]')

figure(11)
loglog(fr, noise./abs(Tdiff).'/3000)
xlabel('Frequency [Hz]')
ylabel('Sensitivity [Hz^{-1/2}]')

%% Quantum noise
%
% We can change the homodyne readout phase by tuning the propagation phase
% in one of spaces that connect the input pick-off with the detectors. The
% optimal phase, i.e. the one that give the lowest quantum noise, is 0.83
% radians. This value has been found by scanning the phase and simulating
% the quantum noise at a fixed frequency.
s.sHomodyne.phi = 0.83;

%% We compute the quantum noise for both photodiodes
noise = s.QuantumNoise({'PD1', 'PD2'}, fr);

%%
% then we sum the two noises in quadrature and normalize with respect to
% standard quantum vacuum (shot-noise):
noise_pd = sqrt(noise(1,:).^2 + noise(2,:).^2);
shot_pd = sqrt(2*s.h_ * s.c_/s.lambda * (s.PD1 + s.PD2));
%%
% Here is the result:
figure(12)
semilogx(fr, 20*log10(noise_pd/shot_pd), 'LineWidth', 2)
xlabel('Frequency [Hz]')
ylabel('Squeezing [db]')