% Frequency-Domain LIGO1 Interferometer Model
%  computes frequency dependent IFO response matrix
% m = freqIFO(p, f);
%
% [AS_Q; REFL_Q; POB_Q; AS_I; REFL_I; POB_I] = ...
%  m * [ETMx; ETMy; ITMx; ITMy; RM; BS; dFreq]
%
% Giving a third argument 'DOF' will cause the response matrix to be
% in the canonical degree-of-freedom basis
% m = freqIFO(p, f, 'DOF');
% [AS_Q; REFL_Q; POB_Q; AS_I; REFL_I; POB_I] = ...
%  m * [DARM; CARM; MICH; PRC; dFreq]
% where DARM = Lx - Ly, CARM = Lx + Ly, MICH = lx - ly, PRC  = lx + ly
% ('DRV' gives a similar matrix, but with a DOF drive choice like H1)..
%
% A third argument of 'COMP' returns m in the computational basis,
% where BS motion is taken to be in the direction of the RM
% rather than along its HR surface normal (which is the convention
% used for other masses).
% m = freqIFO(p, f, 'COMP');
% [AS_Q; REFL_Q; POB_Q; AS_I; REFL_I; POB_I] = ...
%  m * [ETMx; ETMy; ITMx; ITMy; RM; sqrt(2) * BS; dFreq]
%
% A forth argument of 'CM' uses dFreq to force REFL_I to zero.
% ('CM' can also be used as a third argument, in which case the
% default basis is used.)
%
% % ==== Parameter Struct (and possible values)
% % == Optical Path Lengths between HR Surfaces
% p.lArmX = 3995;			% X arm length [m]
% p.lArmY = 3995;			% Y arm length [m]
% p.lRMbs = 4.397;			% RM to BS optical path length [m]
% p.lBSx  = 4.9765;			% BS to ITMX optical path length [m]
% p.lBSy  = 4.5975;			% BS to ITMY optical path length [m]
%
% % == HR Surface Power Transmission
% p.Tix = 0.0275;			% ITMx power transmission
% p.Tiy = 0.028;			% ITMy power transmission
% p.Trm = 0.027;			% RM power transmission 
% p.Tbs = 0.49;				% BS power transmission
%
% % == HR Surface Power Loss
% p.Lex = 110e-6;			% ETMx power loss
% p.Ley =  40e-6;			% ETMy power loss
% p.Lix = 70e-6;			% ITMx power loss
% p.Liy = 70e-6;			% ITMy power loss
% p.Lrm = 100e-6;			% RM power loss
% p.Lbs = 100e-6;			% BS power loss
%
% % == Bulk Media and AR Surfaces Losses (single pass)
% p.LBix = 100e-6;			% ITMx back side power loss
% p.LBiy = 100e-6;			% ITMy back side power loss
% p.LBrm = 100e-6;			% RM back side power loss
% p.LBbs = 100e-6;			% BS back side power loss
%
% % == Input Beam
% p.lambda = 1064;			% wave length [nm]
% p.ampMod = 0.44;			% input modulation depth
% p.freqMod = 24.482e6;			% input modulation frequency (Hz)
% p.pCR = -1;				% (NPTRX + NPTRY) / 2, or -1 for max
% p.pSB = -1;				% NSPOB, or -1 for max
%
% % == Demod Phases
% p.phiAS   = 0;			% AS port demod phase   [degree]
% p.phiREFL = 0;			% REFL port demod phase [degree]
% p.phiPOB  = 0;			% POB port demod phase  [degree]
%
% The normalized powers, pCR and pSB are used to compute
% the input coupling.  A value of -1 uses perfect coupling.
%
% This function can be called with up to 5 output arguments as
% [m, pm, mCR, mP1, mP2] = freqIFO(p, f);
% in which case
% m   - is the total response matrx, as before
% pm  - a modified parameter struct
% mCR - is the carrier signal matrix, in the COMP basis
% mP1 - is the first SB signal matrix, in the COMP basis
% mP2 - is the second SB signal matrix, in the COMP basis

% ================================================================
% Some relationships that may help determine arm
% cavity parameters from measured data:
%  rRT = rITM * rETM			% arm round-trip reflectivity
%  tau = -2 * lArm / (c * log(rRT))	% arm decay time
%  fCP = 1 / (2 * pi * tau)		% arm cavity pole frequency
%
% solving for rRT gives...
%  rRT = exp(-4 * pi * lArm * fCP / c)
% ================================================================

% ================================================================
% To Do?
% Include lower SBs (pM1 and pM2)
%  pM1 = freqProp(p, -p.freqMod);
%  pM2 = freqProp(p, -2 * p.freqMod);
% Include some aIn asymmetry between P1 and M1
%  This will make pP1 ~= conj(pM1), which will result in offsets
%  in output signals, which creates a hell of a mess.  (Mirrors
%  should be offset to zero the signal offsets, this involves
%  control loops...the mess grows.)
% Add higher-modes...at least the even ones
%  This too will introduce offsets, and a nasty mess with them.
% ================================================================

function [m, p, mCR, mP1, mP2] = freqIFO(p, f, varargin)

  % ==== Additional Parameters
  % things with units
  p.c = 299792458;			% speed of light [m/s]
  p.k = 2 * pi / p.lambda;		% wave number [1/nm]

  % front-side amplitude reflectivity
  p.rEX = sqrt(1 - p.Lex);
  p.rEY = sqrt(1 - p.Ley);
  p.rIX = sqrt(1 - p.Tix - p.Lix);
  p.rIY = sqrt(1 - p.Tiy - p.Liy);
  p.rRM = sqrt(1 - p.Trm - p.Lrm);
  p.rBS = sqrt(1 - p.Tbs - p.Lbs);

  % back-side amplitude reflectivity
  p.bIX = -p.rIX * (1 - p.LBix);
  p.bIY = -p.rIY * (1 - p.LBiy);
  p.bRM = -p.rRM * (1 - p.LBrm);
  p.bBS = -p.rBS * (1 - p.LBbs);

  % transmissivity
  p.tIX = sqrt(p.Tix * (1 - p.LBix));
  p.tIY = sqrt(p.Tiy * (1 - p.LBiy));
  p.tRM = sqrt(p.Trm * (1 - p.LBrm));
  p.tBS = sqrt(p.Tbs * (1 - p.LBbs));

  % ==== DC Fields
  pCR = freqProp(p, 0);
  pP1 = freqProp(p,  p.freqMod);
  pP2 = freqProp(p,  2 * p.freqMod);

  % PRC input carrier coupling
  pCRmax = abs(pCR.gP);
  if( p.pCR > 0 )
    pCR.aIn = sqrt(p.pCR) / pCRmax;		% coupling to input beam
    fprintf(1, 'Carrier input coupling is %.1f%%\n', pCR.aIn^2 * 100);
  else
    pCR.aIn = 1;
    fprintf(1, 'Maximum nPtr is %.1f\n', pCRmax^2);
  end

  % PRC input sideband coupling
  pSBmax = abs(pP1.gP);
  if( p.pSB > 0 )
    pP1.aIn = sqrt(p.pSB) / pSBmax;		% coupling to input beam
    fprintf(1, 'Sideband input coupling is %.1f%%\n', pP1.aIn^2 * 100);
  else
    pP1.aIn = 1;
    fprintf(1, 'Maximum nSpob is %.1f\n', pSBmax^2);
  end

  % already too much about 2nd SB, don't worry about coupling
  pP2.aIn = 1;

  % field amplitudes
  pCR = fieldAmp(p, pCR, bessel(0, p.ampMod));
  pP1 = fieldAmp(p, pP1, bessel(1, p.ampMod));
  pP2 = fieldAmp(p, pP2, bessel(2, p.ampMod));

  % ==== Signals from Various Sources
  % add DC to f vector, for reference in demodulation
  ff = [f(:); 0];

  % compute signals
  % (assumes upper and lower audio SBs are equivalent)
  pf = freqProp(p, ff);
  mCR = sigMatrix(pCR, pf);

  pf = freqProp(p, ff + p.freqMod);
  mP1 = sigMatrix(pP1, pf);

  pf = freqProp(p, ff + 2 * p.freqMod);
  mP2 = sigMatrix(pP2, pf);

  % ==== Demodulation
  % mode overlaps
  % This computation assumes that the cavity modes
  % overlap only in the input mode.  There are actually
  % two free overlap parameters and a short paper worth
  % of calculations here...for now I'll stick with the
  % simplest assumption.

  a01 = pCR.aIn * pP1.aIn;			% CR-SB cavity mode overlap
  a01r = pCR.aIn;				% CR-SB reject mode overlap
  a10r = pP1.aIn;				% SB-CR reject mode overlap

  % == local oscillators
  % SB at each port, overlapped with CR-cavity mode
  e01(1) = pP1.eAS * a01;			% at AS port
  e01(2) = [pCR.aIn, a01, a01r] * pP1.eREFL;	% at REFL port
  e01(3) = pP1.ePOB * a01;			% at POB port

  % CR at each port, overlapped with SB-cavity mode
  e10(1) = pCR.eAS * a01;
  e10(2) = [pP1.aIn, a01, a10r] * pCR.eREFL;
  e10(3) = pCR.ePOB * a01;

  % 2nd SB at each port, overlapped with SB-cavity mode
  e12(1) = pP2.eAS * pP1.aIn;
  e12(2) = [pP1.aIn, pP1.aIn, 0] * pP2.eREFL;
  e12(3) = pP2.ePOB * pP1.aIn;

  % SB at each port, overlapped with input mode
  e21(1) = pP1.eAS * pP1.aIn;
  e21(2) = [1, pP1.aIn, 1] * pP1.eREFL;
  e21(3) = pP1.ePOB * pP1.aIn;

  % == demod phases
  phi(1) = (pi/180) * p.phiAS + angle(e01(1)) - pi / 2;
  phi(2) = (pi/180) * p.phiREFL + angle(e01(2));
  phi(3) = (pi/180) * p.phiPOB + angle(e01(3));

  % == demod computation
  % The factor of 4 comes from 2 factors of 2.  The first is from
  % upper and lower RF SBs, the second is the fact that mirrors
  % moving generate phi = (2 * k * z) on reflection.  Note that the
  % "Readout Paper" has a factor of 2 in S_0, which is the first
  % factor listed here, so the factor of 4 appearing in all signal
  % equations actually includes an additional factor of 2 for the
  % demodulation that has not yet been applied to those signals.

  m = zeros(length(f), 6, 7);
  for n = 1:3
    for nn = 1:7
      [aI01, aQ01] = demod(mCR(:, n, nn), e01(n),  1, phi(n));
      [aI10, aQ10] = demod(mP1(:, n, nn), e10(n), -1, phi(n));
      [aI12, aQ12] = demod(mP1(:, n, nn), e12(n),  1, phi(n));
      [aI21, aQ21] = demod(mP2(:, n, nn), e21(n), -1, phi(n));
      m(:, n    , nn) = (aQ01 + aQ10 + aQ12 + aQ21) * (4 * p.k);
      m(:, n + 3, nn) = (aI01 + aI10 + aI12 + aI21) * (4 * p.k);

      % uncomment below to eliminate 2nd order SB signals
      %m(:, n    , nn) = (aQ01 + aQ10) * (4 * p.k);
      %m(:, n + 3, nn) = (aI01 + aI10) * (4 * p.k);
    end

    % == frequency input is a little different...
    % convert from mirror position back to phase (2 * p.k),
    % and integrate phase to get frequency (i * f)
    m(:, n    , 7) = m(:, n    , 7) ./ (i * f(:) * 2 * p.k);
    m(:, n + 3, 7) = m(:, n + 3, 7) ./ (i * f(:) * 2 * p.k);

  end

  % save some info for later
  p.cr = pCR;
  p.p1 = pP1;
  p.p2 = pP2;

  p.e01 = e01;
  p.e10 = e10;
  p.e12 = e12;
  p.e21 = e21;
  p.phi = phi;

  % ==== Basis Transforms
  % mirror positions, with proper BS calibration
  r2 = 1/sqrt(2);
  mMIR = diag([1 1 1 1 1 r2 1]);

  % cannonical DOFs, averages RM, BS and ITM motions for MICH and PRC
  mDOF = [-2 -2 -1 -1  0;  2 -2  1 -1  0
           0  0  1  1  0;  0  0 -1  1  0
           0  0 -1 -1  0;  0  0  1  0  0
           0  0  0  0  4] / 4;

  % cannonical DOFs, uses only RM and BS motion for MICH and PRC
  mDRV = [-1 -1  0  0  0;  1 -1  0  0  0
           0  0  0  0  0;  0  0  0  0  0
           0  0 -1 -1  0;  0  0  1  0  0
           0  0  0  0  2] / 2;

  if( length(varargin) > 0 )
    if( ~isstr(varargin{1}) )
      error('Non-string input basis specification.');
    elseif( strcmp(varargin{1}, 'DOF') )
      m = inputTrans(m, mDOF);
    elseif( strcmp(varargin{1}, 'DRV') )
      m = inputTrans(m, mDRV);
    elseif( strcmp(varargin{1}, 'CM') )
      m = inputTrans(m, mMIR);			% transform to MIR basis
      varargin{2} = 'CM';			% add CM transform
    elseif( strcmp(varargin{1}, 'COMP') )
      % no transform
    else
      error(sprintf('Unknown input basis specification "%s".', varargin{1}));
    end
   else
     m = inputTrans(m, mMIR);
   end

  if( length(varargin) > 1 )
    if( ~isstr(varargin{2}) )
      error('Non-string CM basis specification.');
    elseif( strcmp(varargin{2}, 'CM') )
      % build CM matrix using dFreq to force REFL_I to zero
      nf = size(m, 3);
      mCM = zeros(size(m, 1), nf, nf - 1);
      for n = 1:nf - 1
        mCM(:,nf, n) = -m(:, 5, n) ./ m(:, 5, nf);
        mCM(:, n, n) = 1;
      end

      % transform to CM basis (no dFreq)
      m = inputTrans(m, mCM);
    else
      error(sprintf('Unknown CM basis specification "%s".', varargin{2}));
    end
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% freqProp
% compute frequency dependent IFO properties
%
% the param struct, "p", must contain amplitude
% reflectivities and transmissivities,
%  rRM, rBS, tBS, rIX, tIX, rIY, tIY, rEX, rEY
% as well as optical path lengths (in meters)
%  lPrmX, lPrmY, lArmX, lArmY
% and the speed of light, p.c, (in meters/second)
%
% The return struct contains a 13 x 13 transfer matrix
% for field propagation from point to point in the IFO.
% The 13 field locations are
%  1) ETMx - from ITMx
%  2) ETMy - from ITMy
%  3) ITMx - from ETMx
%  4) ITMx - from BS
%  5) ITMy - from ETMy
%  6) ITMy - from BS
%  7) RM   - from BS
%  8) RM   - from REFL port
%  9) BS   - from ITMx
% 10) BS   - from ITMy
% 11) BS   - from RM
% 12) REFL - output
% 13) AS   - output
%
% A few of the elements are picked out specifically
% as "cavity gains".  They are
%  gX = X-Arm Cavity Gain
%  gY = Y-Arm Cavity Gain
%  gP = Power Recycling Cavity Gain
%
% ps = freqProp(p, f)

function ps = freqProp(p, f)

  % ==== Phase Computations
  % == One-Way Phases
  dphi_dl = -2 * pi * i * f / p.c;

  ps.pArmX = exp(dphi_dl * p.lArmX);
  ps.pArmY = exp(dphi_dl * p.lArmY);
  ps.pRMbs = exp(dphi_dl * p.lRMbs);
  ps.pBSx = exp(dphi_dl * p.lBSx);
  ps.pBSy = exp(dphi_dl * p.lBSy);

  Nf = length(f);
  mc = zeros(13, 13, Nf);

  % via ETMx and ETMy
  mc(3, 1, :) = ps.pArmX * p.rEX;	% from ITMx to ITMx via ETMx
  mc(5, 2, :) = ps.pArmY * p.rEY;	% from ITMy to ITMy via ETMy

  % via ITMx
  mc(1, 3, :) = ps.pArmX * p.rIX;	% from ETMx to ETMx
  mc(9, 3, :) = ps.pBSx * p.tIX;	% from ETMx to BS
  mc(1, 4, :) = ps.pArmX * p.tIX;	% from BS to ETMx
  mc(9, 4, :) = ps.pBSx * p.bIX;	% from BS to BS

  % via ITMy
  mc( 2, 5, :) = ps.pArmY * p.rIY;	% from ETMy to ETMy front
  mc(10, 5, :) = ps.pBSy * p.tIY;	% from ETMy to BS
  mc( 2, 6, :) = ps.pArmY * p.tIY;	% from BS to ETMy front
  mc(10, 6, :) = ps.pBSy * p.bIY;	% from BS to BS

  % via RM
  mc(11, 7, :) = p.rRM * ps.pRMbs;	% from BS to BS
  mc(12, 7, :) = p.tRM;			% from BS to REFL
  mc(11, 8, :) = p.tRM * ps.pRMbs;	% from REFL to BS
  mc(12, 8, :) = p.bRM;			% from REFL to REFL

  % via BS
  mc( 7,  9, :) = p.tBS * ps.pRMbs;	% from ITMx to RM
  mc(13,  9, :) = p.bBS;		% from ITMx to AS
  mc( 7, 10, :) = p.rBS * ps.pRMbs;	% from ITMy to RM
  mc(13, 10, :) = p.tBS;		% from ITMy to AS
  mc( 4, 11, :) = p.tBS * ps.pBSx;	% from RM to ITMx
  mc( 6, 11, :) = p.rBS * ps.pBSy;	% from RM to ITMy

  % ==== Invert Propagation Matrix
  m = invprop(mc);
  ps.mc = mc;
  ps.m = m;

  % cavity gains
  ps.gX = m(1, 1, :) ./ (m(4, 1, :) .* mc(1, 4) + 1);
  ps.gY = m(2, 2, :) ./ (m(6, 2, :) .* mc(2, 6) + 1);
  ps.gP = m(7, 7, :);

  % composite reflectivites
  ps.rX = m( 9, 8, :) ./ m( 4, 8, :);
  ps.rY = m(10, 8, :) ./ m( 6, 8, :);
  ps.rM = m( 7, 8, :) ./ m(11, 8, :);
  ps.tM = m(13, 8, :) ./ m(11, 8, :);
  ps.rIFO = m(12, 8, :) ./ m(8, 8, :);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fieldAmp
% compute field amplitudes
%
% eAS and ePOB are in the cavity mode
% eREFL(1) - prompt reflection in the input mode
% eREFL(2) - leakage field in the cavity mode
% eREFL(3) - reject field in the (input - cavity) mode

function ps = fieldAmp(p, ps, eIn)

  if( size(ps.m, 3) ~= 1 )
    error('fieldAmp only accepts single frequency field matrices.')
  end

  ps.eIFO = ps.m(:, 8) * ps.aIn * eIn;

  % AS and POB
  ps.eAS = ps.eIFO(13);
  ps.ePOB = ps.eIFO(11);

  % REFL
  ps.eREFL = zeros(3, 1);
  ps.eREFL(1) = ps.mc(12, 8) * eIn;		% prompt reflection at REFL
  ps.eREFL(2) = ps.mc(12, 7) * ps.eIFO(7);	% leakage field to REFL
  ps.eREFL(3) = -(p.tRM^2 / (1 + p.rRM)) * ...	% anti-resonant
                 sqrt(1 - ps.aIn^2) * eIn;	% non-modematched


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sigMatrix
% compute signal matrix for some frequency vector
%
% ps is the source field struct (from freqProp, with a single frequency)
% pf is the source field struct (from freqProp, with a frequency vector)
%
% Each matrix element is computed with a sum of terms of the form:
% m  = (propagation throughout IFO) .* (reflection to next field point) *
%      (source field at signal gerating reflector);
%
% This equation should be read from right to left (i.e., source, then
% reflection, then propagation).  Reflection from the "back" of an
% optic is given a negative sign in the summation.
%
% the return matrix is 3 dimentional with
% N frequencies
% 3 ports [AS; REFL; POB]
% 7 sources [ETMx; ETMy; ITMx; ITMy; RM; BS; dFreq]

function ms = sigMatrix(ps, pf)

  % init
  eIFO = ps.eIFO;				% source field vector
  mr = permute(pf.mc, [3 1 2]);			% reflection matrix
  mp = permute(pf.m, [3 1 2]);			% propagation matrix

  np = 13:-1:11;				% port indices, [AS; REFL; POB]
  ms = zeros(size(mp, 1), length(np), 7);	% signal matrix
  
  % compute signal matrix
  for n = 1:length(np)
    p = np(n);
    ms(:, n, 1) = mp(:, p, 3) .* mr(:, 3, 1) * eIFO(1);		% ETMx front
    ms(:, n, 2) = mp(:, p, 5) .* mr(:, 5, 2) * eIFO(2);		% ETMy front
    ms(:, n, 3) = mp(:, p, 1) .* mr(:, 1, 3) * eIFO(3) - ...	% ITMx front
                  mp(:, p, 9) .* mr(:, 9, 4) * eIFO(4);		% ITMx back
    ms(:, n, 4) = mp(:, p, 2) .* mr(:, 2, 5) * eIFO(5) - ...	% ITMy front
                  mp(:, p,10) .* mr(:,10, 6) * eIFO(6);		% ITMy back
    ms(:, n, 5) = mp(:, p,11) .* mr(:,11, 7) * eIFO(7) - ...	% RM front
                  mp(:, p,12) .* mr(:,12, 8) * eIFO(8);		% RM back
    ms(:, n, 6) =-mp(:, p,13) .* mr(:,13, 9) * eIFO(9) + ...	% BS from ITMx
                  mp(:, p, 7) .* mr(:, 7,10) * eIFO(10) + ...	% BS from ITMy
                  mp(:, p, 6) .* mr(:, 6,11) * eIFO(11);	% BS from RM
    ms(:, n, 7) = mp(:, p, 8) * eIFO(8);			% dFreq
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% demod
% demodulation signal between two fields
%
% aSig is the singal carrying field
%  aSig(end) is used as the RF phase reference
%  and must be generated with an audio frequency
%  of 0 (i.e., f(end) = 0 in call to freqProp).
%
% aLO is the local-oscillator
%
% df is sign of the demodulation
%  df = (fLO - fSig) / fDemod, where fLO
%  is the frequency of the local oscillator
%  field, fSig is the frequency of the signal
%  carrier (this does not include the audio
%  signal frequency), and fDemod is the
%  demodulation frequency.  If df is not -1 or
%  1, the demodulation results are zero.
%
% phiDemod is the demod phase in radians
%
% [aI, aQ] = demod(aSig, aLO, df, phi)

function [aI, aQ] = demod(aSig, aLO, df, phiDemod)

  % assuming f(end) == 0, phase of aSig(end) is RF phase,
  % remove this phase from aSig to leave audio signal
  aRF = exp(i * angle(aSig(end)));
  aAF = aSig(1:end-1) / aRF;

  % conjugate upper RF freq
  if( df == 1 )
    aLO = conj(aLO);
  elseif( df == -1 )
    aRF = conj(aRF);
    aAF = -aAF;
  else
    error('Third argument, df, must be -1 or 1.');
  end

  % make I and Q RF phases, then replace AF phase
  a = aLO * aRF * exp(i * phiDemod);
  aI = -real(a) * aAF;
  aQ =  imag(a) * aAF;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% inputTrans
% transform m into input basis mDOF

function mr = inputTrans(m, mDOF);

  if( ndims(mDOF) == 2 )
    mr = zeros(size(m, 1), size(m, 2), size(mDOF, 2));
    for n = 1:size(m, 1)
      mr(n, :, :) = shiftdim(shiftdim(m(n, :, :), 1) * mDOF, -1);
    end
  else
    mr = zeros(size(m, 1), size(m, 2), size(mDOF, 3));
    for n = 1:size(m, 1)
      mr(n, :, :) = shiftdim(shiftdim(m(n, :, :), 1) * ...
                          shiftdim(mDOF(n, :, :), 1), -1);
    end
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% invprop
% invert optical propagation matrix

function m = invprop(m)

  % consolidate warnings, if any
  warnState = warning;
  warning off
  lastwarn('') 

  % invert
  eyeN = eye(size(m, 2));
  for nf = 1:size(m, 3);
    m(:, :, nf) = inv(eyeN - m(:, :, nf));
  end

  % deal with warnings
  if( ~isempty(lastwarn) )
    warning('IFO propagation nearly singular.')
  end
  warning(warnState)
