function Sensitivity2

  % create the model
  clear
  load lcgt_optmat_tuned2.mat
%  opt = lcgt; % opt is changed in extractTF2 so we must not reload lcgt.
    
    % get some drive indexes
  nEX = getDriveIndex(opt, 'EMX');
  nEY = getDriveIndex(opt, 'EMY');

  % compute the DC signals and TFs on resonances
  pos = zeros(opt.Ndrive, 1);

  clear fDC sigDC1 sigAC1 mechAC1 noiseAC1 pos
  
  % compute the DC signals and TFs on resonances
  pos = zeros(opt.Ndrive, 1);
  [fDC, sigDC1, sigAC1, mechAC1, noiseAC1] = tickle(opt, pos, f);

  % make a response plot
  h11 = (0.5*getTF(sigAC1, nCARM_P, nEX) + 0.5*getTF(sigAC1, nCARM_P, nEY))*cos(nCARM_phi) + (0.5*getTF(sigAC1, nCARM_Q, nEX) + 0.5*getTF(sigAC1, nCARM_Q, nEY))*sin(nCARM_phi);
  h12 = (0.5*getTF(sigAC1, nCARM_P, nEX) - 0.5*getTF(sigAC1, nCARM_P, nEY))*cos(nCARM_phi) + (0.5*getTF(sigAC1, nCARM_Q, nEX) - 0.5*getTF(sigAC1, nCARM_Q, nEY))*sin(nCARM_phi);
  h13 = -getTF(sigAC1, nCARM_P, nPRM)*cos(nCARM_phi) - getTF(sigAC1, nCARM_Q, nPRM)*sin(nCARM_phi);
  h14 = ((-1*getTF(sigAC1hp, nCARM_P, nPRM) + 1*getTF(sigAC1hp, nCARM_P, nSRM) + sqrt(2)*getTF(sigAC1hp, nCARM_P, nBS)))*sin(nCARM_phi) + ((-1*getTF(sigAC1hp, nCARM_Q, nPRM) + 1*getTF(sigAC1hp, nCARM_Q, nSRM) + sqrt(2)*getTF(sigAC1hp, nCARM_Q, nBS)))*cos(nCARM_phi);
  h15 = -getTF(sigAC1, nCARM_P, nSRM)*cos(nCARM_phi) - getTF(sigAC1, nCARM_Q, nSRM)*sin(nCARM_phi);
  h21 = 0.5*getTF(sigAC1, nDARM, nEX) + 0.5*getTF(sigAC1, nDARM, nEY);
  h22 = 0.5*getTF(sigAC1, nDARM, nEX) - 0.5*getTF(sigAC1, nDARM, nEY);
  h23 = -getTF(sigAC1, nDARM, nPRM);
  h24 = ((-1*getTF(sigAC1hp, nDARM, nPRM) + 1*getTF(sigAC1hp, nDARM, nSRM) + sqrt(2)*getTF(sigAC1hp, nDARM, nBS)));
  h25 = -getTF(sigAC1, nDARM, nSRM);
  h31 = (0.5*getTF(sigAC1, nPRCL_P, nEX) + 0.5*getTF(sigAC1, nPRCL_P, nEY))*cos(nPRCL_phi) + (0.5*getTF(sigAC1, nPRCL_Q, nEX) + 0.5*getTF(sigAC1, nPRCL_Q, nEY))*sin(nPRCL_phi);
  h32 = (0.5*getTF(sigAC1, nPRCL_P, nEX) - 0.5*getTF(sigAC1, nPRCL_P, nEY))*cos(nPRCL_phi) + (0.5*getTF(sigAC1, nPRCL_Q, nEX) - 0.5*getTF(sigAC1, nPRCL_Q, nEY))*sin(nPRCL_phi);
  h33 = -getTF(sigAC1, nPRCL_P, nPRM)*cos(nPRCL_phi) - getTF(sigAC1, nPRCL_Q, nPRM)*sin(nPRCL_phi);
  h34 = ((-1*getTF(sigAC1hp, nPRCL_P, nPRM) + 1*getTF(sigAC1hp, nPRCL_P, nSRM) + sqrt(2)*getTF(sigAC1hp, nPRCL_P, nBS)))*sin(nPRCL_phi) + ((-1*getTF(sigAC1hp, nPRCL_Q, nPRM) + 1*getTF(sigAC1hp, nPRCL_Q, nSRM) + sqrt(2)*getTF(sigAC1hp, nPRCL_Q, nBS)))*cos(nPRCL_phi);
  h35 = -getTF(sigAC1, nPRCL_P, nSRM)*cos(nPRCL_phi) - getTF(sigAC1, nPRCL_Q, nSRM)*sin(nPRCL_phi);
  h41 = (0.5*getTF(sigAC1, nMICH_P, nEX) + 0.5*getTF(sigAC1, nMICH_P, nEY))*cos(nMICH_phi) + (0.5*getTF(sigAC1, nMICH_Q, nEX) + 0.5*getTF(sigAC1, nMICH_Q, nEY))*sin(nMICH_phi);
  h42 = (0.5*getTF(sigAC1, nMICH_P, nEX) - 0.5*getTF(sigAC1, nMICH_P, nEY))*cos(nMICH_phi) + (0.5*getTF(sigAC1, nMICH_Q, nEX) - 0.5*getTF(sigAC1, nMICH_Q, nEY))*sin(nMICH_phi);
  h43 = -getTF(sigAC1, nMICH_P, nPRM)*cos(nMICH_phi) - getTF(sigAC1, nMICH_Q, nPRM)*sin(nMICH_phi);
  h44 = ((-1*getTF(sigAC1hp, nMICH_P, nPRM) + 1*getTF(sigAC1hp, nMICH_P, nSRM) + sqrt(2)*getTF(sigAC1hp, nMICH_P, nBS)))*sin(nMICH_phi) + ((-1*getTF(sigAC1hp, nMICH_Q, nPRM) + 1*getTF(sigAC1hp, nMICH_Q, nSRM) + sqrt(2)*getTF(sigAC1hp, nMICH_Q, nBS)))*cos(nMICH_phi);
  h45 = -getTF(sigAC1, nMICH_P, nSRM)*cos(nMICH_phi) - getTF(sigAC1, nMICH_Q, nSRM)*sin(nMICH_phi);
  h51 = (0.5*getTF(sigAC1, nSRCL_P, nEX) + 0.5*getTF(sigAC1, nSRCL_P, nEY))*cos(nSRCL_phi) + (0.5*getTF(sigAC1, nSRCL_Q, nEX) + 0.5*getTF(sigAC1, nSRCL_Q, nEY))*sin(nSRCL_phi);
  h52 = (0.5*getTF(sigAC1, nSRCL_P, nEX) - 0.5*getTF(sigAC1, nSRCL_P, nEY))*cos(nSRCL_phi) + (0.5*getTF(sigAC1, nSRCL_Q, nEX) - 0.5*getTF(sigAC1, nSRCL_Q, nEY))*sin(nSRCL_phi);
  h53 = -getTF(sigAC1, nSRCL_P, nPRM)*cos(nSRCL_phi) - getTF(sigAC1, nSRCL_Q, nPRM)*sin(nSRCL_phi);
  h54 = ((-1*getTF(sigAC1hp, nSRCL_P, nPRM) + 1*getTF(sigAC1hp, nSRCL_P, nSRM) + sqrt(2)*getTF(sigAC1hp, nSRCL_P, nBS)))*sin(nSRCL_phi) + ((-1*getTF(sigAC1hp, nSRCL_Q, nPRM) + 1*getTF(sigAC1hp, nSRCL_Q, nSRM) + sqrt(2)*getTF(sigAC1hp, nSRCL_Q, nBS)))*cos(nSRCL_phi);
  h55 = -getTF(sigAC1, nSRCL_P, nSRM)*cos(nSRCL_phi) - getTF(sigAC1, nSRCL_Q, nSRM)*sin(nSRCL_phi);
  
  % make a noise plot
  n1 = noiseAC1(nCARM_P, :)'*cos(nCARM_phi) + noiseAC1(nCARM_Q, :)'*sin(nCARM_phi);
  n2 = noiseAC1(nDARM, :)';
  n3 = noiseAC1(nPRCL_P, :)'*cos(nPRCL_phi) + noiseAC1(nPRCL_Q, :)'*sin(nPRCL_phi);
  n4 = noiseAC1(nMICH_P, :)'*cos(nMICH_phi) + noiseAC1(nMICH_Q, :)'*sin(nMICH_phi);
  n5 = noiseAC1(nSRCL_P, :)'*cos(nSRCL_phi) + noiseAC1(nSRCL_Q, :)'*sin(nSRCL_phi);
  
clear HH GG nn
  HH=zeros(length(f),5,5);
  GG=zeros(length(f),5,5);
  nn=zeros(length(f),5);
  
  HH(:,1,1)=h11(:);
  HH(:,1,2)=h12(:);
  HH(:,1,3)=h13(:);
  HH(:,1,4)=h14(:);
  HH(:,1,5)=h15(:);
  HH(:,2,1)=h21(:);
  HH(:,2,2)=h22(:);
  HH(:,2,3)=h23(:);
  HH(:,2,4)=h24(:);
  HH(:,2,5)=h25(:);
  HH(:,3,1)=h31(:);
  HH(:,3,2)=h32(:);
  HH(:,3,3)=h33(:);
  HH(:,3,4)=h34(:);
  HH(:,3,5)=h35(:);
  HH(:,4,1)=h41(:);
  HH(:,4,2)=h42(:);
  HH(:,4,3)=h43(:);
  HH(:,4,4)=h44(:);
  HH(:,4,5)=h45(:);
  HH(:,5,1)=h51(:);
  HH(:,5,2)=h52(:);
  HH(:,5,3)=h53(:);
  HH(:,5,4)=h54(:);
  HH(:,5,5)=h55(:);

  UGF1 = 10000;
  UGF2 = 100;
  UGF3 = 20;
  
  GG(:,1,1)=1 ./ HH(:,1,1) * UGF1^1.5./f(:).^1.5;
  GG(:,2,2)=1 ./ HH(:,2,2) * UGF2^1.5./f(:).^1.5;
  GG(:,3,3)=1 ./ HH(:,3,3) * UGF3^1.5./f(:).^1.5;
  GG(:,4,4)=1 ./ HH(:,4,4) * UGF3^1.5./f(:).^1.5;
  GG(:,5,5)=1 ./ HH(:,5,5) * UGF3^1.5./f(:).^1.5;
  
  nn(:,1)=n1;
  nn(:,2)=n2;
  nn(:,3)=n3;
  nn(:,4)=n4;
  nn(:,5)=n5;
  
  load hclassic.dat
  disp_totalh = hclassic(:, 2);
  disp_total = interp1(hclassic(:, 1), disp_totalh, f);
  load hclassicPRCL.dat
  disp_totalhPRCL = hclassicPRCL(:, 2);
  disp_totalPRCL = interp1(hclassicPRCL(:, 1), disp_totalhPRCL, f);
  load hclassicSRCL.dat
  disp_totalhSRCL = hclassicSRCL(:, 2);
  disp_totalSRCL = interp1(hclassicSRCL(:, 1), disp_totalhSRCL, f);
  load hclassicMICH.dat
  disp_totalhMICH = hclassicMICH(:, 2);
  disp_totalMICH = interp1(hclassicMICH(:, 1), disp_totalhMICH, f);

  xi = zeros(length(f),5);
  Id = eye(5);
xi(:,1) = disp_total;
xi(:,3) = disp_totalPRCL;
xi(:,4) = disp_totalMICH;
xi(:,5) = disp_totalSRCL;

%  Aij = zeros(5);
  Dij = zeros(5);
%  Greal = zeros(length(f),5,5);
  yt = zeros(length(f),5,5);
  yi = zeros(length(f),5);
  hGWt = zeros(length(f),5,5);
  hGW = zeros(length(f),5);
  
    load fname Greal;

    for ff = 1 : length(f)
%    Aij(:, :) = GG(ff, :, :);    
    Dij(:, :) = HH(ff, :, :);  % Dij at one particular frequency
%    Gopt = inv(Dij) * Aij;
%    Gopt = Dij \ Aij;  % faster way to execute inverse
%    for a = 1 : 5
%      Greal(ff,a,a) = Gopt(a,a);
%    end

%%% Feedforward

    Greal2d(:,:) = Greal(ff,:,:);
    Mij = -Dij * Greal2d;

%    Bij = inv(Id - Mij);       % Bij = (1 - Mij)^-1
%    BD = Bij * Dij;
    Bij = (Id - Mij) \ eye(5);       % Bij = (1 - Mij)^-1
    BD = (Id - Mij) \ Dij;   % faster way to execute inverse
 
    Dcomp=ones(1,5)'; % temporary
    
    for a = 1 : 5
      for b = 1 : 5
        yt(ff, a, b) = sqrt((BD(a, b) * xi(ff, b))^2 + (Bij(a, b) * (nn(ff, b) * Dcomp(b)))^2); % each noise contribution to y_i
        hGWt(ff, a, b) = yt(ff, a, b) / BD(a, a); % each of signal to noise ratio
      end
      for b = 1 : 5
        yi(ff, a) = sqrt(yi(ff, a)^2 + abs(yt(ff, a, b))^2); %total y_i
      end
      hGW(ff, a) = yi(ff, a) / BD(a, a); % signal to noise from total y_i
    end
  end  %end of ff for-loop

  
%  figure(1)
%  loglog(f, abs(hGW(:, 2))/Lcav, f, abs(hGWt(:, 2, 2))/Lcav, ':m', hclassic(:,1), hclassic(:,2), ':k', 'LineWidth', 3) % quantum total
%  legend('total', 'Lm only', 'classic noise');
%  xlabel('Frequency [Hz]')
%  ylabel('m/rHz')
%  axis([5 2000 1e-25 1e-19])
%  grid on

%  hold on
load qn.dat
  figure(1)
%  loglog(f, abs(hGWt(:, 2, 1))/Lcav, f, abs(hGWt(:, 2, 2))/Lcav, f, abs(hGWt(:, 2, 3))/Lcav, f, abs(hGWt(:, 2, 4))/Lcav, f, abs(hGWt(:, 2, 5))/Lcav, qn(:,1),qn(:,2),hclassic(:,1), hclassic(:,2), ':k', 'LineWidth', 3) % quantum total
  loglog(f, abs(hGWt(:, 2, 1))/Lcav, f, abs(hGWt(:, 2, 2))/Lcav, f, abs(hGWt(:, 2, 3))/Lcav, f, abs(hGWt(:, 2, 4))/Lcav, f, abs(hGWt(:, 2, 5))/Lcav, hclassic(:,1), hclassic(:,2), f, abs(hGW(:, 2))/Lcav,  ':k', 'LineWidth', 3) % quantum total
  legend('Lp', 'Lm', 'slp', 'slm', 'sls', 'classic noise', 'total');
  xlabel('Frequency [Hz]')
  ylabel('Loop noise [1/rHz]')
  axis([5 2000 1e-25 1e-19])
  grid on
  
  clear dat
   dat = [f', abs(hGW(:, 2))/Lcav, abs(hGWt(:, 2, 1))/Lcav, abs(hGWt(:, 2, 2))/Lcav, abs(hGWt(:, 2, 3))/Lcav, abs(hGWt(:, 2, 4))/Lcav, abs(hGWt(:, 2, 5))/Lcav]; 
   save('results2.dat', 'dat', '-ascii');