%% KAGRA QN calculation
% This file is depeloped based on the MIST example file, QN example part 2.
%   Keiko, June 2014
%
clear all
close all
addpath(genpath('/Users/keiko/!Work/MIST'))
%% Simulation of quantum noise with MIST - part 2
%
% This example is the second part of a series demonstrating how to simulate
% quantum noise with MIST. In this example we compute the quantum noise
% limited sensitivity of advanced LIGO and then simualte the effect of
% squeezing and anti-squeezing.

%% The configuration file

MIST('KAGRA_bb.mist');
s = DualRecycledBroadBand(0);
%%
% Since radiation pressure is taken into account for the four test masses,
% we have to define their mechanical response (simple pendulum resonance at
% 0.8 Hz with an moderately high quality factor)
% taken from Somiya Optickle model

w = 2 * pi * 0.8;
dampRes = [0.001 + 1i, 0.001 - 1i];
mass = 30;
s.ITMX.setMechanics(zpk([], -w*dampRes, 1/mass));
s.ETMX.setMechanics(zpk([], -w*dampRes, 1/mass));
s.ITMY.setMechanics(zpk([], -w*dampRes, 1/mass));
s.ETMY.setMechanics(zpk([], -w*dampRes, 1/mass));

%%
% As explained in the previous example, MIST does not add automatically 
% quantum vacuum sources (yet). So we have to add one source manually. We
% already know that the dominant contribution to quantum noise comes from
% the vacuum entering from the anti symmetric port. Therefore the
% configuration file contains the following command:
%
%   vacuum qv nOMC2 into=OMC
%
% where we named the source |qv| and specified that vacuum is added in
% trasmission of the OMC and it actually propagate into the OMC.

%% The DARM transfer function
%
% The first step is to compute the calibration of the anti symmetric port
% photodetetor, or in other words the transfer function from DARM to the AS
% signal (placed after an ideal OMC). In the configuration file a driver
% called |mir| has been added. Using this driver we can compute the
% response to a GW signal, simulated by modulating the two arm propagation
% in a differential way:
s.z_drv_sXCAV = 0.5;
s.z_drv_sYCAV = -0.5;
%%
% Having defined the motion we are interested in, we can use the
% |TransferFunction| method to compute the calibration we need.
fr = logspace(-1,4,100);
T = s.TransferFunction('drv', 'AP', fr);
%%
% Remember that when a space is driven, the amplitude motion is measured in
% meters, therefore this transfer function is in units of W/m.

%% Quantum noise
%
% The computation of quantum noise is quite straightforward, using the
% |QuantumNoise| method:
noise = s.QuantumNoise('AP', fr);
%%
% And finally here is the results
figure()
loglog(fr, noise./abs(T)' / 3000, 'r', 'LineWidth', 2)
xlim([0.1, 10000])
xlabel('Frequency [Hz]')
ylabel('Sensitivity [Hz^{-1/2}]')
title('Dual recycled -  zero detuning')


%% Squeezing and anti-squeezing
%
% Including frequency independent squeezing is simple, and we can do it in
% the same way we used for the previous example:
s.qv.db = 6;
noise2 = s.QuantumNoise('AP', fr);
%%
% Similarly, we can introduce some anti squeezing:
s.qv.db = -6;
noise3 = s.QuantumNoise('AP', fr);
%%
% Here is a comparison of the various results
figure()
loglog(fr, noise./abs(T)' / 3000, fr, noise2./abs(T)' / 3000, fr, noise3./abs(T)' / 3000, 'LineWidth', 2)
legend('Quantum noise', '6db squeezing', '6db anti-squeezing')
xlim([0.1, 1000])
xlabel('Frequency [Hz]')
ylabel('Sensitivity [Hz^{-1/2}]')
title('Dual recycled - zero detuning')
%%
% As expected, introducing squeezing improves the high frequency part, but
% the radiation pressure dominated low frequency region gets worse. Anti
% squeezing have exactly the opposite effect.
%
% Before it's too late, we save these results to a file, to compare them
% with the frequency dependent squeezing case lter on.