function [ObjectPlane_x,ObjectPlane_y,Magnification_x,Magnification_y,MagMat_x,MagMat_y]=FACET_EspecImaging(B_Vector,Gamma_Beam,ImagePlaneLocation)
% function [ObjectPlane_x,ObjectPlane_y,Magnification_x,Magnification_y,MagMat_x,MagMat_y]=FACET_EspecImaging(B_Vector,Gamma_Beam,ImagePlaneLocation)
% Script gives imaginng condition in x and y and corresponding matrizes and
% magnification between two z locations.


% Storage of different locations around the beamline. 
% These are not used in the code and is merely a lookup-table
    
    z_Filament =1993.2737;      %m z location of gas jet
    z_Q0D      =1996.98249;     %m  z location of Q0 (center)
    z_Q1D      =1999.206615;    %m z location of Q1 (center)
    z_Q2D      =2001.431049;    %m z location of Q2 (center)
    z_B5D36    =2005.940051;    %m Forgot what item what was... tbd.
    z_DTOTR    =2015.2599;      %m z location  of DT OTR screen
    z_LFOV     =2015.63;        %m z location of LFOV
    z_PRDump   =2017.53;        %m z location of PR dump
    z_PIC_CENT =1992.82;        %m z location of E320 IP 
% ------------------------------------------------------
    
% -------Define constants ----
    qE=1.6022e-19;
    c0=299792458;
    mE=9.1094e-31;
    kGauss2Telsa=10^-1;
% --------------------------------

    %Quad effective lengths
    L_Q0=1; %m
    L_Q1=1; %m
    L_Q2=1; %m
    L_Quad_Vector=[L_Q0,L_Q1,L_Q2];
%     ------------------------

%----Drift lengths----------
    LD2=(z_Q1D-L_Q1/2)-(z_Q0D+L_Q0/2); % Q0 to Q1
    LD3=(z_Q2D-L_Q2/2)-(z_Q1D+L_Q1/2); % Q2 to Q1
    LD4=ImagePlaneLocation-(z_Q2D+L_Q2/2); % Q2 to screen

    
% ---- Calculate  K -----------
    B_L_Vector=B_Vector*kGauss2Telsa./L_Quad_Vector; % field gradient in T/m
    K_Vector=B_L_Vector*qE/(Gamma_Beam*mE*c0);
% -----------------------------


% ---- Calculate beam transport matrices ------
    M2=Rmat(L_Q0,K_Vector(1)); % Q0
    M3=Rmat(LD2); % Q0 to Q1
    M4=Rmat(L_Q1,K_Vector(2)); % Q1
    M5=Rmat(LD3); % Q1 to Q2
    M6=Rmat(L_Q2,K_Vector(3)); % Q2
    M7=Rmat(LD4); % Q2 to screen
%   Mtot=M7*M6*M5*M4*M3*M2*M1; % This is the total matrix of the espec.
%   Just for reference, not needed.
    MatEspecPlusImageDrift=M7*M6*M5*M4*M3*M2;
%  ---- Calulate output --------
    ObjectPlane_x=-MatEspecPlusImageDrift(1,2)./MatEspecPlusImageDrift(1,1);
    MagMat_x=MatEspecPlusImageDrift*Rmat(ObjectPlane_x);


    ObjectPlane_y=-MatEspecPlusImageDrift(3,4)./MatEspecPlusImageDrift(3,3);
    MagMat_y=MatEspecPlusImageDrift.*Rmat(ObjectPlane_y);

    Magnification_x=abs(MagMat_x(1,1));
    Magnification_y=abs(MagMat_y(3,3));



end


function [matrix] = Rmat(l,k)
    %% GENERAL TRANSFER MATRIX (DRIFT AND QUADRUPOLE)
%     [matrix] = Rmat(l,k)
    % drift matrix (if k==0) or quad matrix
    if ~exist('k','var') || k == 0
        matrix = [1,l,0,0;0,1,0,0;0,0,1,l;0,0,0,1];
    else
        matrix = real([cos(sqrt(k)*l),sin(sqrt(k)*l)/sqrt(k),0,0;-sin(sqrt(k)*l)*sqrt(k),cos(sqrt(k)*l),0,0;0,0,cosh(sqrt(k)*l),sinh(sqrt(k)*l)/sqrt(k);0,0,sinh(sqrt(k)*l)*sqrt(k),cosh(sqrt(k)*l)]);
    end

end