function p=parambKAGRAp1()
% Parameter set for bKAGRA Phase 1 (Mar 2018 configuration)
% 3-km PRMI with sapphire mirrors
% see http://gwwiki.icrr.u-tokyo.ac.jp/JGWwiki/LCGT/subgroup/ifo/MIF/OptParam
% see also http://gwwiki.icrr.u-tokyo.ac.jp/JGWwiki/KAGRA/Subgroups/IOO/OptParam
% Initiated by: Yuta Michimura

%% MODEL NAME
p.modelName = 'bKAGRAp1';

%% UNITS
ppm=1e-6;
MHz=1e6;
lambda=1.064e-06;
pm=1e-12;

%% LASER
p.Pin=2;

%% RF MODULATION
p.fmod1 = 16.880962*MHz;
p.fmod2 = 45.015898*MHz;
p.fmod3 = 56.269873*MHz;
p.vMod = generateRFFrequencyVector([p.fmod1],1);    % only f1 for bKAGRA phase 1
p.g1 = 0.1i;  % modulation depth (imaginary number for phase modulation)

%% FPMI MIRRORS
Nsilica = 1.44967;
Nsapphire = 1.754;  % refractive index of test masses

% BS
p.BS.aio = 45; % incident angle
p.BS.Chr = 0; % 1/ROC of BS
p.BS.Thr = 0.5; % BS transmission (ignore asymmetry for now)
p.BS.Lhr = 0; % BS HR Loss
p.BS.Rar = 0; % BS AR Reflection (for POB)
p.BS.Lmd = 0; % ignore the substrate loss
p.BS.Nmd = Nsilica; % refractive index

% ITMs
% no ITMs for bKAGRA Phase 1
p.ITMX.aio = 0; % incident angle
p.ITMX.Chr = 1.0/1900; % 1/ROC of ITMX
p.ITMX.Thr = 0.004; % ITMX transmission
p.ITMX.Lhr = 45*ppm; % ITMX HR Loss
p.ITMX.Rar = 0; % ITMX AR Reflection (for POX; ignore)
p.ITMX.Lmd = 0; % ignore the substrate loss
p.ITMX.Nmd = Nsapphire; % refractive index

p.ITMY.aio = 0;
p.ITMY.Chr = p.ITMX.Chr; % 1/ROC of ITMY
p.ITMY.Thr = 0.004; % ITMY transmission
p.ITMY.Lhr = 45*ppm; % ITMY HR Loss
p.ITMY.Rar = 0; %ITMY AR Reflection (for POY; ignore)
p.ITMY.Lmd = 0; % Ignore the substrate loss.
p.ITMY.Nmd = Nsapphire; % refractive index

% ETMs
p.ETMX.aio = 0; % incident angle
p.ETMX.Chr = 1.0/1900; % 1/ROC of ETMX
p.ETMX.Thr = 10*ppm; % ETMX transmission (for TRX)
p.ETMX.Lhr = 45*ppm; % ETMX HR Loss (ignore asymmetry for now)
p.ETMX.Rar = 0; % ETMX AR Reflection (ignore)
p.ETMX.Lmd = 0; % ignore the substrate loss
p.ETMX.Nmd = Nsapphire;  % refractive index
p.ETMX.dWaist = 1726.8; % this value determines the mode of incident beam

p.ETMY.aio = 0; % incident angle
p.ETMY.Chr = p.ETMX.Chr; % 1/ROC of ETMY
p.ETMY.Thr = 10*ppm; % ETMY transmission (for TRY)
p.ETMY.Lhr = 45*ppm; % ETMY HR Loss (ignore asymmetry for now)
p.ETMY.Rar = 0; % ETMY AR Reflection (ignore)
p.ETMY.Lmd = 0; % ignore the substrate loss
p.ETMY.Nmd = Nsapphire; % refractive index

%% RECYCLING CAVITY MIRRORS
% PRM
p.PRM.aio = 0; % incident angle
p.PRM.Chr = 1/458.128519465; % 1/ROC of PRM
p.PRM.Thr = 0.1; % PRM transmission
p.PRM.Lhr = 45*ppm; % PRM HR Loss
p.PRM.Rar = 0; % PRM AR Reflection (ignore)
p.PRM.Lmd = 0; % ignore the substrate loss
p.PRM.Nmd = 1.45; % refractive index

% PR2
p.PR2.aio = 0.6860; % incident angle
p.PR2.Chr = 1.0/-3.0764084715; % 1/ROC of PR2
p.PR2.Thr = 500*ppm; % PR2 transmission (for POP)
p.PR2.Lhr = 0; % PR2 HR Loss
p.PR2.Rar = 0; % PR2 AR Reflection (ignore)
p.PR2.Lmd = 0; % ignore the substrate loss
p.PR2.Nmd = Nsilica; % refractive index

% PR3
p.PR3.aio = 0.6860; % incident angle
p.PR3.Chr = 1.0/24.9164838708; % 1/ROC of PR3
p.PR3.Thr = 0; % PR3 transmission
p.PR3.Lhr = 0; % PR3 HR Loss
p.PR3.Rar = 0; % PR3 AR Reflection (ignore)
p.PR3.Lmd = 0; % ignore the substrate loss
p.PR3.Nmd = Nsilica; % refractive index

% SRM
% no SRM for bKAGRA Phase 1
p.SRM.aio = 0; % incident angle
p.SRM.Chr = 1.0/458.1285; % 1/ROC of SRM
p.SRM.Thr = 0.1536; % SRM transmission
p.SRM.Lhr = 0; % SRM HR Loss
p.SRM.Rar = 0; % SRM AR Reflection (ignore)
p.SRM.Lmd = 0; % ignore the substrate loss
p.SRM.Nmd = Nsilica; % refractive index

% SR2
p.SR2.aio = 0.6860; % incident angle
p.SR2.Chr = 1.0/-2.98718007727; % 1/ROC of SR2
p.SR2.Thr = 0; % SR2 transmission
p.SR2.Lhr = 0; % SR2 HR Loss
p.SR2.Rar = 0; % SR2 AR Reflection (ignore)
p.SR2.Lmd = 0; % ignore the substrate loss
p.SR2.Nmd = Nsilica; % refractive index

% SR3
p.SR3.aio = 0.6860;
p.SR3.Chr = 1.0/24.9164838708; % 1/ROC of SR3
p.SR3.Thr = 0; % SR3 transmission
p.SR3.Lhr = 0; % SR3 HR Loss
p.SR3.Rar = 0; % SR3 AR Reflection (ignore)
p.SR3.Lmd = 0; % ignore the substrate loss
p.SR3.Nmd = Nsilica; % refractive index

%% OTHER MIRRORS
% input modematching telescope
p.IMMT1.aio = 0.3;
p.IMMT1.Chr = 1/-8.953;
p.IMMT1.Thr = 0;
p.IMMT1.Lhr = 0;
p.IMMT1.Rar = 0;
p.IMMT1.Lmd = 0;
p.IMMT1.Nmd = Nsilica;

p.IMMT2.aio = 0.3;
p.IMMT2.Chr = 1/-13.910;
p.IMMT2.Thr = 0;
p.IMMT2.Lhr = 0;
p.IMMT2.Rar = 0;
p.IMMT2.Lmd = 0;
p.IMMT2.Nmd = Nsilica;

% beam splitters for splitting A and B
p.HALF.aio = 45;
p.HALF.Chr = 0;
p.HALF.Thr = 0.5;
p.HALF.Lhr = 0;
p.HALF.Rar = 0;
p.HALF.Lmd = 0;
p.HALF.Nmd = Nsilica;

%% LENGTHS 
% Michelson part
p.Las = 3.3299;  %Schnupp asymmetry
p.LMIavg = 25;  %Average length of the Michelson arms 
p.LBS_ITMX = p.LMIavg+p.Las/2; % Michelson X arm
p.LBS_ITMY = p.LMIavg-p.Las/2; % Michelson Y arm

% arm cavity length
p.Larm = 3000;

% PRC lengths
p.dLP2 = -0.024;%-0.0308; % move PR2 and PR3 each by p.dLP2/2(m) to stabilize 3-km PRMI cavity
p.LPRM_PR2 = 14.7614883609 + p.dLP2/2; %Distance between PRM and PR2
p.LPR2_PR3 = 11.0660636313 + p.dLP2; %Distance between PR2 and PR3
p.LPR3_BS = 15.7637759961 + p.dLP2/2; %Distance between PR3 and BS

% SRC lengths
p.LSRM_SR2 = 14.7412236675; %Distance between SRM and SR2
p.LSR2_SR3 = 11.1115048847; %Distance between SR2 and SR3
p.LSR3_BS = 15.7385994361; %Distance between SR3 and BS

% other lengths are all 0

%% ABCD MATRIX for PRMI
% ABCD matrix of each optical element (pitch)
p.PRM.ABCD = [1,0;
    -2*p.PRM.Chr,1];
p.PR2.ABCD = [1,0;
    -2*p.PR2.Chr*cosd(p.PR2.aio),1]; 
p.PR3.ABCD = [1,0;
    -2*p.PR3.Chr*cosd(p.PR3.aio),1];
p.ETMX.ABCD = [1,0;
    -2*p.ETMX.Chr,1];
p.Lp1.ABCD = [1,p.LPRM_PR2;
    0,1];
p.Lp2.ABCD = [1,p.LPR2_PR3;
    0,1];
p.Lp3.ABCD = [1,p.LPR3_BS+p.LBS_ITMX+p.Larm;
    0,1];

% roundtrip ABCD matrix
p.roundtrip.ABCD = p.Lp3.ABCD * p.PR3.ABCD * p.Lp2.ABCD * p.PR2.ABCD...
    * p.Lp1.ABCD * p.PRM.ABCD * p.Lp1.ABCD * p.PR2.ABCD * p.Lp2.ABCD...
    * p.PR3.ABCD * p.Lp3.ABCD * p.ETMX.ABCD;
p.A = p.roundtrip.ABCD(1,1);
p.B = p.roundtrip.ABCD(1,2);
p.C = p.roundtrip.ABCD(2,1);
p.D = p.roundtrip.ABCD(2,2);
p.criteria = (p.A+p.D)/2;   % abs(p.criteria) < 1  <=> stable 
p.z = real((p.A-p.D+sqrt((p.A+p.D)^2-4))/2/p.C);    % eigenmode's z at ETMX
p.z0 = imag((p.A-p.D+sqrt((p.A+p.D)^2-4))/2/p.C); 

% derivative of roundtrip Gouy phase
p.PR3.ABCDd = [0,0;
    2*p.PR3.Chr,0];
p.Lp1.ABCDd = [0,p.LPRM_PR2;
    0,0];

p.roundtrip.ABCDd = p.Lp3.ABCD * p.PR3.ABCD * p.Lp2.ABCD * p.PR2.ABCD...
    * p.Lp1.ABCDd * p.PRM.ABCD * p.Lp1.ABCD * p.PR2.ABCD * p.Lp2.ABCD...
    * p.PR3.ABCD * p.Lp3.ABCD * p.ETMX.ABCD...
    + p.Lp3.ABCD * p.PR3.ABCD * p.Lp2.ABCD * p.PR2.ABCD...
    * p.Lp1.ABCD * p.PRM.ABCD * p.Lp1.ABCDd * p.PR2.ABCD * p.Lp2.ABCD...
    * p.PR3.ABCD * p.Lp3.ABCD * p.ETMX.ABCD;
p.Ad = p.roundtrip.ABCDd(1,1);
p.Dd = p.roundtrip.ABCDd(2,2);
p.difcriteria = (p.Ad+p.Dd)/2; 

% roundtrip ABCD matrix (to obtain q of incident beam at PRM)
p.roundtrip.ABCD_PRM = p.Lp1.ABCD * p.PR2.ABCD * p.Lp2.ABCD * p.PR3.ABCD...
    * p.Lp3.ABCD * p.ETMX.ABCD * p.Lp3.ABCD * p.PR3.ABCD * p.Lp2.ABCD...
    * p.PR2.ABCD * p.Lp1.ABCD * p.PRM.ABCD;
p.A_PRM = p.roundtrip.ABCD_PRM(1,1);
p.B_PRM = p.roundtrip.ABCD_PRM(1,2);
p.C_PRM = p.roundtrip.ABCD_PRM(2,1);
p.D_PRM = p.roundtrip.ABCD_PRM(2,2);
p.z_PRM = real((p.A_PRM-p.D_PRM+sqrt((p.A_PRM+p.D_PRM)^2-4))/2/p.C_PRM);    % eigenmode's z at ETMX
p.z0_PRM = imag((p.A_PRM-p.D_PRM+sqrt((p.A_PRM+p.D_PRM)^2-4))/2/p.C_PRM); 
%p.q_PRM = -p.z_PRM + 1i*p.z0_PRM;
p.q_PRM = -9.1273 + 1i*64.0168;

p.profile = p.Lp3.ABCD * p.PR3.ABCD * p.Lp2.ABCD...
    * p.PR2.ABCD * p.Lp1.ABCD...
    *[p.q_PRM; 1];
p.profile_normalized = p.profile(1)/p.profile(2);

%% ATTENUATORS
% tuned later (about 50mW at each probe (except AS))
p.AttPOP = 0;
p.AttREFL = 0;
p.AttAS = 0;
p.AttTRX = 0;
p.AttTRY = 0;
p.AttREFLQPD = 0;
p.AttASQPD = 0;
p.AttPOPQPD = 0;
p.AttTRXQPD = 0;
p.AttTRYQPD = 0;

%% OPERATING POINTS
p.posOffsetPRM = 0;
p.armOffset = 0;        % DC readout

%% MECHANICAL TRANSFER FUNCTIONS (force/torque on TM to displacement/angle of TM)
% TO BE REPLACED WITH REAL SUSPENSION MODEL LATER
dampRes = [-0.1+1i, -0.1-1i];

% BS (Type-B)
p.BS.mass = 18.9;
p.BS.inertia = 1/12*p.BS.mass*(3*(0.37/2)^2+0.08^2);
p.BS.mechTFpos = zpk([],1*dampRes,1/p.BS.mass);
p.BS.mechTFpit = zpk([],1*dampRes,1/p.BS.inertia);
p.BS.mechTFyaw = zpk([],1*dampRes,1/p.BS.inertia);

% ITMs, ETMs (Type-A)
p.TM.mass = 22.7;
p.TM.inertia = 1/12*p.TM.mass*(3*(0.22/2)^2+0.15^2);
p.ITMX.mechTFpos = zpk([],1*dampRes,1/p.TM.mass);
p.ITMX.mechTFpit = zpk([],1*dampRes,1/p.TM.inertia);
p.ITMX.mechTFyaw = zpk([],1*dampRes,1/p.TM.inertia);
p.ITMY.mechTFpos = zpk([],1*dampRes,1/p.TM.mass);
p.ITMY.mechTFpit = zpk([],1*dampRes,1/p.TM.inertia);
p.ITMY.mechTFyaw = zpk([],1*dampRes,1/p.TM.inertia);
p.ETMX.mechTFpos = zpk([],1*dampRes,1/p.TM.mass);
p.ETMX.mechTFpit = zpk([],1*dampRes,1/p.TM.inertia);
p.ETMX.mechTFyaw = zpk([],1*dampRes,1/p.TM.inertia);
p.ETMY.mechTFpos = zpk([],1*dampRes,1/p.TM.mass);
p.ETMY.mechTFpit = zpk([],1*dampRes,1/p.TM.inertia);
p.ETMY.mechTFyaw = zpk([],1*dampRes,1/p.TM.inertia);

% PRs (Type-Bp)
p.PR.mass = 10.8;
p.PR.inertia = 1/12*p.PR.mass*(3*(0.25/2)^2+0.1^2);
p.PRM.mechTFpos = zpk([],1*dampRes,1/p.PR.mass);
p.PRM.mechTFpit = zpk([],1*dampRes,1/p.PR.inertia);
p.PRM.mechTFyaw = zpk([],1*dampRes,1/p.PR.inertia);
p.PR2.mechTFpos = zpk([],1*dampRes,1/p.PR.mass);
p.PR2.mechTFpit = zpk([],1*dampRes,1/p.PR.inertia);
p.PR2.mechTFyaw = zpk([],1*dampRes,1/p.PR.inertia);
p.PR3.mechTFpos = zpk([],1*dampRes,1/p.PR.mass);
p.PR3.mechTFpit = zpk([],1*dampRes,1/p.PR.inertia);
p.PR3.mechTFyaw = zpk([],1*dampRes,1/p.PR.inertia);

% SRs (Type-B)
p.SR.mass = 10.8;
p.SR.inertia = 1/12*p.SR.mass*(3*(0.25/2)^2+0.1^2);
p.SRM.mechTFpos = zpk([],1*dampRes,1/p.SR.mass);
p.SRM.mechTFpit = zpk([],1*dampRes,1/p.SR.inertia);
p.SRM.mechTFyaw = zpk([],1*dampRes,1/p.SR.inertia);
p.SR2.mechTFpos = zpk([],1*dampRes,1/p.SR.mass);
p.SR2.mechTFpit = zpk([],1*dampRes,1/p.SR.inertia);
p.SR2.mechTFyaw = zpk([],1*dampRes,1/p.SR.inertia);
p.SR3.mechTFpos = zpk([],1*dampRes,1/p.SR.mass);
p.SR3.mechTFpit = zpk([],1*dampRes,1/p.SR.inertia);
p.SR3.mechTFyaw = zpk([],1*dampRes,1/p.SR.inertia);

%% Frequency FOR tickle, tickle01
p.ftickle = 100; % frequency used for tickle
p.ftickle01 = 10; % frequency used for tickle01

%% DEMOD PHASES
% tuned later
% Some constants for calculating demodulation phases
c = 299792458;
omegamod1 = 2*pi*p.fmod1;
omegamod2 = 2*pi*p.fmod2;

% demod phases (degrees); TUNED LATER
p.demodphasePOP1 = 0;      % demod phase for POP f1 demodulation I phase
p.demodphaseREFL1 = -1.1; %-26.3;
p.demodphaseAS1 = 21.4;

%% GOUY PHASES  
% Optimal Gouy phases
p.GouyREFLB = -5.57*pi/180;%-3.86*pi/180;%-5.633*pi/180;
p.GouyREFLA = p.GouyREFLB + pi/2;%-19.90*pi/180;    % Gouy phase for QPDA at REFL
%p.GouyREFLB = p.GouyREFLA + acos(p.criteria)/2;
p.GouyASA = 27.7*pi/180;%34.5*pi/180;%27.24*pi/180;      % Gouy phase for QPDA at AS
p.GouyASB = p.GouyASA + pi/2;


p.GouyPOPA = -82.5*pi/180;%-79.0*pi/180;%-82.77*pi/180;%-83.01*pi/180;%-80.42*pi/180;     % Gouy phase for QPDA at POP
p.GouyPOPB = p.GouyPOPA + pi/2;

p.GouyTRXA = 84.5*pi/180;%86.1*pi/180;%83.97*pi/180;%85.53*pi/180;     % Gouy phase for QPDA at TRX
p.GouyTRXB = p.GouyTRXA + pi/2;
p.GouyTRYA = 84.5*pi/180;%86.1*pi/180;%83.97*pi/180;%85.53*pi/180;     % Gouy phase for QPDA at TRY
p.GouyTRYB = p.GouyTRYA + pi/2;