function dat=lcgt_param(Rf)

% Parameters of LCGT Interferometer 
%      ( conventional configulation ,  finite contrast)
%     Original :  Masaki Ando  June 18, 2001 for Mathematica
%     Modified :            October 02, 2002 for Matlab

% Constants
format short e

%if nargin==0; Rf=0.995; end;
if nargin==0; Rf=0.996; end;

lambda = 1064e-9;              % Wavelength of a laser source [m]
c = 2.99792458e8;              % Speed of light [m/s]
ee = 1.60217733e-19;           % Elementary charge [C]
L = 3000;                      % Baseline length [m]
%P = 100;                       % Lase power [W]
P = 75; 
eta = 1;                       % Quantum efficiency [A/W]
gdet = 1e3 * sqrt(2);          % Photo detector effective impedance [Ohm]
El = sqrt(P * eta);            % Effective laser power [A]
Omega = (2 * pi * c)/lambda;   % Angular frequency of laser [Hz]

% Losses

contrast = 0.995; 
CavRefLoss  = 1000 * 1e-6; 
loss = 10 * 1e-6;
BSloss = 100 * 1e-6; 
ARloss = 0.001;
BSAR  = ARloss;
Rst = 1; 
Tcont = (1 - contrast)/(1 + contrast); 
Rcont = (2 * contrast)/(1 + contrast); 
tcont = sqrt(Tcont); 
rcont  = sqrt(Rcont); 
tcavloss = sqrt(1 - CavRefLoss/2); 
rst = sqrt(Rst); 

% Mirrors

%rf = sqrt(0.994);       % Front mirror
rf = sqrt(Rf); 
re = sqrt(0.99995);     % End mirror
%rr = sqrt(0.76);        % PRM
%rs = sqrt(0.68);        % SEM
rr = sqrt(0.80);        
rs = sqrt(0.77);
rp = sqrt(0.001);       % Pick-off mirror

tf = sqrt(1 - rf^2 - loss); 
te = sqrt(1 - re^2 - loss);
tr = sqrt(1 - rr^2 - loss); 
ts = sqrt(1 - rs^2 - loss); 
tp = sqrt((1 - rp^2) * (1 - loss) * (1 - ARloss)); 

tar = sqrt(1 - ARloss);
tbs = sqrt((1 - BSAR) * (1 - BSloss)); 
tst = sqrt(1 - Rst);

% Modulation

wm = 15e6;             % Modulation frequency [rad]
delta = 0.25;          % Asymmetry [m]
m = 0.9;               % Modulation index [rad]
j0 = BesselJ(0, m);    % Bessel function
j1 = BesselJ(1, m);    
Con = El^2 * j0 * j1; 
alpha = (2 *pi * wm * delta)/c; 
lplus = 5;            % Recycling cavity length

% Caliculated parameters
%% Parameters of arm cavity reflectivity etc.

rreso = tcavloss^2 * tar^2 * (-rf + tf^2 * re/(1 - rf * re) );
ranti = tcavloss^2 * tar^2 * (-rf - tf^2 * re/(1 + rf * re) );
rdreso =tcavloss^2 * tar^2 * tf^2 * re / (1 - rf * re)^2;
rdanti =tcavloss^2 * tar^2 * tf^2 * re / (1 + rf * re)^2;
gcav =  tcavloss * tar * tf / (1 - rf * re);

%% Finesse etc.

fine = pi * sqrt(rf * re) / (1 - rf * re);
tau = 2 * L * fine / (pi * c);
cutoff = 1 / (tau * 2 * pi);
NFP = 2 * fine / pi;

%% parameters of Fabry-Perot-Michelson inteferometer reflectivity

rfpm0 =  tbs^2 * rreso * rcont;
rfpm1 = -tbs^2 * ranti * cos(alpha) * rcont;
rfpmp0 = tbs^2 * tp^2 * rreso * rcont * rst^4;
rfpmp1 =-tbs^2 * tp^2 * ranti * cos(alpha) * rcont * rst^4;
% 
% % Optimal recycling gain case
% rr = abs(rfpmp0);
% tr = sqrt(1 - rr^2 - loss);
% ((rr^2 + tr^2) * rfpmp1)^2;

% Parameters of power recycling 
%% Gain

g0 = tar * tr / (1 - rr * rfpmp0);
g1 = tar * tr / (1 - rr * rfpmp1);
Gp = g0^2;
Gs = g0 * g1;
Pbs = g0^2 * P;
Pcav = tbs * tp * gcav^2 * Pbs;

%% Reflectivity

rrfpm0 = tar^2 * (-rr + tr^2 * rfpmp0 / (1 - rr * rfpmp0) );
rrfpm1 = tar^2 * (-rr + tr^2 * rfpmp1 / (1 - rr * rfpmp1) );

% Parameters of signal recycling 

ers = tcavloss^2 * tar^2 * tbs^2 * rcont * rs;  
ets = tcavloss * tar * tbs * rcont * ts;

erf = rf - tf^2 * ers/(1 - rf * ers);
etf= tf * ets /(1 - rf * ers);

efine = pi * sqrt(erf * re) / (1 - erf * re);
etau = 2 * L * efine / (pi * c);
ecutoff = 1 / (etau * 2 * pi);
eNFP = 2 * efine / pi;

sbg = fine/efine;                % Signal bandwidth gain
sigloss = 1 - etf^2 / (1 - (erf * re)^2);    % Signal loss

% Signal
%% Signal size

%%% v1
v1Lm = Con * g0 * g1 * abs(rdreso) * ranti * sin(alpha) * tp^2 * tbs^4 * Rcont * rst^4;
v1ellm =  Con * g0 * g1 * rreso * ranti * sin(alpha) * tp^2 * tbs^4 * Rcont * rst^4;

%%% v2
v2Lm = Con * rrfpm0 * g1^2 * abs(rdanti) * sin(alpha) * tp^2 * tbs^2 * rcont * rst^4;
v2ellm = -Con * rrfpm0 * g1^2 * ranti * sin(alpha) * tp^2 * tbs^2 * rcont * rst^4;

%%% v3
v3Lp = Con * g0^2 * g1 * abs(rdreso) * ranti * cos(alpha) ...
   * tp^2 * rp^2 * tbs^4/tr * Rcont * rst^4;
v3ellp = Con * g0 * g1 * rreso * ranti * (g0 - g1) * cos(alpha) ... 
   * tp^2 * rp^2 * tbs^4/ tr * Rcont * rst^4;

%%% v4
v4Lp = Con * ( -rrfpm1 * g0^2 * abs(rdreso) + rrfpm0 * g1^2 * abs(rdanti)) ...
   * tp^2 * tbs^2 * rcont * rst^4;
v4ellp = -Con * (rrfpm1 * g0^2 * rreso + g1^2 * rrfpm0 * ranti * cos(alpha)) ... 
   * tp^2 * tbs^2 * rcont * rst^4;

%% Calibration

alpha1 = (v1Lm   * gdet * 4 * pi)/lambda;
alpha2 = (v2ellm * gdet * 4 * pi)/lambda;
alpha4 = (v4Lp   * gdet * 4 * pi)/lambda;
alpha3 = (v3ellp * gdet * 4 * pi)/lambda;

%% Signal ratio

sig = (1/v1Lm) * [[v1Lm, v1ellm,    0,      0]; ...
       [v2Lm, v2ellm,    0,      0]; ...
       [   0,      0, v3Lp, v3ellp]; ...
       [   0,      0, v4Lp, v4ellp]];
sigr = [(1/v1Lm) * [v1Lm, v1ellm, 0, 0]; ... 
       (1/v2ellm)* [v2Lm, v2ellm, 0, 0]; ...
       (1/v3Lp)  * [0, 0, v3Lp, v3ellp]; ... 
       (1/v4ellp)* [0, 0, v4Lp, v4ellp]];

% Shot noise
%% Shot noise

v1shot = sqrt( ee * El^2 * tp^2 * tbs^4 * Rcont * rst^2)...
   * sqrt(        2 * ( j1 * g1 * ranti * sin(alpha) )^2 ...
   + Tcont * ( ( j0 * g0 * rreso )^2  ...
   + 2 * ( j1 * g1 * ranti * cos(alpha) )^2 ) );
v2shot = sqrt( ee * ( (j0 * rrfpm0)^2 + 2 * (j1 * rrfpm1)^2 ) * El^2);
v3shot = sqrt( ee * ( (j0 * g0 * rreso)^2 + 2 * ( j1 * g1 * ranti * cos(alpha) )^2) ...
   * rp^2 * tp^2 * tbs^4 * rst^2 * El^2 * Rcont);
v4shot =  v2shot;

%% Shot noise ( V/rt[Hz] )
v1shot * gdet ;
v2shot * gdet ;
v3shot * gdet ;
v4shot * gdet ;

%% Shot noise level ( m/rt[Hz] )
shotlevel = lambda/(4 * pi) ... 
   * [[abs(v1shot/v1Lm), abs(v1shot/v1ellm), 0, 0]; ...
       [abs(v2shot/v2Lm), abs(v2shot/v2ellm), 0, 0]; ...
       [ 0, 0, abs(v3shot/v3Lp), abs(v3shot/v3ellp)]; ... 
       [ 0, 0, abs(v4shot/v4Lp), abs(v4shot/v4ellp)]];

% Frequency response
%% Cavity cut-off

omegac = (1 - rf * re)/(rf * re * 2 * L) * c;
nuc = omegac/(2 * pi);

%% Coupled cavity cut-off

omegacc = (1 - rr * rfpmp0) * omegac / (1 + rr * rfpmp0); 
nucc =omegacc/(2 *pi);

%% Recycling cavity cut-off

omegarec = - (1 + rr * tp^2 * ranti * cos(alpha) * rcont) ...
   / ( (rr * tp^2 * ranti * cos(alpha) * lplus)/c) ;
nurec = omegarec/(2 *pi);

dat=[Rf,Pcav, Pcav/g0^2 , fine, g0^2];

if nargin==0,

   % Parameters
   %% Print Out
   disp(' ')
   disp(['Laser Power                 : ', num2str(P), ' W'] ); 
   disp(['Wave Length                 : ', num2str(lambda * 1e9),' nm']); 
   disp(['Reflectivity of the Mirrors']); 
   disp(['   Front Mirror             : ', num2str(100 * rf^2), ' %']); 
   disp(['   End Mirror               : ', num2str(100 * re^2), ' %']); 
   disp(['   Power Recycling Mirror   : ', num2str(100 * rr^2), ' %']); 
   disp(['   Signal Extraction Mirror : ', num2str(100 * rs^2), ' %']); 
   disp(['   Pick-off Mirror          : ', num2str(100 * rp^2), ' %']);
   disp(['   Steering Mirror          : ', num2str(100 * Rst),  ' %']); 
   disp(['Modulation Frequency        : ', num2str(wm/1e6), ' MHz']); 
   disp(['Mudulation Index            : ', num2str(m)]); 
   disp(['Asymmetry                   : ', num2str(delta), ' m']); 
   disp(['The ƒ¿ parameter            : ', num2str(alpha)]); 
   disp(['Effective Mudulation Index  : ', num2str(m*sin(alpha))]); 
   disp(['Recycing Cavity Length      : ', num2str(lplus), ' m']); 
   disp(['Loss in an Optical Component: ', num2str(loss * 1e6), ' ppm']); 
   disp(['Reflectivity of an AR coat  : ', num2str(100 * ARloss), ' %']); 
   disp(['Beam Splitter AR coat       : ', num2str(100 * BSAR), ' %']); 
   disp(['Contrast of the Fringe      : ', num2str(100 * contrast), ' %']); 
   disp(['Arm Cavity']);
   disp(['    Length                  : ', num2str(L), ' m']); 
   disp(['    Reflectivity']); 
   disp(['          for the Carrier   : ', num2str(100 * rreso^2), ' %']); 
   disp(['          for the Sidebands : ', num2str(100 * ranti^2), ' %']); 
   disp(['    Power in a cavity       : ', num2str(Pcav/1e3), ' kW']); 
   disp(['    Phase Change Enhancement ']); 
   disp(['          for the Carrier   : ', num2str(rdreso)]); 
   disp(['          for the Sidebands : ', num2str(rdanti)]); 
   disp(['    Finesse                 : ', num2str(fine)]); 
   disp(['    Cut-off Frequency       : ', num2str(cutoff), ' Hz']); 
   disp(['    Loss on the Reflection  : ', num2str(100 * CavRefLoss),' %']); 
   disp(['Fabry-Perot-Michelson interferometer']);
   disp(['    Reflectivity']); 
   disp(['          for the Carrier   : ', num2str(100 * rfpmp0^2),' %']); 
   disp(['          for the Sidebands : ', num2str(100 * rfpmp1^2),' %']); 
   disp(['Power Recycled Interferometer']); 
   disp(['    Recycling Gain']); 
   disp(['          for the Carrier   : ', num2str(g0^2)]); 
   disp(['          for the Sidebands : ', num2str(g1^2)]); 
   disp(['    Power on BS             : ', num2str(Pbs),' W']); 
   disp(['    Reflectivity']); 
   disp(['          for the Carrier   : ', num2str(100 * rrfpm0^2), ' %']); 
   disp(['          for the Sidebands : ', num2str(100 * rrfpm1^2),  ' %']); 
   disp(['    Coupled Cavity cutoff   : ', num2str(nucc), ' Hz']); 
   disp(['    Recycling Cavity cutoff : ', num2str(nurec/1e6),  ' MHz']); 
   disp(['Signal Extraction Cavity']); 
   disp(['    Compound Front Mirror   : ', num2str(100 * erf^2), ' %']); 
   disp(['    Signal Band Gain        : ', num2str(sbg)]);
   disp(['    Signal band             : ', num2str(ecutoff), ' Hz']); 
   disp(['    Signal loss             : ', num2str(100 * sigloss ), ' %']); 

%     disp(['Signals']); 
%     disp(['    Signal Ratio Matrix : ']); 
%     disp(sig);
%     disp(['    Signal Ratio Matrix with nomalization: ']); 
%     disp(sigr); 
%     disp(['    Shot-noise Level : ', ' m/sqrt(Hz)']); 
%     disp(shotlevel); 
end;
