function Noise = make_Noise( freq, Servo_0dB, unable_noise_list, varargin )
%UNTITLED22 Summary of this function goes here
%   Detailed explanation goes here
Noise = containers.Map;
Freq = containers.Map;
Abs = containers.Map;
Phase = containers.Map;
Servo = containers.Map;
TTFSS_Elec = TTFSS_filters(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NPRO frequency noise [Hz/rHz]
Noise('NPRO') = 10000./freq;
Noise('freq_IMCpath') = zeros(size(freq));
Noise('freq_FSSpath') = zeros(size(freq));
% KOACH on
if mnismember(varargin, 'KOACH_off')
for ii = 1:5
    mnimportdata(['N1226_' num2str(ii) '.txt'],Freq,Abs,Phase,['NOISE_RefCavPD' num2str(ii)],'sp');
end
%combine data
mncombine(Freq,Abs,'NOISE_RefCavPD',5);
% dvide by SR560 gain
Abs('NOISE_RefCavPD') = Abs('NOISE_RefCavPD')/-TTFSS_Elec('common_gain1');
else
% KOACH off
for ii = 1:5
    mnimportdata(['N1226_' num2str(ii+5) '.txt'],Freq,Abs,Phase,['NOISE_RefCavPD' num2str(ii)],'sp');
end
%combine data
mncombine(Freq,Abs,'NOISE_RefCavPD',5);
% dvide by SR560 gain
Abs('NOISE_RefCavPD') = Abs('NOISE_RefCavPD')/-TTFSS_Elec('common_gain1');
end

% load('NPRO_noise.mat')
% Noise('NPRO') = mnNaN2zero(interp1(NPRO_noise('ff'),NPRO_noise('Abs'),freq));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ADC noise [cnts/rHz]
% data = importdata('N1102_16.txt');
% Freq('ADCnoise') = data(:,1);
% Abs('ADCnoise') = data(:,2);
% loglog(Freq('ADCnoise'),Abs('ADCnoise'),Freq('ADCnoise'),Freq('ADCnoise').^(-3.3/5)*1e-2+5e-3)
Noise('ADC') = freq.^(-3.3/5)*1e-2+5e-3;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Servo Noise
% Common mode servo NOISE
% IN1 terminated
% same configuration with the locking secuence
% used SR560 as pre-AMP
% SR560 : AC couple, gain 100, BPF
%importdata
for ii = 1:5
    mnimportdata(['N0901_6' num2str(ii) '.txt'],Freq,Abs,Phase,['NOISE_CMS' num2str(6-ii)],'sp');
end
%combine data
mncombine(Freq,Abs,'NOISE_CMS',5);
% dvide by SR560 gain
Abs('NOISE_CMS') = Abs('NOISE_CMS')/100;
% servo model
IN1GAIN = mn_CMSgain(-10);
Servo('CMS_AOM') = Servo_0dB('CMS_AOM')*10^(IN1GAIN/20);
servo_abs = mnbode_list(Servo('CMS_AOM'),Freq('NOISE_CMS'));
inputnoise_cms = Abs('NOISE_CMS')./servo_abs;
Noise('CMS_IMC') = mnNaN2zero(interp1(Freq('NOISE_CMS'),inputnoise_cms,freq));
Noise('CMS_CARM') = mnNaN2zero(interp1(Freq('NOISE_CMS'),inputnoise_cms,freq));

% TTFSS noise
for ii = 1:5
    mnimportdata(['N0728_' num2str(ii+6) '.txt'],Freq,Abs,Phase,['NOISE_TTFSS' num2str(ii)],'sp');
end
%combine data
mncombine(Freq,Abs,'NOISE_TTFSS',5);
% dvide by SR560 gain
Abs('NOISE_TTFSS') = Abs('NOISE_TTFSS')/100;
% servo model
fast_pol = 1;
CG = 30;
FG = 30;
Servo('TTFSS_PZT') = Servo_0dB('TTFSS_PZT')*10^(CG/20)*10^(FG/20)*fast_pol;
servo_abs = mnbode_list(Servo('TTFSS_PZT'),Freq('NOISE_TTFSS'));
inputnoise_ttfss = Abs('NOISE_TTFSS')./servo_abs;
Noise('TTFSS') = mnNaN2zero(interp1(Freq('NOISE_TTFSS'),inputnoise_ttfss,freq));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Actuator noise
% AOM actuator noise (driver input reffered)

Noise('act_AOM') = 3e-8*ones(size(freq));
if mnismember(varargin,'LIGO_VCO')
    Noise('act_AOM') = 3e-9*ones(size(freq));
end
Noise('act_PZT') = zeros(size(freq));
Noise('act_EOM') = zeros(size(freq));
Noise('act_TEMP') = zeros(size(freq));
Noise('act_IMCsus') = zeros(size(freq));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cavity noise
% IMC seismic motion [m/rHz]
% Noise('seismic_IMC') = 1e-9./(freq.^2+1);
% mnimportdata('kamiokaNoisy.txt',Freq,Abs,Phase,'seis','sp');
% Noise('Seis') = mnNaN2zero(interp1(Freq('seis'),Abs('seis'),freq));

% CARM seismic motion [m/rHz]
Noise('seismic_CARM') = 0;

% Refcav motion [m/rHz]
% if ~exist('RefCav_Noise.mat','file') || mnismember(varargin, 'RefCav_unknown')
%     Noise('motion_RefCav') = zeros(size(freq));
% elseif mnismember(varargin, 'KOACH_off')
%     load('RefCav_Noise.mat')
%     Noise('motion_RefCav') =...
%         mnNaN2zero(interp1(RefCav_Noise('ff_KOACH_off'),RefCav_Noise('KOACH_off'),freq));
% else
%     load('RefCav_Noise.mat')
%     Noise('motion_RefCav') =...
%         mnNaN2zero(interp1(RefCav_Noise('ff_KOACH_on'),RefCav_Noise('KOACH_on'),freq));
% end
data = importdata('kamiokaMNM.dat');
temp_ff = data(:,1);
temp_noise = data(:,2);
Noise('motion_RefCav') = mnNaN2zero(interp1(temp_ff,temp_noise,freq))/1000;

% IMC seismic [m/rHz]
if ~exist('IMC_Seismic.mat','file')
    Noise('seismic_IMC') = zeros(size(freq));
elseif mnismember(varargin,'kamioka_spectrum')
    data = importdata('kamiokaMNM.dat');
    temp_ff = data(:,1);
    temp_noise = data(:,2);
    Noise('seismic_IMC') = mnNaN2zero(interp1(temp_ff,temp_noise,freq));
else
    load('IMC_Seismic.mat')
    temp_noise =...
        (interp1(Seismic_IMC('ff'),Seismic_IMC('Abs'),freq));
    for ii = 1:length(freq)
        if isnan(temp_noise(ii))
            temp_noise(ii) = 2e-10./(freq(ii)/+1).^2;
        end
    end
    Noise('seismic_IMC') = temp_noise;
end

% residual gas noise
% see handbook
V0 = 2.24e-3;%volume of 1mol air
A0 = 6.02e23;%Avogadro's number
kb = 1.38e-23;%Boltzmann constant
mO2 = 32e-3/A0;%mass of O2
mN2 = 28e-3/A0;%mass of N2


if mnismember(varargin,'Pressure_RefCav')
    p = varargin{mnismember(varargin,'Pressure_RefCav')+1};
else
    p = 1;
end

if mnismember(varargin,'T_PSL')
    T = varargin{mnismember(varargin,'T_PSL')+1};
else
    T = 300;
end

u0O2 = sqrt(3*kb*T/2/mO2);%average velosity of O2
u0N2 = sqrt(3*kb*T/2/mN2);%average velosity of N2
n0N2 = 1.000297;%refraction index
n0O2 = 1.000272;%refraction index
L = 10e-2;%length of RefCav
p0 = 101325;%standard pressure (1atm)
T0 = 273.15;%standard temperature


dx_O2 = sqrt(8*sqrt(2*L)/sqrt(pi)*(n0O2-1)^2/((A0/V0)*u0O2*sqrt(1064e-9))*(p/p0)*(T/T0)^(3/2));
dx_N2 = sqrt(8*sqrt(2*L)/sqrt(pi)*(n0N2-1)^2/((A0/V0)*u0N2*sqrt(1064e-9))*(p/p0)*(T/T0)^(3/2));

Noise('residual_gas_RefCav') = sqrt(dx_O2^2*0.2+dx_N2^2*0.8)*ones(size(freq));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cavity shot noise (see LIGO-T980045)
if mnismember(varargin,'P_RefCav_REFL')
    P_RefCav = varargin{mnismember(varargin,'P_RefCav_REFL')+1};
else
    P_RefCav = 1e-3;
end
if mnismember(varargin,'P_IMC_REFL')
    P_IMC = varargin{mnismember(varargin,'P_IMC_REFL')+1};
else
    P_IMC = 1e-3;
end
if mnismember(varargin,'P_CARM_REFL')
    P_CARM = varargin{mnismember(varargin,'P_CARM_REFL')+1};
else
    P_CARM = 1e-3;
end
if mnismember(varargin,'Opgain_RefCav')
    H_RefCav = varargin{mnismember(varargin,'Opgain_RefCav')+1};
else
    H_RefCav = 2e-6;
end
if mnismember(varargin,'Opgain_IMC')
    H_IMC = varargin{mnismember(varargin,'Opgain_IMC')+1};
else
    H_IMC = 7e-6;
end
if mnismember(varargin,'Opgain_CARM')
    H_CARM = varargin{mnismember(varargin,'Opgain_CARM')+1};
else
    H_CARM = 7e-6;
end
finesse_RefCav = 1.1e4;
finesse_IMC = 460;
finesse_CARM = 1e4;
L_RefCav = 10e-2;
L_IMC = 53;
L_CARM = 3e3;

Noise('shot_IMC') = 1.2e-5*H_IMC*(1e4/finesse_IMC)*(20e-2/L_IMC)*sqrt(500e-3/P_IMC)*ones(size(freq));
Noise('shot_RefCav') = 1.2e-5*H_RefCav*(1e4/finesse_RefCav)*(20e-2/L_RefCav)*sqrt(500e-3/P_RefCav)*ones(size(freq));
Noise('shot_CARM') = 1.2e-5*H_CARM*(1e4/finesse_CARM)*(20e-2/L_CARM)*sqrt(500e-3/P_CARM)*ones(size(freq));


Noise('scat') = mnNaN2zero(interp1(Freq('NOISE_RefCavPD'),Abs('NOISE_RefCavPD'),freq))/H_RefCav;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sensor noise
% Noise('sensor_IMC') = zeros(size(freq));
data = importdata('N1206_2.txt');
temp_freq = data(:,1);
shutter = 'open';
if strcmp(shutter , 'open')
    temp_sp = data(:,2);
else
    temp_sp = data(:,3);
end

[temp_ff1, temp_abs1] = mntrim(0,10,temp_freq,temp_sp);
[temp_ff2, temp_abs2] = mntrim(10.1,1e9,temp_freq,temp_sp); 
temp_ff2 = downsample(temp_ff2,50);
temp_abs2 = downsample(temp_abs2,50);
ff = [temp_ff1; temp_ff2];
sp = [temp_abs1; temp_abs2];

Servo = containers.Map;
% Servo model
IN1GAIN = mn_CMSgain(25);
K1IMC_MCL_SERVO_GAIN = -40;
Servo('CMS_SUS') = Servo_0dB('CMS_SUS_analog_susonly')*Servo_0dB('CMS_SUS_digital_susonly')...
    *10^(IN1GAIN/20)*K1IMC_MCL_SERVO_GAIN;

servo_model = mnbode_list(Servo('CMS_SUS'),ff,'c');
tempdata = sp./abs(servo_model);
[fitfreq,fitdata] = mntrim(1e-1,1e3,ff,tempdata);
% fit
[xData, yData] = prepareCurveData( fitfreq, fitdata );

% Set up fittype and options.
ft = fittype( 'a/((f/fc)^c+1)+b', 'independent', 'f', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Robust = 'Bisquare';
opts.StartPoint = [6.284e-6 1.1e-7 1.8 8];

% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
if mnismember(varargin,'IMC_sensor_fit')

Noise('sensor_IMC') = fitresult.a./((freq/fitresult.fc).^fitresult.c+1)+fitresult.b;

else
    Noise('sensor_IMC') = mnNaN2zero(interp1(ff,tempdata,freq));
end


if ~exist('Sensor_Noise_RefCav.mat','file') || mnismember(varargin, 'sensor_unknown')
    Noise('sensor_RefCav') = zeros(size(freq));
elseif mnismember(varargin, 'KOACH_off')
    load('Sensor_Noise_RefCav.mat')
    temp_noise =...
        interp1(Sensor_Noise_RefCav('ff_KOACH_off'),Sensor_Noise_RefCav('KOACH_off'),freq);
    for ii = 1:length(freq)
        if isnan(temp_noise(ii))
            temp_noise(ii) = 8.071e-5*(15./(freq(ii)/0.1+1).^0.6./(freq(ii)/2000+1)./(freq(ii)/3000+1).*(freq(ii)/700+1))+fitresult.b;
        end
    end
    Noise('sensor_RefCav') = temp_noise;
else
    load('Sensor_Noise_RefCav.mat')
    temp_noise =...
        interp1(Sensor_Noise_RefCav('ff_KOACH_on'),Sensor_Noise_RefCav('KOACH_on'),freq);
    for ii = 1:length(freq)
        if isnan(temp_noise(ii))
            temp_noise(ii) = 8.761e-4.*(1./(freq(ii)/8+1)./(freq(ii)/8e3+1))+fitresult.b;
        end
    end
    Noise('sensor_RefCav') = temp_noise;
end
if mnismember(varargin, 'KOACH_off')
[ff,sp] = import_1206data('','sensor RefCav','KOACH_off');
else
[ff,sp] = import_1206data('','sensor RefCav','KOACH_on');
end
Noise('sensor_RefCav') = mnNaN2zero(interp1(ff,sp,freq));
% Noise('sensor_RefCav') = zeros(size(freq));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unknown noises
% DAC noise
Noise('DAC') = 1e-7*2^16/20*ones(size(freq));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% angular motion of MC suspensions
% oplev spectrum
mnimportdata('N1130_1.txt',Freq,Abs,Phase,'Oplev_MCE_p','sp','diaggui')
[temp_ff1, temp_abs1] = mntrim(2.885,6.791,Freq('Oplev_MCE_p'),Abs('Oplev_MCE_p'));

Noise('Oplev_MCEp') = mnNaN2zero(interp1(temp_ff1,temp_abs1,freq));

mnimportdata('N1130_2.txt',Freq,Abs,Phase,'Oplev_MCE_y','sp','diaggui')
[temp_ff1, temp_abs1] = mntrim(0.1584,1.135,Freq('Oplev_MCE_y'),Abs('Oplev_MCE_y'));

Noise('Oplev_MCEy') = mnNaN2zero(interp1(temp_ff1,temp_abs1,freq));

Noise('Oplev_MCIp') = zeros(size(freq));
Noise('Oplev_MCOp') = zeros(size(freq));
Noise('Oplev_MCIy') = zeros(size(freq));
Noise('Oplev_MCOy') = zeros(size(freq));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if mnismember(unable_noise_list,'seismic_IMC')
    Noise('seismic_IMC') = zeros(size(freq));
end

if mnismember(unable_noise_list,'sensor_IMC')
    Noise('sensor_IMC') = zeros(size(freq));
end
if mnismember(unable_noise_list,'motion_RefCav')
    Noise('motion_RefCav') = zeros(size(freq));
end
if mnismember(unable_noise_list,'shot_IMC')
    Noise('motion_RefCav') = zeros(size(freq));
end
if mnismember(unable_noise_list,'A2f')
    Noise('Oplev_MCEp') = zeros(size(freq));
    Noise('Oplev_MCEy') = zeros(size(freq));
end
if mnismember(unable_noise_list,'act_AOM')
    Noise('act_AOM') = zeros(size(freq));
end