1144% Noise budget script for iKAGRA Michelson
%
% Author: Yuta Michimura

%% Add paths
clear all
close all
findNbSVNroot;  % find the root of Simulink NB
addpath(genpath([NbSVNroot 'Common/Utils']));
addpath(genpath([NbSVNroot 'Dev/Utils/']));

figdir = './results/';  % directory to save figures


%% Define some parameters, filters, and noises
% see the Simulink file (noiseModel) for each parameter
freq = logspace(-1,4,1000);
noiseModel = 'Michelson_aco';
IFO = 'K1';
site = 'kamioka';

% get online data (see http://klog.icrr.u-tokyo.ac.jp/osl/index.php?r=1357)
OnlineData = GWData;
GPStime = 1144857600;   % time to get online data 
duration = 55;  % data get time duration
freqsamp = 16384;   % sampling frequency
freqmin = 0.1;  % minimum frequency (= frequency reslution)
Ndata = ceil(freqsamp/freqmin);
ndata = 2^nextpow2(Ndata);
duration2 = 505;  % data get time duration
freqsamp2 = 2048;   % sampling frequency
freqmin2 = 0.1;  % minimum frequency (= frequency reslution)
Ndata2 = ceil(freqsamp2/freqmin2);
ndata2 = 2^nextpow2(Ndata2);

% get LiveParts parameters
liveParts(noiseModel,GPStime,duration,freq);


% Constants
c = 299792458.; % speed of light in vacuum (m/s)
hP = 6.62606957e-34;	% Planck's Constant [J*s]
lamb = 1064e-9; % laser wavelength (m)
nu = c/lamb;    % laser frequency (Hz)

% Dictionaries
Filt = containers.Map;  % digital filter transfer functions [cnts/cnts]
Susp = containers.Map;  % suspension actuation transfer functions [m/N]
Noise = containers.Map; % all sorts of noises
PdResp = containers.Map;    % PD rensponse [A/W]
TransImp = containers.Map;  % PD trans impedance [V/A]
Elec = containers.Map;  % electronics gain [V/V]
Coil = containers.Map;  % coil gain [N/V]


% ADC and DAC
ADC_V2C = 2^16/40;  % ADC factor [cnts/V] 
DAC_C2V = 20/2^16;  % DAC factor [V/cnts]


% PD signal chain and noises
for pd = {'REFL'}
    pd = char(pd);
    PdResp([pd,'_DC']) = 0.3; % A/W (Thorlabs PDA100A)
    TransImp([pd,'_DC']) = 0.75e3; % V/A (0dB setting; see https://www.thorlabs.co.jp/thorcat/13000/PDA100A-Manual.pdf)
    PdResp([pd,'_RF']) = 0.8; % A/W (PD is PerkinElmer C30642G (see JGW-D1201280-v2) / )
    TransImp([pd,'_RF']) = 200/0.8; % V/A (see http://granite.phys.s.u-tokyo.ac.jp/internal/2012_m_ishigaki.pdf)
    noise = load('./Noises/20160402/LSC-REFL_PDA1_DC_OUT_DQ_dark.txt');  % dark noise measured when PR3 misaligned
    Noise([pd,'_DC_PD']) = interp1(noise(:,1), noise(:,2), freq, 'linear');  % PD dark noise in [counts/rtHz];
    noise = load('./Noises/20160407/LSC-REFL_PDA1_RF17_Q_dark.txt');  % dark noise measured when MI misaligned
    Noise([pd,'_RF_PD']) = interp1(noise(:,1), noise(:,2), freq, 'linear');  % RF PD dark noise in [counts/rtHz];
end
DeModGain = 10.9; % Demodulation gain (see LIGO-F1100004-v4)
noise = load('./Noises/adc_noise.txt');
%[data,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-ADCNOISE_GND_OUT_DQ'); % use online data 
%[noisep,noisef] = pwelch(data,hanning(ndata),[],ndata,freqsamp); % calculate power spectral density
%noise = [noisef sqrt(noisep)]; 
Noise('ADC') = interp1(noise(:,1), noise(:,2), freq, 'linear');   % ADC noise [cnts/rtHz]
Noise('ADC') = 0; % since ADC noise is included in PD dark noise now


% Intensity noise (see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1209)
[data,time] = OnlineData.fetch(GPStime,duration,'K1:IMC-MCL_TRANS_OUT_DQ'); % use online data 
[noisep,noisef] = pwelch(data,hanning(ndata),[],ndata,freqsamp); % calculate power spectral density
noise = [noisef sqrt(noisep)]; 
Noise('RIN') = interp1(noise(:,1), noise(:,2), freq, 'linear')/0.4;   % relative intensity noise [/rtHz]
INT_NOISE_COUP = 0.0096; % intensity noise coupling to REFL DC PD [W] fitted to fit with measured spectrum


% IFO Optical responses and noises
las = 3.3299;	% Schnupp asymmetry [m]
P0 = 250e-3;	% input power to BS	[W]
Ppd = 8e-3;	% power at REFL PD	[W] ~3000 counts
Freq2MICH = las/nu; % Laser frequency noise to MICH coupling [m/Hz]
noise = load('./Noises/IMCOOL.dat');
Noise('LaserFreq') = interp1(noise(:,1), noise(:,2), freq, 'linear'); % Laser frequency noise [Hz/rtHz]; calculated with /kagranoisebudget/FSS/KAGRA_FSS.mdl model
%MICHOpticalGain = -2*pi*nu/c*P0;  % W/m (theoretical)
MICHOpticalGainDC = 3.4e10/ADC_V2C/PdResp('REFL_DC')/TransImp('REFL_DC'); % measured value (see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1340); this is nominal value, adjusted later using calibration peak
Noise('REFL_DC_SHOT') = sqrt(2*hP*nu*Ppd); % W/rtHz
MICHOpticalGainRF = 5.3e10/ADC_V2C/PdResp('REFL_RF')/TransImp('REFL_RF')/DeModGain; % measured value (see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1390)
Noise('REFL_RF_SHOT') = sqrt(2*hP*nu*Ppd); % W/rtHz  THIS IS A FAKE (power at PD should be measured / Shot noise should be calculated in a case of dark fringe) 


% Suspension chain and noises
for mir = {'BS','ETMX','ETMY'}
    mir = char(mir);
    Coil([mir,'_TM']) = 1e-4;   % N/V see email from K. Yamamoto on Mar 22 (to be confirmed)
    Noise([mir,'_TM_Coil']) = 57e-9./freq+2.1e-9; % high power coil driver noise [V/rtHz]; See Figure 3 of https://dcc.ligo.org/LIGO-T080014
    Elec([mir,'_TM_Whitening']) = 1; % no whitening for iKAGRA
    Elec([mir,'_TM_DeWhitening']) = 1; % no dewhitening for iKAGRA
end
Filt('BS_ISCINF') = LvFilt_BS_ISCINF.gain;  % TO BE REPLACED WITH LIVE PARTS
Filt('ETMX_ISCINF') = LvFilt_ETMX_ISCINF.gain;  % TO BE REPLACED WITH LIVE PARTS
Filt('ETMY_ISCINF') = LvFilt_ETMY_ISCINF.gain;  % TO BE REPLACED WITH LIVE PARTS

% TM act to TM [m/N] see http://klog.icrr.u-tokyo.ac.jp/osl/?r=1340
Susp('BS_TM_TM') = myzpk('zpk',[],[0.87 3.5],1e-2);  % THIS IS FAKE
Susp('ETMX_TM_TM') = myzpk('zpk',[],[0.92 5.1],5.1e-6/1e-4);
Susp('ETMY_TM_TM') = myzpk('zpk',[],[0.94 4.6],3.3e-6/1e-4);
noise = load('./Noises/adc_noise.txt');
Noise('DAC') = interp1(noise(:,1), noise(:,2), freq, 'linear')*DAC_C2V/2;   % DAC noise [V/rtHz]; noise level is OK, but fake (see http://gwwiki.icrr.u-tokyo.ac.jp/JGWwiki/CLIO/Tasks/DigitalControl/Caltech_setup?action=AttachFile&do=get&target=analog_system_investigation.pdf)

%Seis_TypeC = 8e-9*freq.^-2./(1+(2./freq).^-8);
%Seis_TypeB = 1.5e-11*freq.^-9;
%Seis_TypeBp=2.3e-10*freq.^-7;
%Seis_TypeBpfixed=2e-10*freq.^-5;

% seismic motion to mirror motion [m/m]
Susp('BS_SEIS_TM') = myzpk('zpk',[],[0.93 5.0; 2.8 3.0],1);    % rough caculation
Susp('ETMX_SEIS_TM') = myzpk('zpk',[],[0.93 5.0; 2.8 3.0],1);    
Susp('ETMY_SEIS_TM') = myzpk('zpk',[],[0.93 5.0; 2.8 3.0],1);    

% seismic vibration data from seismometer 
[data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-BS_SEIS_NS_SENSINF_OUT_DQ');    % in unit of [m/s]
[noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density
noise = [noisef sqrt(noisep)./(2*pi*noisef)];  
Noise('BS_SEIS') = interp1(noise(:,1), noise(:,2), freq, 'linear');

%[data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-EX_SEIS_WE_SENSINF_OUT_DQ');    % in unit of [m/s]
[data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-BS_SEIS_WE_SENSINF_OUT_DQ');    % in unit of [m/s]
[noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density
noise = [noisef sqrt(noisep)./(2*pi*noisef)];  
Noise('ETMX_SEIS') = interp1(noise(:,1), noise(:,2), freq, 'linear');

%[data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-EY_SEIS_NS_SENSINF_OUT_DQ');    % in unit of [m/s]
[data,time] = OnlineData.fetch(GPStime,duration,'K1:PEM-BS_SEIS_NS_SENSINF_OUT_DQ');    % in unit of [m/s]
[noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density
noise = [noisef sqrt(noisep)./(2*pi*noisef)]; 
Noise('ETMY_SEIS') = interp1(noise(:,1), noise(:,2), freq, 'linear');


% coupling from mirror tilt(Pitch/Yaw) fluctuation
BeamMissAlignP = 1e-2;  % [m] beamspot offset from the center of mirror (fitted to fit measured spectrum by using peak height) 
[data,time] = OnlineData.fetch(GPStime,duration2,'K1:VIS-BS_TM_OPLEV_PIT_OUT_DQ'); % mirror tilt fluctuation (OpLev signal; [urad])
[noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density
noise = [noisef sqrt(noisep)*1e-6]; % [rad/rtHz]
Noise('BS_PIT') = interp1(noise(:,1), noise(:,2), freq, 'linear');

BeamMissAlignY = 1e-2;  % [m] beamspot offset from the center of mirror (fitted to fit measured spectrum by using peak height) 
[data,time] = OnlineData.fetch(GPStime,duration2,'K1:VIS-BS_TM_OPLEV_YAW_OUT_DQ'); % mirror tilt fluctuation (OpLev signal; [urad])
[noisep,noisef] = pwelch(data,hanning(ndata2),[],ndata2,freqsamp2); % calculate power spectral density
noise = [noisef sqrt(noisep)*1e-6]; % [rad/rtHz]
Noise('BS_YAW') = interp1(noise(:,1), noise(:,2), freq, 'linear');


% LSC filters and matrices
%FiltGain = 50;   % filter gain
%Filt('LSC_MICH') = FiltGain*myzpk('zpk',[70;0.2],[2000;2000;0.02],10); % TO BE REPLACED WITH LIVE PARTS
Filt('LSC_MICH') = LvFilt_LSC_MICH.gain * LvFilt_LSC_MICH.fm1 * LvFilt_LSC_MICH.fm2 * (-1); % factor -1 is for adjusting +/- (I don't know where this factor comes from. Please find it.)  


% AA and AI filters
% Anaglog AA/AI: 3rd-order Butterworth at 10kHz, notch at 65535Hz (LIGO-T070038, LIGO-D070081)
% Digital AA/AI: /opt/rtcds/rtscore/release/src/fe/controller.c
[z,p,k]=butter(3,2*pi*1e4,'low','s');
Elec('Digital_AA') = zpk(z,p,k);
Elec('Analog_AA') = zpk(z,p,k);
Elec('Digital_AI') = zpk(z,p,k);
Elec('Analog_AI') = zpk(z,p,k);


% Whitening and Dewhitening filters
Elec('REFL_DC_Whitening') = 1; % no whitening for iKAGRA
Elec('REFL_DC_DeWhitening') = 1; % no dewhitening for iKAGRA
Elec('REFL_RF_Whitening') = 1; % no whitening for iKAGRA
Elec('REFL_RF_DeWhitening') = 1; % no dewhitening for iKAGRA


% Measured noise at K1:LSC_MICH_OUT_DQ
%noise = load('./Noises/20160407/LSC-MICH_OUT_DQ.txt');
[data,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-MICH_OUT_DQ'); % use online data 
[noisep,noisef] = pwelch(data,hanning(ndata),[],ndata,freqsamp); % calculate power spectral density
noise = [noisef sqrt(noisep)]; 
MeasuredNoise = interp1(noise(:,1), noise(:,2), freq, 'linear');


% UGF servo
[data,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-UGF_SERVO_OUT_DQ');
UGFSERVO_GAIN = 1/mean(data);

% Accorstic noise
data = importdata('./Noises/20160430aco/excitation_test_spectrum.txt');
ff = data(:,1);
acc = data(:,3);%m/s^2
data = importdata('./Noises/20160430aco/excitation_test_tf2MICH_OUT.txt');
ff2 = data(:,1);
tf = data(:,2);
tf = tf;
Noise('BS_ACOUSTIC') = interp1(ff,tf.*acc,freq,'linear');

% calibration 
fcal = 80;   % calibration line frequency
UGFAdjustFactor = OLTCALIB_MTRX(1);  % get value from Live parts
[dataAr,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_A1_REAL_OUT_DQ'); % demodulated calibration line after injection (real part)
[dataAi,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_A1_IMAG_OUT_DQ'); % demodulated calibration line after injection (imaginary part)
[dataBr,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_B1_REAL_OUT_DQ'); % demodulated calibration line before injection (real part)
[dataBi,time] = OnlineData.fetch(GPStime,duration,'K1:LSC-CAL_B1_IMAG_OUT_DQ'); % demodulated calibration line before injection (imaginary part)
FBAfterInj = sqrt(mean(dataAr)^2+mean(dataAi)^2);     % use LSC-LKIN_CAL1_A_REAL and IMAG
FBBeforeInj = sqrt(mean(dataBr)^2+mean(dataBi)^2);    % use LSC-LKIN_CAL1_B_REAL and IMAG
OLTF_cal = FBBeforeInj/FBAfterInj;   % OpenLoop gain / calibration factor

% DOF Switches
MICH_SW = 1;

%% Plot nominal openloop transfer function with measured one
% set(0, 'DefaultAxesFontSize',10);
% MICH_SW = 0;    % turn off the loop to measure OLTF
% [A,B,C,D]=linmod(noiseModel);
% systm=ss(A,B,C,D);
% OLTF=systm(1,1,:);
% meas=load('./TransferFunctions/20160409/MICHOLTF-UGFservo.txt');
% meas(:,3)=angle(-exp(i*meas(:,3)/180*pi))/pi*180;
% figure(100)
% plotbode(freq,OLTF);
% subplot(2,1,1)
% loglog(meas(:,1),meas(:,2),'.');
% legend({'NB model','Measured 20160409'})
% loglog([freq(1),freq(end)],[1,1],'k')
% ylim([1e-6,1e6]);
% xlim([freq(1),freq(end)]);
% set(gca,'YTick',logspace(-6,6,7));
% set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1));
% title(sprintf('MICH OLTF (uncalibrated, optical gain: %10.2e counts/m)',abs(MICHOpticalGainRF)*ADC_V2C*PdResp('REFL_RF')*TransImp('REFL_RF')*DeModGain));
% subplot(2,1,2)
% loglog(meas(:,1),meas(:,3),'.');
% xlim([freq(1),freq(end)]);
% set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1));
% saveas(gcf,[figdir,'MICHOLTF.png'])
% 
% MICH_SW = 1;    % turn on the loop

%% Plot actual openloop transfer function using calibration line
meas=load('./TransferFunctions/20160409/MICHOLTF-UGFservo.txt');    % measured OLTF data
meas(:,3)=angle(-exp(i*meas(:,3)/180*pi))/pi*180;

MICH_SW = 0;    % turn off the loop to measure OLTF
[A,B,C,D]=linmod(noiseModel);
systm=ss(A,B,C,D);
OLTF=systm(1,1,:);
[OLTF_uncal,~]=bode(OLTF,2*pi*fcal);    % theoretical openloop transfer function gain at fcal
MICHOpticalGainRF = MICHOpticalGainRF * OLTF_cal/(OLTF_uncal/UGFSERVO_GAIN);    % calibrate optical gain using calibration line

MICH_SW = 0;    % turn off the loop to measure OLTF
[A,B,C,D]=linmod(noiseModel);
systm=ss(A,B,C,D);
OLTF=systm(1,1,:);

set(0, 'DefaultAxesFontSize',10);
figure(101)
plotbode(freq,OLTF);
subplot(2,1,1)
loglog(meas(:,1),meas(:,2),'.');
legend({'NB model','Measured 20160409'})
loglog([fcal],[1/UGFAdjustFactor],'k+')  % calibration point
loglog([freq(1),freq(end)],[1,1],'k')
ylim([1e-6,1e6]);
xlim([freq(1),freq(end)]);
set(gca,'YTick',logspace(-6,6,7));
set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1));
title(sprintf('MICH OLTF (calibrated, optical gain: %10.2e counts/m)',abs(MICHOpticalGainRF)*ADC_V2C*PdResp('REFL_RF')*TransImp('REFL_RF')*DeModGain));
subplot(2,1,2)
loglog(meas(:,1),meas(:,3),'.');
xlim([freq(1),freq(end)]);
set(gca,'XTick',logspace(log10(freq(1)),log10(freq(end)),log10(freq(end))-log10(freq(1))+1));
saveas(gcf,[figdir,'MICHOLTF_calibrated.png'])

MICH_SW = 1;    % turn on the loop

%% Compute noises and save cache
% Compute noises
[noises, sys] = nbFromSimulink(noiseModel, freq, 'dof', 'MICH');

% save cached outputs
saveFunctionCache();

%% Plot noise budget
set(0, 'DefaultAxesFontSize',10);
disp('Plotting noises')
nb = nbGroupNoises(noiseModel, noises, sys);
nb.sortModel();
matlabNoisePlot(nb);

figure(2)
ylim([1e-17,1e-5]);
xlim([1e-1,1e4]);
set(gca,'YTick',logspace(-17,-5,13));
set(gca,'XTick',logspace(-1,4,6));
saveas(gcf,[figdir,'MICHNoiseBudgetActuator_aco.png'])

figure(4)
ylim([1e-17,1e-5]);
xlim([1e-1,1e4]);
set(gca,'YTick',logspace(-17,-5,13));
set(gca,'XTick',logspace(-1,4,6));
saveas(gcf,[figdir,'MICHNoiseBudgetSensor_aco.png'])
%%
figure(1)
ylim([1e-17,1e-5]);
xlim([1e-1,1e4]);
set(gca,'YTick',logspace(-17,-5,13));
set(gca,'XTick',logspace(-1,4,6));
% legend('Location','notheast')
saveas(gcf,[figdir,'MICHNoiseBudget_aco.png'])