function runMI()
% Example code for simple Michelson interferometer
% Initiated by: Yuta Michimura

%% Clear everything
close all
clear all

%% Define some parameters
p.modelName = 'exampleMI'; % name of this model
p.Pin = 1;           % input power [W]
p.fmod1 = 15e6;      % modulation frequency [Hz]
p.g1 = 0.1i;         % modulation depth (imaginary number for phase modulation)
p.ETMChr = 1/1900;   % 1/RoC radius of curvature for ETM [1/m]
p.Lx = 3e3;          % x-arm length [m]
p.Ly = 3e3-3;        % y-arm length [m]

%% Construct optickle model and set some parameters, lists
opt = Optickle([-p.fmod1,0,p.fmod1]);

% Add laser source
%  [opt, sn] = addSource(opt, name, vArf, z0, z)
%  vArf - amplitudes of each RF component (Nrf x 1)
%  z0 - beam range = (waist size)^2 * pi / lambda
%  z - distance to waist (negative if beam is converging)
opt.addSource('Laser', [0,p.Pin,0]);

% Add RFmodulator (EOM)
%  [opt, sn] = addRFmodulator(opt, name, fMod, aMod)
%  name - name of this optic
%  fMod - modulation frequency
%  aMod - modulation index (imaginary for phase, real for amplitude)
opt.addRFmodulator('Mod1', p.fmod1, p.g1);

% Add mirrors
%  [opt, sn] = addMirror(opt, name, aio, Chr, Thr, Lhr, Rar, Lmd, Nmd)
%  aio - angle of incidence (in degrees)
%  Chr - curvature of HR surface (Chr = 1 / radius of curvature)
%  Thr - power transmission of HR suface
%  Lhr - power loss on reflection from HR surface
%  Rar - power reflection of AR surface
%  Nmd - refractive index of medium (1.45 for fused silica, SiO2)
%  Lmd - power loss in medium (one pass)
opt.addMirror('ETMX',0, p.ETMChr, 0, 0, 0, 0, 1.45);
opt.addMirror('ETMY',0, p.ETMChr, 0, 0, 0, 0, 1.45);

% Add beam splitter
%  [opt, sn] = addBeamSplitter(opt, name, aio, Chr, Thr, Lhr, Rar, Lmd, Nmd)
%  aio - angle of incidence (in degrees)
%  Chr - curvature of HR surface (Chr = 1 / radius of curvature)
%  Thr - power transmission of HR surface
%  Lhr - power loss on reflection from HR surface
%  Rar - power reflection of AR surface
%  Nmd - refractive index of medium (1.45 for fused silica, SiO2)
%  Lmd - power loss in medium (one pass)
opt.addBeamSplitter('BS',45, 0, 0.5, 0, 0, 0, 1.45);

% Add links
opt.addLink('Laser', 'out', 'Mod1', 'in',0);
opt.addLink('Mod1', 'out', 'BS', 'frA',0);
opt.addLink('BS', 'bkA', 'ETMX', 'fr', p.Lx);
opt.addLink('ETMX', 'fr', 'BS', 'bkB', p.Lx);
opt.addLink('BS', 'frA', 'ETMY', 'fr', p.Ly);
opt.addLink('ETMY', 'fr', 'BS', 'frB', p.Ly);

% Add REFL probes
opt.addSink('REFL',0);
opt.addLink('BS', 'frB', 'REFL','in',0);
%  opt = addReadout(opt, name, fphi, names)
%  name - sink name (also used as base name for probes in the readout)
%  fphi - demod frequency and phase (Nx2)
%  names - optional cell array of probe suffixes (instead of 1:N)
opt.addReadout('REFL', [p.fmod1,0]);

% Add AS probes
opt.addSink('AS',0);
opt.addLink('BS', 'bkB', 'AS','in',0);
opt.addReadout('AS', [p.fmod1,0]);

% some parameters
lambda = opt.lambda(1);

% DOF definitions
probeNames={'REFL_DC','REFL_I', 'REFL_Q','AS_DC','AS_I','AS_Q'};  % probes of interest
driveNames = {'ETMX','ETMY','BS'};  % drives of interest

%% Draw a graph for optickle model
opt2dot(opt, [p.modelName,'.pdf']); % GraphViz needs to be installed
    
%% Check power at each field, each probe
[fDC, sigDC] = opt.tickle([], []);
fprintf('DC fields (fDC matrix):\n');
showfDC(opt, fDC);
fprintf('\nProbes (sigDC matrix):\n');
showsigDC(opt, sigDC);

%% Sweep ETM in length (sweep function)
% define which mirror to sweep by how much
sweepDOF = 'ETMX';
range = lambda/2;
pos = zeros(opt.Ndrive, 1001);
x = linspace(range,-range,1001);
pos(opt.getDriveNum(sweepDOF),:) = x;

% do the sweep
[fDC, sigDC] = opt.sweep(pos);

% plot
figure()
for ii=1:length(probeNames)
    plot(x,sigDC(opt.getProbeNum(probeNames(ii)),:));
    hold on
end
legend(probeNames,'interpreter', 'none')
xlabel([sweepDOF,' displacement [m]']);
ylabel('signal [W]');
xlim([-range,range])
grid on

saveas(gcf,[p.modelName,sweepDOF,'sweep.png'])

%% Length sensing matrix (tickle function)
[fDC, sigDC, sigAC] = opt.tickle([], 10);
% [fDC, sigDC, sigAC] = opt.tickle([0,0,lambda/16,-lambda/16,0], 10);  % detuned
iDrv = cellfun(@opt.getDriveNum,driveNames);
iPrb = cellfun(@opt.getProbeNum,probeNames);
mSens = sigAC(iPrb,iDrv)

% use KAGRA tools to show sensing matrix better
fprintf('Length Sensing Matrix [W/m]:\n');
showSensingMatrix(opt,sigAC,probeNames,driveNames);