function ScatteredLight()
% Construct and simulate length sensing and control of LCGT
%
% == Usage ==
%
% Before running this function, review and edit the following files.
%
% (1) paramLCGT.m
%
% This function returns a structure p, which contains all the parameters
% necessary for the LSC simulation.
% Read the comments at the beginning of the file for details.
% Usually you only have to edit this file and files included herein.
%
% (2) modelLCGT.m
%
% This function constructs an Optickle model of LCGT taking the parameter
% structure p as an argument.
% If you want to modify the optical configuration of LCGT, edit this file.
%
% (3) determineDemodPhase.m
%
% This function is used to optimize the demodulation phases for various
% probes. It also generates vector sensing matrix for each probe.
% Edit this file if you want to change how the demodulation phases are
% optimized.
%

%{{{ Description

%
% This part adds scattered light fields to various parts of
% the interferometer.
%
% Scattering model:
% Assume that the scattering surface (baffles, tube walls etc)
% are moving at a particular frequency f by an amplitude a [m].
% Some stray light coming out of the interferometer is reflected by this
% surface.
% Then it will gain phase modulation of 
% dphi = 2*pi*a/lambda = k*a, 
% where lambda is the wavelength and k is wave number of the laser.
% A portion of this scattered light will be coupled back to the
% interferometer. The DC amplitude (not power) of the scattered 
% light coupled back to the interferometer is denoted Escl in the following.
% 
% We implement the above scattering model into the Optickle simulation by
% the following steps. First, we add a light source with amplitude
% Escl. The beam coming out of this light source is then phase modulated
% by k*a [rad]. After that, this light is injected to various parts of
% the interferometer, for example, the HR surface of ETMX.
%

%}}}

%{{{ ==== Initialization ====

%{{{ Path setup

%% Path setup

baseDir = pwd;
addpath(genpath([baseDir,'/Optickle']));
addpath(genpath([baseDir,'/tools']));
addpath(genpath([baseDir,'/Tests']));

%}}}

%{{{ Read parameters

%% Read parameters
addpath(genpath([pwd,'/parameters']));
ScatteredLightTF;

%}}}

%{{{ Create result directory

%% Create result directory

if p.DRSE
    resultDir = 'DRSE';
else
    resultDir = 'BRSE';
end
resultDir = [baseDir,'/results/',p.Name,'_',num2str(round(rem(now*86400,100000))),'/'];

mkdir(resultDir);

%}}}

%{{{ Flags

%% Flags

%Whether to generate auxiliary plots
osaplots = p.osaplots;
sweepplots = p.sweepplots;

%}}}

%{{{ Units constants

%% Units constants
nm=1e-9;
pm=1e-12;
MHz=1e6;
ppm=1e-6;

%}}}

%{{{ Close all figures

%% Close all figures
close all

%}}}

%}}}

%{{{ ===== Adjust IFO parameters =====

%{{{ Construct an Optickle model

%% Construct an Optickle model
[opt,n,l,pr,dr] = modelLCGT(p);

% == Extract some useful model information ==
drvNames=getDriveNames(opt);
Ndrive=opt.Ndrive;
Nlink=opt.Nlink;
lambda=opt.lambda;
c=opt.c;

%}}}

%{{{ *********** Check DC fields ***************

%{{{ Sweep

%% Sweep
pos=zeros(Ndrive,1);
[fDC, sigDC]=sweep(opt,pos);

%}}}

%{{{ Scale the input laser power to make the carrier power at the BS constant

%% Scale the input laser power to make the carrier power at the BS constant

%Carrier power at the BS
P_BS = abs(fDC(l.PR3toBS,find(p.vMod==0)))^2;

%Scale the input laser power
p.Pin = p.Pin * p.Pbs/P_BS;

%Reconstruct the Optickle Model
[opt,n,l,pr,dr] = modelLCGT(p);

%Recalculate the power
[fDC, sigDC]=sweep(opt,pos);
P_BS = abs(fDC(l.PR3toBS,find(p.vMod==0)))^2;

%Arm power
Parm = abs(fDC(l.ITMXtoETMX,find(p.vMod==0)))^2;

%}}}

%{{{ Adjust Homodyne Phase

%% Adjust Homodyne Phase

%{{{ Measure the AS power by arm mismatch

% No offset
opt = setPosOffset(opt, n.ETMX, 0); 
opt = setPosOffset(opt, n.ETMY, 0);
pos=zeros(Ndrive,1);
[fDC, sigDC]=sweep(opt,pos);
% Carrier power leaking out to the AS port by reflectivity mismatch
Pmismatch = abs(fDC(l.OMCTR, find(p.vMod==0))).^2;

%}}}

%{{{ Measure the AS power with arm offset

% Add offset
testArmOffset = 0.5*pm;
opt = setPosOffset(opt, n.ETMX, testArmOffset); 
opt = setPosOffset(opt, n.ETMY, -testArmOffset);

pos=zeros(Ndrive,1);
[fDC, sigDC]=sweep(opt,pos);

% Carrier power leaking out to the AS by the DARM offset
Pdc = abs(fDC(l.OMCTR, find(p.vMod==0))).^2 - Pmismatch;

%}}}

%{{{ Adjust the arm offset

% Necessary Pdc to achieve the target HD phase
PdcTarget = Pmismatch * tan(p.HDphase*pi/180)^2;
%PdcTarget = 100e-3;

% Arm offset necessary to realize the target HD phase
p.armOffset = p.DCReadout * sqrt(PdcTarget/(Pdc/(testArmOffset^2)));

%}}}

%}}}

%{{{ Check HD Phase

%% Check HD Phase
% Add offset
opt = setPosOffset(opt, n.ETMX, p.armOffset); 
opt = setPosOffset(opt, n.ETMY, -p.armOffset);

pos=zeros(Ndrive,1);
[fDC, sigDC]=sweep(opt,pos);

% Carrier power leaking out to the AS by the DARM offset
Pdc = abs(fDC(l.OMCTR, find(p.vMod==0))).^2 - Pmismatch;
% Homodyne phase
HDphase = atan(sqrt(Pdc/Pmismatch))*180/pi;

%}}}

%{{{ Power at detection ports

%% Power at detection ports

P_REFL1 = abs(fDC(l.REFL1,:)*conj(fDC(l.REFL1,:))');
P_REFL2 = abs(fDC(l.REFL2,:)*conj(fDC(l.REFL2,:))');
P_REFL3 = abs(fDC(l.REFL3,:)*conj(fDC(l.REFL3,:))');
P_REFL4 = abs(fDC(l.REFL4,:)*conj(fDC(l.REFL4,:))');

P_POP1 = abs(fDC(l.POP1,:)*conj(fDC(l.POP1,:))');
P_POP2 = abs(fDC(l.POP2,:)*conj(fDC(l.POP2,:))');

P_ASPO = abs(fDC(l.ASPO,:)*conj(fDC(l.ASPO,:))');
P_OMCREFL = abs(fDC(l.OMCREFL,:)*conj(fDC(l.OMCREFL,:))');
P_AS = abs(fDC(l.OMCTR,:)*conj(fDC(l.OMCTR,:))');
P_POX = abs(fDC(l.POX,:)*conj(fDC(l.POX,:))');
P_POY = abs(fDC(l.POY,:)*conj(fDC(l.POY,:))');

%}}}

%{{{ Adjust attenuators

%% Adjust attenuators

P_REFL_Max = max([P_REFL1,P_REFL2,P_REFL3,P_REFL4]);

if P_REFL_Max > p.PpdMax
   p.REFL_ATTN = 1 - p.PpdMax/P_REFL_Max; 
end

P_POP_Max = max([P_POP1,P_POP2]);
if P_POP_Max > p.PpdMax
   p.POP_ATTN = 1 - p.PpdMax/P_POP_Max; 
end

%Reconstruct the Optickle Model
[opt,n,l,pr,dr] = modelLCGT(p);

%}}}

%{{{ Check power at detection ports again

%% Check power at detection ports again

pos=zeros(Ndrive,1);
[fDC, sigDC]=sweep(opt,pos);

P_REFL1 = abs(fDC(l.REFL1,:)*conj(fDC(l.REFL1,:))');
P_REFL2 = abs(fDC(l.REFL2,:)*conj(fDC(l.REFL2,:))');
P_REFL3 = abs(fDC(l.REFL3,:)*conj(fDC(l.REFL3,:))');
P_REFL4 = abs(fDC(l.REFL4,:)*conj(fDC(l.REFL4,:))');

P_POP1 = abs(fDC(l.POP1,:)*conj(fDC(l.POP1,:))');
P_POP2 = abs(fDC(l.POP2,:)*conj(fDC(l.POP2,:))');

P_ASPO = abs(fDC(l.ASPO,:)*conj(fDC(l.ASPO,:))');
P_OMCREFL = abs(fDC(l.OMCREFL,:)*conj(fDC(l.OMCREFL,:))');
P_AS = abs(fDC(l.OMCTR,:)*conj(fDC(l.OMCTR,:))');
P_POX = abs(fDC(l.POX,:)*conj(fDC(l.POX,:))');
P_POY = abs(fDC(l.POY,:)*conj(fDC(l.POY,:))');

%}}}

%{{{ PRG

%% PRG

PRG=abs(fDC(l.PRMtoPR2, find(p.vMod==0))).^2/(abs(fDC(l.PRMbk, find(p.vMod==0))).^2);
Rarm=abs(fDC(l.ITMXtoBS, find(p.vMod==0))).^2/(abs(fDC(l.BStoITMX, find(p.vMod==0))).^2);

%}}}

%}}}

%{{{ ********** Optimize the demodulation phase **********

%{{{ Do optimization

%% Do optimization
mkdir([resultDir,'VectorSensingMatrices']);
determineDemodPhase(0, [resultDir,'VectorSensingMatrices/'], p);

%}}}

%{{{ Load the saved demodulation phases 

%% Load the saved demodulation phases 
demodPhaseLCGT;

%}}}

%}}}

%{{{ Reconstruct the Optickle model

%% Reconstruct the Optickle model
% % Construct the model with the optimized demodulation phase.
[opt,n,l,pr,dr] = modelLCGT(p);

% == Extract some useful model information ==
drvNames=getDriveNames(opt);
Ndrive=opt.Ndrive;
Nlink=opt.Nlink;
lambda=opt.lambda;
c=opt.c;

%}}}

%}}}

%{{{ === Default configuration without scattered light injection ===

%{{{ ******** Tickle ************

%{{{ tickle

%% Tickle

display('Start Tickling ...')
f=logspace(0,4,500);
tic
[fDC, sigDC, sigAC, mMech, noiseAC, noiseMech] = ptickle(opt, [], f);
toc
display('Finished.')
%}}}

%}}}

%{{{ Do DARM Calibration

%sigDARM is how much DARM signal is generated by moving DARM by 1m.
sigDARM0 = [0.5,-0.5]*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'ETMX','ETMY'}),:));
sigDARM0 = sigDARM0(:);

%}}}

%{{{ Plot Pure Quantum Sensitivity

%% Plot Pure Quantum Sensitivity
sigDARM = [0.5,-0.5]*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'ETMX','ETMY'}),:));
sigDARM = sigDARM(:);

noiseDARM = noiseAC(pr.AS_DC, :);
noiseDARM = noiseDARM(:);

pureQN = abs(noiseDARM./sigDARM/p.Larm);

h1=figure(1);
h=loglog(f, abs(noiseDARM./sigDARM/p.Larm));
set(gca,'FontSize',16);
grid on
title('Pure Quantum Strain Sensitivity','FontSize',16);
xlabel('Hz','FontSize',16);
h=ylabel('1/sqrt(Hz)','FontSize',16);

%}}}

%}}}

%{{{ === ETMX Scattered Light ===

%{{{ Construct Optickle model for SCL 

%Scattered light amplitude in sqrt(Watt)
p.EsETMX = sqrt(1e-7); 
p.EsPR2 = 0;
p.EsSR2 = 0; 
p.EsMICHX_PO = 0; 
k = 2*pi/opt.lambda; %Wave number

% % Construct the model with the optimized demodulation phase.
[opt,n,l,pr,dr] = modelLCGTSCL(p);

% == Extract some useful model information ==
drvNames=getDriveNames(opt);
Ndrive=opt.Ndrive;
Nlink=opt.Nlink;
lambda=opt.lambda;
c=opt.c;

%}}}

%{{{ ******** Tickle Again ************

%{{{ tickle

%% Tickle

display('Start Tickling ...')
f=logspace(0,4,500);
tic
[fDC, sigDC, sigAC, mMech, noiseAC, noiseMech] = ptickle(opt, [], f);
toc
display('Finished.')
%}}}

%}}}

%{{{ Do DARM Calibration

%sigDARM is how much DARM signal is generated by moving DARM by 1m.
sigDARM = [0.5,-0.5]*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'ETMX','ETMY'}),:));
sigDARM = sigDARM(:);
loglog(f,abs(sigDARM./sigDARM0));
grid on;
%}}}

%{{{ Plot Scattered Light Coupling

sigSCL = 2*k/p.EsETMX*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'PM_SCL_ETMX'}),:));
sigSCL = sigSCL(:);

SCLCoupling = abs(sigSCL./sigDARM0);

h1=figure(3);
h=loglog(f, SCLCoupling);
set(gca,'FontSize',16);
grid on
title('Scattered Light Coupling from ETMX','FontSize',16);
xlabel('Frequency [Hz]','FontSize',16);
h=ylabel('|G| [1/sqrt(W)]','FontSize',16);

M = [f(:), SCLCoupling(:)];
dlmwrite([resultDir, 'ETMX_SCL.dat'],M, '-append');

%}}}

%}}}

%{{{ === PR2 Scattered Light ===

%{{{ Construct Optickle model for SCL 

%Scattered light amplitude in sqrt(Watt)
p.EsETMX = 0; 
p.EsPR2 = sqrt(1e-7); 
p.EsSR2 = 0; 
p.EsMICHX_PO = 0; 
k = 2*pi/opt.lambda; %Wave number

% % Construct the model with the optimized demodulation phase.
[opt,n,l,pr,dr] = modelLCGTSCL(p);

% == Extract some useful model information ==
drvNames=getDriveNames(opt);
Ndrive=opt.Ndrive;
Nlink=opt.Nlink;
lambda=opt.lambda;
c=opt.c;

%}}}

%{{{ ******** Tickle Again ************

%{{{ tickle

%% Tickle

display('Start Tickling ...')
f=logspace(0,4,500);
tic
[fDC, sigDC, sigAC, mMech, noiseAC, noiseMech] = ptickle(opt, [], f);
toc
display('Finished.')
%}}}

%}}}

%{{{ Do DARM Calibration

%sigDARM is how much DARM signal is generated by moving DARM by 1m.
sigDARM = [0.5,-0.5]*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'ETMX','ETMY'}),:));
sigDARM = sigDARM(:);
loglog(f,abs(sigDARM./sigDARM0));
grid on;
%}}}

%{{{ Plot Scattered Light Coupling

sigSCL = 2*k/p.EsPR2*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'PM_SCL_PR2'}),:));
sigSCL = sigSCL(:);

SCLCoupling = abs(sigSCL./sigDARM0);

h1=figure(3);
h=loglog(f, SCLCoupling);
set(gca,'FontSize',16);
grid on
title('Scattered Light Coupling from PR2','FontSize',16);
xlabel('Frequency[Hz]','FontSize',16);
h=ylabel('|G| [1/sqrt(W)]','FontSize',16);

M = [f(:), SCLCoupling(:)];
dlmwrite([resultDir, 'PR2_SCL.dat'],M, '-append');

%}}}

%}}}

%{{{ === SR2 Scattered Light ===

%{{{ Construct Optickle model for SCL 

%Scattered light amplitude in sqrt(Watt)
p.EsETMX = 0; 
p.EsPR2 = 0;
p.EsSR2 = sqrt(1e-7); 
p.EsMICHX_PO = 0; 
k = 2*pi/opt.lambda; %Wave number

% % Construct the model with the optimized demodulation phase.
[opt,n,l,pr,dr] = modelLCGTSCL(p);

% == Extract some useful model information ==
drvNames=getDriveNames(opt);
Ndrive=opt.Ndrive;
Nlink=opt.Nlink;
lambda=opt.lambda;
c=opt.c;

%}}}

%{{{ ******** Tickle Again ************

%{{{ tickle

%% Tickle

display('Start Tickling ...')
f=logspace(0,4,500);
tic
[fDC, sigDC, sigAC, mMech, noiseAC, noiseMech] = ptickle(opt, [], f);
toc
display('Finished.')
%}}}

%}}}

%{{{ Do DARM Calibration

%sigDARM is how much DARM signal is generated by moving DARM by 1m.
sigDARM = [0.5,-0.5]*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'ETMX','ETMY'}),:));
sigDARM = sigDARM(:);
figure(1)
loglog(f,abs(sigDARM./sigDARM0));
grid on;
%}}}

%{{{ Plot Scattered Light Coupling

sigSCL = 2*k/p.EsSR2*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'PM_SCL_SR2'}),:));
sigSCL = sigSCL(:);

SCLCoupling = abs(sigSCL./sigDARM0);

h1=figure(2);
h=loglog(f, SCLCoupling);
set(gca,'FontSize',16);
grid on
title('Scattered Light Coupling from SR2','FontSize',16);
xlabel('Frequency[Hz]','FontSize',16);
h=ylabel('|G| [1/sqrt(W)]','FontSize',16);

M = [f(:), SCLCoupling(:)];
dlmwrite([resultDir, 'SR2_SCL.dat'],M, '-append');

%}}}

%}}}

%{{{ === MICHX_PO Scattered Light ===

%{{{ Construct Optickle model for SCL 

%Scattered light amplitude in sqrt(Watt)
p.EsETMX = 0; 
p.EsPR2 = 0;
p.EsSR2 = 0; 
p.EsMICHX_PO = sqrt(1e-7); 
k = 2*pi/opt.lambda; %Wave number

% % Construct the model with the optimized demodulation phase.
[opt,n,l,pr,dr] = modelLCGTSCL(p);

% == Extract some useful model information ==
drvNames=getDriveNames(opt);
Ndrive=opt.Ndrive;
Nlink=opt.Nlink;
lambda=opt.lambda;
c=opt.c;

%}}}

%{{{ ******** Tickle Again ************

%{{{ tickle

%% Tickle

display('Start Tickling ...')
f=logspace(0,4,500);
tic
[fDC, sigDC, sigAC, mMech, noiseAC, noiseMech] = ptickle(opt, [], f);
toc
display('Finished.')
%}}}

%}}}

%{{{ Do DARM Calibration

%sigDARM is how much DARM signal is generated by moving DARM by 1m.
sigDARM = [0.5,-0.5]*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'ETMX','ETMY'}),:));
sigDARM = sigDARM(:);
figure(1)
loglog(f,abs(sigDARM./sigDARM0));
grid on;
%}}}

%{{{ Plot Scattered Light Coupling

sigSCL = 2*k/p.EsMICHX_PO*squeeze(sigAC(pr.AS_DC, getDriveNumbers(opt, {'PM_SCL_MICHX_PO'}),:));
sigSCL = sigSCL(:);

SCLCoupling = abs(sigSCL./sigDARM0);

h1=figure(2);
h=loglog(f, SCLCoupling);
set(gca,'FontSize',16);
grid on
title('Scattered Light Coupling from MICHX','FontSize',16);
xlabel('Frequency[Hz]','FontSize',16);
h=ylabel('|G| [1/sqrt(W)]','FontSize',16);

M = [f(:), SCLCoupling(:)];
dlmwrite([resultDir, 'MICHX_SCL.dat'],M, '-append');

%}}}

%}}}