function LCGT_ASCnoise()
% noise calculation for LCGT ASC
% needs (1) LCGT_ASCfilterbank.m  for servo filter setting
%       (2) parambLCGT.m        for  acquiring suspension data file names
%       (3) plotbode.m          for bode plot
%       (4) plotspectrum.m      for plotting spectra 
%       (5) plotspectrumRMS.m   for plotting cumulative RMS of spectra
%       (6) plotA2DARM.m        for plotting angle to DARM coupling
% Author: Yuta Michimura

%% Clear everything, set save directory and aquire some data file names
%close all
clear all
cd('C:\Users\yuta\Lab\LCGT\IFOmodel\bLCGT_ASC')
modelName='bLCGT-BRSE-positive-yaw';
if strfind(modelName,'pitch')>0
    angleName='pitch';
elseif strfind(modelName,'yaw')>0
    angleName='yaw';
end
if strfind(modelName,'positive')>0
    gfactor='positive';
elseif strfind(modelName,'negative')>0
    gfactor='negative';
end
if strfind(modelName,'BRSE')>0
    rse='BRSE';
elseif strfind(modelName,'DRSE')>0
    rse='DRSE';
end
saveDir=['./results/',modelName,'-',date,'/'];
mkdir(saveDir);

sensitivitydatafile='./DARMsensitivity/bLCGTDARMsensitivity_VDRSE.txt';

if strfind(modelName,'bLCGT')>0
p=parambLCGT(rse,gfactor,angleName);  % doesn't matter 'positive' or 'negative'
colororder = [0.0   0.0   1.0
              1.0   0.0   0.0
              0.0   0.7   1.0
              1.0   0.7   0.0
              0.3   0.3   0.3
              0.0   1.0   0.0
              0.0   0.5   0.0
              0.5   0.5   0.0
              1.0   0.0   1.0
              0.7   0.0   0.5
              0.5   0.1   0.1];
elseif strfind(modelName,'iLCGT')>0
colororder = [0.0   0.0   1.0
              1.0   0.0   0.0
              0.0   0.7   1.0
              1.0   0.7   0.0
              0.3   0.3   0.3
              0.0   1.0   0.0
              0.0   0.5   0.0];
end

%% Load matrices and vectors (output of LCGT_ASC.m)
loadDir=['./results/',modelName,'-',date,'/'];
load([loadDir,'mOpt']);          % matrix of IFO response [W/rad] (Nprb x Ndrv x Nfreq)
load([loadDir,'mSens']);         % angular sensing matrix for calculating input matrix (Nprb x Ndof) 
load([loadDir,'mDOFtoMIRROR']);  % matrix that transforms DOF base to mirror base (Ndrv x Ndof)
load([loadDir,'mMech']);         % matrix of modification of actuator TF under radiation pressure (Ndrv x Ndrv x Nfreq)
load([loadDir,'mBSM']);          % matrix TF from each mirror angular motion to beam spot motion on each mirror [m/rad] (Ndrv x Ndrv x Nfreq)
load([loadDir,'freq']);          % vector of frequency (1 x Nfreq)
load([loadDir,'vProbePower']);   % vector of powers at each probes (Nprb x 1)
load([loadDir,'vMirrorNames']);  % vector of mirror names (1 x Ndrv)
load([loadDir,'vDOFNames']);     % vector of DOF names (1 x Ndof)
load([loadDir,'vProbeNames']);   % vector of selected probe names (1 x Nprb)

% Number of things
Nfreq=length(freq);
Ndrv=length(vMirrorNames);
Ndof=length(vDOFNames);
Nprb=length(vProbeNames);
    
% DOFs not to control by WFS
if strfind(modelName,'bLCGT')>0
    vUncontrolledDOFNames={'SR2'};
elseif strfind(modelName,'iLCGT')>0
    vUncontrolledDOFNames={'PR2'};
end
Nuncontrolled=length(vUncontrolledDOFNames);
iUncontrolledDOF=length(Nuncontrolled);  % index of uncontrolled DOF
for kk=1:Ndof
    for jj=1:Nuncontrolled
        if strfind(vDOFNames{kk},vUncontrolledDOFNames{jj})>0
            iUncontrolledDOF(jj)=kk;
        end
    end
end
 
% Constants
hP=6.6261e-34;   % Planck's constant
c=299792458;     % speed of light
lambda=1.064e-6; % wavelength of laser
nu=c/lambda;     % frequency of laser

%% Calculate input matrix
mSens_ignoredsmallsignal=mSens;

for ii=1:Nprb
    maxsignal=max(abs(mSens(ii,:)));
    for kk=1:Ndof
        if abs(mSens(ii,kk)) < maxsignal*0.10
            mSens_ignoredsmallsignal(ii,kk)=0;   % ignore small signal (smaller than 10% of maximum signal)
        end
    end
end

for kk=1:Nuncontrolled
    mSens_ignoredsmallsignal(:,iUncontrolledDOF(kk))=0.*mSens_ignoredsmallsignal(:,iUncontrolledDOF(kk));
end

mInput=pinv(mSens_ignoredsmallsignal);
mat=mInput*mSens;
for kk=1:Ndof
    if mat(kk,kk) ~= 0
        mInput(kk,:)=1/mat(kk,kk)*mInput(kk,:);   % gain balancing for convenience
    end
end
mInput
save([saveDir,'mInput'],'mInput');  % matrix that reconstruct DOF erorr signals from probe signals (Ndof x Nprb)

%% Actuator TFs
% torque on RM to mirror angle
if strfind(modelName,'bLCGT')>0
    H_data=importdata(p.TM_Hact_TFdata);
    H=zeros(Nfreq,3);
    H(:,1:2)=change_sampling([H_data(:,1),H_data(:,2)],freq);
    H(:,1:2:3)=change_sampling([H_data(:,1),H_data(:,3)],freq);
    HTMact=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
    H_data=importdata(p.BS_Hact_TFdata);
    H=zeros(Nfreq,3);
    H(:,1:2)=change_sampling([H_data(:,1),H_data(:,2)],freq);
    H(:,1:2:3)=change_sampling([H_data(:,1),H_data(:,3)],freq);
    HBSact=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
    H_data=importdata(p.RC_Hact_TFdata);
    H=zeros(Nfreq,3);
    H(:,1:2)=change_sampling([H_data(:,1),H_data(:,2)],freq);
    H(:,1:2:3)=change_sampling([H_data(:,1),H_data(:,3)],freq);
    HRCact=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
elseif strfind(modelName,'iLCGT')>0
    H_data=importdata(['../iLCGT_ASC/seism_angle_110513/TypeA/typeA_tf_wRM_',angleName,'_110513.txt']);
    H_data(1,1)=round(H_data(1,1)*1000)/1000;
    H_data(2001,1)=round(H_data(2001,1)*1000)/1000;
    H=zeros(Nfreq,3);
    H(:,1:2)=change_sampling([H_data(:,1),H_data(:,2)],freq);
    H(:,1:2:3)=change_sampling([H_data(:,1),H_data(:,3)],freq);
    HTMact=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
    H_data=importdata(['../iLCGT_ASC/seism_angle_110513/TypeB_BS/BS_tf_wRM_',angleName,'_110513.txt']);
    H_data(1,1)=round(H_data(1,1)*1000)/1000;
    H_data(2001,1)=round(H_data(2001,1)*1000)/1000;
    H=zeros(Nfreq,3);
    H(:,1:2)=change_sampling([H_data(:,1),H_data(:,2)],freq);
    H(:,1:2:3)=change_sampling([H_data(:,1),H_data(:,3)],freq);
    HBSact=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
    H_data=importdata(['../iLCGT_ASC/seism_angle_110513/TypeB_RM/RM_tf_wRM_',angleName,'_110513.txt']);
    H_data(1,1)=round(H_data(1,1)*1000)/1000;
    H_data(2001,1)=round(H_data(2001,1)*1000)/1000;
    H=zeros(Nfreq,3);
    H(:,1:2)=change_sampling([H_data(:,1),H_data(:,2)],freq);
    H(:,1:2:3)=change_sampling([H_data(:,1),H_data(:,3)],freq);
    HRCact=squeeze(H(:,2).*exp(i.*H(:,3)/180*pi));
end

mAct=zeros(Ndrv,Ndrv,Nfreq);
for kk=1:Ndrv
    if strfind(vMirrorNames{kk},'TM')>0
        mAct(kk,kk,:)=HTMact;
    elseif strfind(vMirrorNames{kk},'BS')>0
        mAct(kk,kk,:)=HBSact;
    elseif strfind(vMirrorNames{kk},'PR')>0
        mAct(kk,kk,:)=HRCact;
    elseif strfind(vMirrorNames{kk},'SR')>0
        mAct(kk,kk,:)=HRCact;
    end
end
figure();
plotbode(freq,[HTMact,HBSact,HRCact]);
subplot(2,1,1);
title(['Actuator TF (',angleName,', torque from RM to angle of TM)'])
ylabel('Abs [rad/(Nm)]');
legend('TMs','BS','RC mirrors');

saveas(gcf,[saveDir,'ActuatorTF.png'])

%% designing filters (skip)
    ZQ=[1 0
        1 0
        0.15 10
        0.87 10];
    PQ=[0.01 0
        0.02 10
        0.34 10
        10 1
        14 2
        20 5];
    F=filterdesign(freq,ZQ,PQ,1e12);
%    F=freq'./freq';

dofNo=1;
G_DOF=zeros(Ndof,Ndof,Nfreq);
for kk=1:Nfreq
    G_DOF(:,:,kk)=mInput*mOpt(:,:,kk)*mAct(:,:,kk)*mDOFtoMIRROR;
    G_DOF(dofNo,dofNo,kk)=G_DOF(dofNo,dofNo,kk)*F(kk);
end

figure();
plotbode(freq,[squeeze(G_DOF(dofNo,dofNo,:)),F]);
subplot(2,1,1);
title(['Openloop for ',vDOFNames{dofNo}]);

%% Set control filter matrix
mFilt=zeros(Ndof,Ndof,Nfreq);
filt=LCGT_ASCfilterbank(freq,modelName);

for kk=1:Ndof
    if strfind(vDOFNames{kk},'CS')>0
        mFilt(kk,kk,:)=filt.FS;
    elseif strfind(vDOFNames{kk},'DS')>0
        mFilt(kk,kk,:)=filt.FS;
    elseif strfind(vDOFNames{kk},'CH')>0
        mFilt(kk,kk,:)=filt.FH;
    elseif strfind(vDOFNames{kk},'DH')>0
        mFilt(kk,kk,:)=filt.FH;
    elseif strfind(vDOFNames{kk},'BS')>0
        mFilt(kk,kk,:)=filt.FBS;
    elseif strfind(vDOFNames{kk},'PR')>0
        mFilt(kk,kk,:)=filt.FPRC;
    elseif strfind(vDOFNames{kk},'SR')>0
        mFilt(kk,kk,:)=filt.FSRC;
    end
end

for kk=1:Nuncontrolled
    mFilt(iUncontrolledDOF(kk),iUncontrolledDOF(kk),:)=filt.FPRC*0;
end

figure();
plotbode(freq,[filt.FS,filt.FH,filt.FBS,filt.FPRC,filt.FSRC]);
subplot(2,1,1);
title('Servo Filters');
ylim([1e-1,1e7]);
set(gca,'YTick',logspace(-1,7,9));
legend('SOFT','HARD','BS','PRC mirrors','SRC mirrors');

saveas(gcf,[saveDir,'ServoFilters.png'])

%% Plot openloop TF for each DOF
G_DOF=zeros(Ndof,Ndof,Nfreq);
for kk=1:Nfreq
    G_DOF(:,:,kk)=mFilt(:,:,kk)*mInput*mOpt(:,:,kk)*mAct(:,:,kk)*mDOFtoMIRROR;
end

G_DOFdiag=[];  % diagonal element of G_DOF
for kk=1:Ndof
    G_DOFdiag=[G_DOFdiag,squeeze(G_DOF(kk,kk,:))];
end

figure();
set(gcf,'DefaultAxesColorOrder',colororder)

%set(gcf,'DefaultLineLineWidth',3);
%set(gcf,'defaultAxesFontSize',16);
%set(gcf,'defaultTextFontSize',16);

plotbode(freq,G_DOFdiag);
subplot(2,1,1);
title(['Openloop TF for each DOF (',angleName,')']);
ylim([1e-5,1e5]);
set(gca,'YTick',logspace(-5,5,11));
legend(vDOFNames);

saveas(gcf,[saveDir,'openloopTF_DOF.png'])
saveas(gcf,[saveDir,'openloopTF_DOF.eps'],'epsc')

%% Plot openloop TF for each mirrors
G_mirror=zeros(Ndrv,Ndrv,Nfreq);
for kk=1:Nfreq
    G_mirror(:,:,kk)=mAct(:,:,kk)*mDOFtoMIRROR*mFilt(:,:,kk)*mInput*mOpt(:,:,kk);
end

G_mirrordiag=[];  % diagonal element of G_mirror
for kk=1:Ndrv
    G_mirrordiag=[G_mirrordiag,squeeze(G_mirror(kk,kk,:))];
end

figure();
set(gcf,'DefaultAxesColorOrder',colororder)
plotbode(freq,G_mirrordiag);
subplot(2,1,1);
title(['Openloop TF for each mirror (',angleName,')']);
ylim([1e-5,1e5]);
set(gca,'YTick',logspace(-5,5,11));
legend(vMirrorNames);

saveas(gcf,[saveDir,'openloopTF_mirror.png'])

%% Seismic noise

if strfind(modelName,'bLCGT')>0
    TypeA=importdata(p.TM_seism_TFdata);
    TypeB_BS=importdata(p.BS_seism_TFdata);
    TypeB_RC=importdata(p.RC_seism_TFdata);
elseif strfind(modelName,'iLCGT')>0
    TypeA=importdata(['../iLCGT_ASC/seism_angle_110513/TypeA/typeA_seism_',angleName,'_110513.txt']);
    TypeB_BS=importdata(['../iLCGT_ASC/seism_angle_110513/TypeB_BS/BS_seism_',angleName,'_110513.txt']);
    TypeB_RC=importdata(['../iLCGT_ASC/seism_angle_110513/TypeB_RM/RM_seism_',angleName,'_110513.txt']);
    TypeA(1,1)=round(TypeA(1,1)*1000)/1000;
    TypeA(2001,1)=round(TypeA(2001,1)*1000)/1000;
    TypeB_BS(1,1)=round(TypeB_BS(1,1)*1000)/1000;
    TypeB_BS(2001,1)=round(TypeB_BS(2001,1)*1000)/1000;
    TypeB_RC(1,1)=round(TypeB_RC(1,1)*1000)/1000;
    TypeB_RC(2001,1)=round(TypeB_RC(2001,1)*1000)/1000;
end

SeismicTM=change_sampling(TypeA,freq);
SeismicBS=change_sampling(TypeB_BS,freq);
SeismicRC=change_sampling(TypeB_RC,freq);

vSeismic=zeros(Ndrv,Nfreq);
for kk=1:Ndrv
    if strfind(vMirrorNames{kk},'TM')>0
        vSeismic(kk,:)=SeismicTM(:,2);
    elseif strfind(vMirrorNames{kk},'BS')>0
        vSeismic(kk,:)=SeismicBS(:,2);
    elseif strfind(vMirrorNames{kk},'PR')>0
        vSeismic(kk,:)=SeismicRC(:,2);
    elseif strfind(vMirrorNames{kk},'SR')>0
        vSeismic(kk,:)=SeismicRC(:,2);
    end
end

seismtitle=['Input Seismic Noise (',angleName,')'];
figure();
set(gcf,'DefaultAxesColorOrder',[0 0 1;0 0.5 0;1 0 0])
plotspectrumRMS(freq,[SeismicTM(:,2)';SeismicBS(:,2)';SeismicRC(:,2)']);
title(seismtitle);
ylabel('Angular noise [rad/rtHz], RMS [rad]');
legend('Type-A (TMs)','Type-B (BS)','Type-B (RC mirrors)');
ylim([1e-18,1e-4]);
set(gca,'YTick',logspace(-18,-4,15));
saveas(gcf,[saveDir,'InputSeismicNoise.png'])
saveas(gcf,[saveDir,'InputSeismicNoise.eps'],'epsc')

%% Residual angular motion 
vShot=sqrt(2*hP*nu*vProbePower);  % shot noise of each probe [W/rtHz] (when quantum efficiency = 1)
theta_residual=zeros(Ndrv,Nfreq);
for kk=1:Nfreq
    theta_residual(:,kk)=abs(mMech(:,:,kk))*sqrt((abs(inv(eye(Ndrv)+G_mirror(:,:,kk)))*vSeismic(:,kk)).^2+(abs(inv(eye(Ndrv)+G_mirror(:,:,kk))*G_mirror(:,:,kk)*pinv(mOpt(:,:,kk)))*vShot).^2);
end
if strfind(modelName,'pitch')>0
    restitle='Residual Pitch Motion (with cumulative RMS)';
elseif strfind(modelName,'yaw')>0
    restitle='Residual Yaw Motion (with cumulative RMS)';
end

figure();
set(gcf,'DefaultAxesColorOrder',colororder)

%set(gcf,'DefaultLineLineWidth',3);
%set(gcf,'defaultAxesFontSize',16);
%set(gcf,'defaultTextFontSize',16);

RMStheta_residual=plotspectrumRMS(freq,theta_residual);
title(restitle);
ylabel('Angular noise [rad/rtHz], RMS [rad]');
legend(vMirrorNames);
ylim([1e-18,1e-4]);
set(gca,'YTick',logspace(-18,-4,15));

saveas(gcf,[saveDir,'ResidualAngularMotionRMS.png'])
saveas(gcf,[saveDir,'ResidualAngularMotionRMS.eps'],'epsc')

%% Beam spot motion
beamspotmotion=zeros(Ndrv,Nfreq);
for kk=1:Nfreq
    beamspotmotion(:,kk)=abs(mBSM(:,:,kk))*inv(abs(mMech(:,:,kk)))*theta_residual(:,kk);
end
if strfind(modelName,'pitch')>0
    bsmtitle='Beam Spot Pitch Motion (with cumulative RMS)';
elseif strfind(modelName,'yaw')>0
    bsmtitle='Beam Spot Yaw Motion (with cumulative RMS)';
end

figure();
set(gcf,'DefaultAxesColorOrder',colororder)

%set(gcf,'DefaultLineLineWidth',3);
%set(gcf,'defaultAxesFontSize',16);
%set(gcf,'defaultTextFontSize',16);

RMSbeamspotmotion=plotspectrumRMS(freq,beamspotmotion);
title(bsmtitle);
ylabel('Displacement Noise [m/rtHz], RMS [m]');
legend(vMirrorNames);
ylim([1e-10,1e-0]);
set(gca,'YTick',logspace(-10,-0,11));

saveas(gcf,[saveDir,'BeamSpotMotionRMS.png'])
saveas(gcf,[saveDir,'BeamSpotMotionRMS.eps'],'epsc')

%% Coupling factor to DARM
% see LIGO-T0900511 section 5.2
Finesse=1546;  % Finesse of the arm cavity
coupling=zeros(1,Ndrv);
for kk=1:Ndrv
    if strfind(vMirrorNames{kk},'TM')>0
        coupling(kk)=1;
    elseif strfind(vMirrorNames{kk},'BS')>0
        coupling(kk)=pi/(2*Finesse);
    elseif strfind(vMirrorNames{kk},'PR')>0
        coupling(kk)=1/100;
    elseif strfind(vMirrorNames{kk},'SR')>0
        coupling(kk)=1/100;
    end
end

%% DARM sensitivity
% for now, use bLCGT sensitivity
DARMsensitivity_data=importdata(sensitivitydatafile);
DARMsensitivity=change_sampling(DARMsensitivity_data,freq);

figure();
plotspectrum(freq,DARMsensitivity(:,2));
ylabel('Displacement noise [m/rtHz]');
title('bLCGT Sensitivity');

%% Angular noise coupling to DARM
A2DARM=zeros(Ndrv,Nfreq);
for kk=1:Ndrv
    A2DARM(kk,:)=coupling(kk)*(RMSbeamspotmotion(kk)*theta_residual(kk,:)+RMStheta_residual(kk)*beamspotmotion(kk,:));  % see LIGO-T0900511 eq(5.2)
end
totalA2DARM=zeros(1,Nfreq);
for kk=1:Nfreq
    totalA2DARM(kk)=sqrt(sum(A2DARM(:,kk).^2));
end
totalA2DARM=[freq',totalA2DARM'];
dlmwrite([saveDir,'totalA2DARM.txt'],totalA2DARM,'delimiter','\t','precision',6);

A2DARM=[A2DARM;DARMsensitivity(:,2)'];
figure();
set(gcf,'DefaultAxesColorOrder',colororder)

%set(gcf,'DefaultLineLineWidth',3);
%set(gcf,'defaultAxesFontSize',16);
%set(gcf,'defaultTextFontSize',16);

plotA2DARM(freq,A2DARM);
legend(vMirrorNames,1);
ylim([1e-21,1e-8]);
set(gca,'YTick',logspace(-21,-8,14));
a2darmtitle=['Angular Noise Coupling to DARM (',angleName,')'];
title(a2darmtitle);

saveas(gcf,[saveDir,'A2DARM.png'])
saveas(gcf,[saveDir,'A2DARM.eps'],'epsc')

%% Noise contribution from each probe
vProbeNoise=zeros(Nprb,Nfreq);
for jj=1:Nprb
    vShot=zeros(Nprb,1);
    vShot(jj)=sqrt(2*hP*nu*vProbePower(jj));
    theta_residual=zeros(Ndrv,Nfreq);
    for kk=1:Nfreq
        theta_residual(:,kk)=abs(mMech(:,:,kk))*sqrt((abs(inv(eye(Ndrv)+G_mirror(:,:,kk)))*vSeismic(:,kk)).^2+(abs(inv(eye(Ndrv)+G_mirror(:,:,kk))*G_mirror(:,:,kk)*pinv(mOpt(:,:,kk)))*vShot).^2);
    end
    RMStheta_residual=plotspectrumRMS(freq,theta_residual);
    beamspotmotion=zeros(Ndrv,Nfreq);
    for kk=1:Nfreq
        beamspotmotion(:,kk)=abs(mBSM(:,:,kk))*inv(abs(mMech(:,:,kk)))*theta_residual(:,kk);
    end
    RMSbeamspotmotion=plotspectrumRMS(freq,beamspotmotion);
    A2DARM=zeros(Ndrv,Nfreq);
    for kk=1:Ndrv
        A2DARM(kk,:)=coupling(kk)*(RMSbeamspotmotion(kk)*theta_residual(kk,:)+RMStheta_residual(kk)*beamspotmotion(kk,:));  % see LIGO-T0900511 eq(5.2)
    end
    for kk=1:Nfreq
        vProbeNoise(jj,kk)=sqrt(sum(A2DARM(:,kk).^2));
    end
end
close();
vProbeNoise=[vProbeNoise;DARMsensitivity(:,2)'];
figure();
set(gcf,'DefaultAxesColorOrder',colororder)
plotA2DARM(freq,vProbeNoise);
legend([vProbeNames,'bLCGT sensitivity'],'interpreter','none');
ylim([1e-21,1e-8]);
set(gca,'YTick',logspace(-21,-8,14));
title(['Noise Contribution of each Probe (',angleName,')']);

saveas(gcf,[saveDir,'ProbeNoise.png'])

%% Noise contribution from each Mirror
vMirrorNoise=zeros(Ndrv,Nfreq);
vShot=sqrt(2*hP*nu*vProbePower);  % shot noise of each probe [W/rtHz]
for jj=1:Ndrv
    theta_residual=zeros(Ndrv,Nfreq);
    vSeismic0=zeros(Ndrv,Nfreq);
    vSeismic0(jj,:)=vSeismic(jj,:);
    for kk=1:Nfreq
        theta_residual(:,kk)=abs(mMech(:,:,kk))*sqrt((abs(inv(eye(Ndrv)+G_mirror(:,:,kk)))*vSeismic0(:,kk)).^2+(abs(inv(eye(Ndrv)+G_mirror(:,:,kk))*G_mirror(:,:,kk)*pinv(mOpt(:,:,kk)))*vShot).^2);
    end
    RMStheta_residual=plotspectrumRMS(freq,theta_residual);
    beamspotmotion=zeros(Ndrv,Nfreq);
    for kk=1:Nfreq
        beamspotmotion(:,kk)=abs(mBSM(:,:,kk))*inv(abs(mMech(:,:,kk)))*theta_residual(:,kk);
    end
    RMSbeamspotmotion=plotspectrumRMS(freq,beamspotmotion);
    A2DARM=zeros(Ndrv,Nfreq);
    for kk=1:Ndrv
        A2DARM(kk,:)=coupling(kk)*(RMSbeamspotmotion(kk)*theta_residual(kk,:)+RMStheta_residual(kk)*beamspotmotion(kk,:));  % see LIGO-T0900511 eq(5.2)
    end
    for kk=1:Nfreq
        vMirrorNoise(jj,kk)=sqrt(sum(A2DARM(:,kk).^2));
    end
end
close();
vMirrorNoise=[vMirrorNoise;DARMsensitivity(:,2)'];
figure();
set(gcf,'DefaultAxesColorOrder',colororder)
plotA2DARM(freq,vMirrorNoise);
legend([vMirrorNames,'bLCGT sensitivity'],'interpreter','none');
ylim([1e-21,1e-8]);
set(gca,'YTick',logspace(-21,-8,14));
title(['Noise Contribution of each Mirror (',angleName,')']);

saveas(gcf,[saveDir,'MirrorNoise.png'])

%% Angular noise coupling to DARM without control loop
theta_residual=zeros(Ndrv,Nfreq);
for kk=1:Nfreq
    theta_residual(:,kk)=abs(mMech(:,:,kk))*vSeismic(:,kk);
end
figure();
set(gcf,'DefaultAxesColorOrder',colororder)

%set(gcf,'DefaultLineLineWidth',3);
%set(gcf,'defaultAxesFontSize',16);
%set(gcf,'defaultTextFontSize',16);

RMStheta_residual=plotspectrumRMS(freq,theta_residual);
title(restitle);
ylabel('Angular noise [rad/rtHz], RMS [rad]');
legend(vMirrorNames);
ylim([1e-18,1e-4]);
set(gca,'YTick',logspace(-18,-4,15));

saveas(gcf,[saveDir,'ResidualAngularMotionRMSnocontrol.png'])

beamspotmotion=zeros(Ndrv,Nfreq);
for kk=1:Nfreq
    beamspotmotion(:,kk)=abs(mBSM(:,:,kk))*inv(abs(mMech(:,:,kk)))*theta_residual(:,kk);
end
figure();
set(gcf,'DefaultAxesColorOrder',colororder)
RMSbeamspotmotion=plotspectrumRMS(freq,beamspotmotion);
title(bsmtitle);
ylabel('Displacement Noise [m/rtHz], RMS [m]');
legend(vMirrorNames);
ylim([1e-10,1e-0]);
set(gca,'YTick',logspace(-10,-0,11));

saveas(gcf,[saveDir,'BeamSpotMotionRMSnocontrol.png'])

A2DARM=zeros(Ndrv,Nfreq);
for kk=1:Ndrv
    A2DARM(kk,:)=coupling(kk)*(RMSbeamspotmotion(kk)*theta_residual(kk,:)+RMStheta_residual(kk)*beamspotmotion(kk,:));  % see LIGO-T0900511 eq(5.2)
end
A2DARM=[A2DARM;DARMsensitivity(:,2)'];
figure();
set(gcf,'DefaultAxesColorOrder',colororder)
plotA2DARM(freq,A2DARM);
legend(vMirrorNames,'bLCGT sensitivity');
ylim([1e-21,1e-8]);
set(gca,'YTick',logspace(-21,-8,14));
title([a2darmtitle,' (without ASC)']);

saveas(gcf,[saveDir,'A2DARMnocontrol.png'])