function dat=noise_lcgt(parm,bw,prg,sbg);

%  noise.m
%  calculate LCGT sensitivity
%           June 24, 2001
%           Masaki Ando

if nargin==0
%    parm=802.249e3; bw=202.9999; prg=10; sbg=10; 
   parm=771.7754e3;
   bw=236.9261;
   prg=10.6285;
   sbg=14.6838; 

%else 
%    prg=1;
end;

p0=75*(parm/771.775e3)*(bw/236.9261*14.6838/sbg)*(10.6285/prg);
finesse=249.827*100./bw*sbg;

t0=300;

fr_i=[-1:0.01:5]';
fr=10.^fr_i;
%fr=[1:5:5e3]';

[num,c]=size(fr);

%

kb=1.381e-23;                % Boltzmann's constant [J/K]
g=9.8;                       % Acceleration of the gravity [m/s^2]
c=2.99792458e8;              % Speed of light [m/s]
hb=1.05457266912510183e-34;  % Reduced Planck constant [J s]
eta=0.9;                     % Quantum efficiency [electron/photon]

% LCGT parameters

%m=50;                       % Mass of a mirror [kg]
l=3000;                      % Length of an arm [m]
lambda=1064e-9;              % Laser wavelength [m]
%p0=100;                      % Laser power [W]
%prg=50;                        % Power recycling gain  
w0=3e-2;                     % Beam radius on a mirror [m]
%finesse=100;                 % Finesse

% Seismic noise

Gseismic=1e-9./fr.^2;        % Seismic motion [m/Hz^(1/2)]
%Hiso=10^(-3)./fr.^10;        % total isolation ratio 
Hiso=10^(-3)./fr.^6;  
hseis=2./l.*Gseismic.*Hiso;


% Mirror Thermal noise

qm=1e8;                       % Q-value of a mirror mass
rho=4.0*10^3;                 % density [kg/m^3]
%radius=15*10^(-2);            % radius [m]
%height=18*10^(-2);            % height [m] 
radius=12.5*10^(-2);            
height=15*10^(-2);            
mass= rho*pi*radius^2*height; % Mass [kg] 
E0=4.0e11;               % Young's modulus [Pa] 
sigma = 0.29;                 % Poisson ratio 
dcoa=5e-6;             % Thickness of coating [m] 
phicoa=4e-4;               % Loss angle of coating
Tm=20;                        % Temperature of mirror [K] 

phimir=1/qm;

hmirrorhomo=2/l*sqrt(4*kb*Tm*(1-sigma^2)*phimir./(sqrt(pi)*E0*w0*2*pi.*fr));     

hthermoe = 9.4e-25.*(Tm/20).*(rho./4.0e3).^(1/4).*(3e3./l).*(100./fr).^(1/4);

hmirrorcoa= 1.2e-24.*(phicoa./4e-4).^(1/2).*(dcoa./5e-6).^(1/2) ...
  .*(Tm./20).^(1/2).*(3e-2./w0).*...
   (4e11/E0).^(1/2).*(3e3/l).*(100./fr).^(1/2);

hmirror=sqrt(hmirrorhomo.^2+ hthermoe.^2 + hmirrorcoa.^2);      


% Pendulum Thermal noise

qp=2e8;                      % Q-value of a pendulum
%qv=qp/2;
lsus=0.40;                   % Length of wire [m]
fp=1/(2*pi)*sqrt(g/lsus);    % Resonant frequency of suspension [Hz]
tp=13;                       % Temperature of a pendulum [K]

hpen=2/l*sqrt(4*kb*tp./mass./qp./(2*pi)^3./fp.^2./fr./...
                abs(1+i./qp-(fr./fp).^2).^2);

%for j=1:1000
%   fn=514*j;
%   mn=1.034/2*(fn/fp)^2;
%   th_p=sqrt(th_p.^2+4*kb*t./mn./qv./(2*pi)^3./fn.^2./fr./...
%      abs(1+i./qv-(fr./fn).^2).^2);
%end

% Shot noise

pin=p0*prg;
fcut=c/4/l/finesse*sbg;
tau=2*l*finesse/(pi*c)/sbg;

%hshot=sqrt(hb*lambda/(4*pi*c*pin*tau^2).*(1+(fr./fcav).^2));
%hshot=sqrt(hb*lambda/(4*pi*c*pin*tau^2).*(1+(fr./fcav).^2))*sqrt(3/2);
hshot=sqrt(hb*lambda*fcut/(2*parm*l).*(1+(fr./fcut).^2))*sqrt(3/2)/sqrt(eta);

% Radiation pressure noise

%hradi=4/l*(finesse/sbg/pi)./(mass*(2*pi.*fr).^2)...
%    .*sqrt(16.*pi.*hb.*pin./(c.*lambda)./(1+(fr./fcut).^2));

hradi=4/l./(mass*(2*pi.*fr).^2)...
   .*sqrt(2*hb.*parm./(lambda*l*fcut)./(1+(fr./fcut).^2))/sqrt(eta);

% Total noise

tot=sqrt(hmirror.^2+hpen.^2+hshot.^2+hradi.^2+hseis.^2);
%tot=sqrt(hshot.^2);

dat=[fr,hseis,hpen,hmirror,hradi,hshot,tot];

if nargin==0
%   Save and plot data

   loglog(fr,hmirrorhomo,fr,hpen,fr,hradi,fr,hshot,fr,hseis,fr,hmirror,fr,tot,'k.')
   axis([1,1e5,1e-25,1e-18])
   grid on

   save lcgt050704.spe dat -ascii
end;



% %%%%%%%%%%%
% 
% T=2*pi/finesse*sbg;
% gam=T*c/(4*l);
% %I0= 1000*T/0.0511;
% I0= 10000*T/0.0511;
% 
% ISQL= mass*l^2*gam^4/4/(c/lambda*2*pi);
% Omega=2*pi*fr;
% beta=atan(Omega/gam);
% kappa=2*(I0/ISQL)*gam^4./(Omega.^2.*(gam.^2+Omega.^2));
% 
% hSQL=sqrt(8*hb./(mass*(Omega*l).^2));
% 
% Phi=0;
% phi=pi/2;  %%;
% 
% rho0=0.; %%
% tau0=sqrt(1-rho0);
% 
% c11=(1+rho0^2)*(cos(2*phi)+kappa/2*sin(2*phi))-2*rho0*cos(2*(beta+Phi));
% c12=-tau0^2*(sin(2*phi)+kappa*sin(phi)^2);
% c21= tau0^2*(sin(2*phi)-kappa*cos(phi)^2);
% c22=c11;
% 
% D1=-( 1+rho0*exp(2*i * (beta+Phi))) *sin(phi);
% D2=-(-1+rho0*exp(2*i * (beta+Phi))) *cos(phi);
% 
% zeta=pi/2;
% 
% Sh= hSQL.^2./(2* kappa) ...
%      .* ( (c11 * sin(zeta) + c21 * cos(zeta)).^2 + ...
%           (c12 * sin(zeta) + c22 * cos(zeta)).^2 ) ...
%           ./(tau0^2 * abs(D1 * sin(zeta) + D2 * cos(zeta)).^2);
% 
% if nargin==0
%    hold on
%    loglog(fr,sqrt(Sh)); grid on
%    hold off
% end;