%
%
%

function make_mystery_plot

% frequency vector for power-law curves
f = logspace(1, 4, 1000)';

% some parameters
cav_pole = 340;
cav_cal = abs(1 + 1i * f / cav_pole);

% make basic curves
nCTN = 2e-20 * sqrt(30 ./ f);
nRPN = 13e-20 * (10 ./ f).^2;
nShot = 2.3e-20 * (cav_pole / 340) .* cav_cal;
nOscAM = 0.75e-20 .* cav_cal;

% other curves
nMystery = 0.9e-20 * (100 ./ f).^2;

% totals
nAll = [sqrt(nShot.^2 + nRPN.^2), nCTN, nOscAM, nMystery];
nAllUncorr = nShot;
nAllCorr = [nRPN, nCTN, nOscAM, nMystery];
nTotal = sqrt(sum(nAll.^2, 2));
nTotalCorr = sqrt(sum(nAllCorr.^2, 2));
nTotalUncorr = sqrt(sum(nAllUncorr.^2, 2));

% load DARM noises (from alog 25799)
nData = load('DARM_noise');
fD = nData.fDARM;
fp = 215; fz = 193;
data_cal = (fz/fp) * abs((1 + 1i * fD / fz)./(1 + 1i * fD / fp));
nData.dlDARM(:, 1) = nData.dlDARM(:, 1) .* data_cal;
nData.dlDARM(:, 2) = nData.dlDARM(:, 2) .* data_cal;
cav_cal_DARM = interpLogLog(f, cav_cal, nData.fDARM);

nUncorrDARM = interpLogLog(f, nTotalUncorr, nData.fDARM);
nTotalDARM = interpLogLog(f, nTotal, nData.fDARM);
nTotalCorrDARM = interpLogLog(f, nTotalCorr, nData.fDARM);
nCorr = asd_diff(nData.dlDARM(:, 1), nUncorrDARM);
nRes = asd_diff(nData.dlDARM(:, 1), nTotalDARM);
nResCorr = asd_diff(nData.dlDARM(:, 2), nTotalCorrDARM);

% plot
figure(1)
loglog(nData.fDARM, [nData.dlDARM, nRes], ...
f, [nAll, nTotal, nTotalCorr])
legend('Total', 'Correlated', 'Residual', ...
'Quantum', 'Coating', 'Osc AM', 'Mystery (1/f^2)', ...
'Sum', 'Sum Corr')
axis([20, 1000, 3e-21, 1e-18])
grid on
xlabel('frequency [Hz]')
ylabel('Noise [m / \surd Hz]')
title('Mystery Noise')

figure(2)
loglog(nData.fDARM, [nData.dlDARM, ...
asd_diff(nData.dlDARM(:, 1), nData.dlDARM(:, 2))./ cav_cal_DARM], ...
f, nShot ./ cav_cal)
legend('Total', 'Correlated', 'Residual', 'Shot')
axis([20, 1000, 1e-20, 1e-19])
grid on
xlabel('frequency [Hz]')
ylabel('Noise [m / \surd Hz]')
title('Mystery Noise')
end

function asd = asd_sum(varargin)
asd = sqrt(sum([varargin{:}].^2, 2));
end

function asd = asd_diff(a, b)
asd = sqrt(max(asd_sum(a).^2 - asd_sum(b).^2, 1e-30^2));
end