function Sensitivity

  % create the model
  clear
  load lcgt_optmat_tuned.mat
  %opt = lcgt; % opt can be changed in extractTF 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
%  f = logspace(0, 3, 10)';
  pos = zeros(opt.Ndrive, 1);
  %pos(nEX) = 0.1e-9;
  [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);
  
  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
    
  FF = zeros(5);
  FF(2,3) = -1.01; % FF gain = 100
  FF(2,4) = -1.01; % FF gain = 100
  FF(2,5) = -1.01; % FF gain = 100
%  FF(2,3) = -1.03; % FF gain = 30
%  FF(2,4) = -1.03; % FF gain = 30
%  FF(2,5) = -1.03; % FF gain = 30
  
%%% natural feed-forward terms
clear F1 F2 B1 B2
F1=0;
F2=0;
B1=0;
B2=0;
%    F1 = 74/3000;
%    F2 = 74/3000;
%    B1 = 1/2;
%    B2 = 1/2;
    Greal(ff,3,1) = Greal(ff,1,1) * F1; %FSS
    Greal(ff,5,1) = Greal(ff,1,1) * F2; %FSS
    Greal(ff,3,4) = - Greal(ff,4,4) * B1; %slm feedback to BS
    Greal(ff,5,4) = Greal(ff,4,4) * B2; %slm feedback to BS
%%%

%%% determnants
%%% invDijk=((invDii,invDij,invDik),(invDji,invDjj,invDjk),(invDki,invDkj,invDkk))

invD0 = inv(Dij);
invD0(:,1)=[];
invD0(:,1)=[];

invD123 = invD0;
invD123(5,:)=[];
invD123(4,:)=[];

invD124 = invD0;
invD124(5,:)=[];
invD124(3,:)=[];

invD125 = invD0;
invD125(4,:)=[];
invD125(3,:)=[];

invD134 = invD0;
invD134(5,:)=[];
invD134(2,:)=[];

invD145 = invD0;
invD145(2,:)=[];
invD145(2,:)=[];

invD234 = invD0;
invD234(5,:)=[];
invD234(1,:)=[];

invD235 = invD0;
invD235(4,:)=[];
invD235(1,:)=[];

invD245 = invD0;
invD245(3,:)=[];
invD245(1,:)=[];

invD345 = invD0;
invD345(1,:)=[];
invD345(1,:)=[];

%%% feed-forward

Greal(ff,2,3) = - Greal(ff,3,3) * FF(2, 3)* ( F2*det(invD124) - det(invD245) ) / ( F2*det(invD134) - det(invD345) + F1*det(invD145));
Greal(ff,2,4) = - Greal(ff,4,4) * FF(2, 4)* ( (F1*B2-F2*B1)*det(invD124) + F2*det(invD123) - F1*det(invD125) + B2*det(invD234) - det(invD235) + B1*det(invD245)) / ( -F2*det(invD134) + det(invD345) -F1*det(invD145) );
Greal(ff,2,5) = - Greal(ff,5,5) * FF(2, 5)* ( -F1*det(invD124) - det(invD234) ) / ( F2*det(invD134) - det(invD345) + F1*det(invD145));

    Greal2d = zeros(5);
    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
  
  save fname Greal
  
  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(2)
%  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), ':k', 'LineWidth', 3) % quantum total
  legend('Lp', 'Lm', 'slp', 'slm', 'sls', 'classic noise');
  xlabel('Frequency [Hz]')
  ylabel('Loop noise [1/rHz]')
  axis([5 2000 1e-25 1e-19])
  grid on
  
  clear dat
   dat = [f', 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('results.dat', 'dat', '-ascii'); 
  clear dat2