% Fj   The function Fj(x) appearing in extreme spin squeezing 
%      of Sorensen and Molmer.
%      See: http://arxiv.org/abs/quant-ph/0011035
%      Only for integer j and x<1.
%      Fj(x,j) gives F_j(x). The algorithm looks for the ground state 
%
%         H=Jz^2-lambdax*Jx;
%
%      It searches for the lambdax corresponding to the state 
%      with ex(Jx)=xj. Then, it gives back var(Jz) for that state.
%
%      An addtional argument can be given
%      to set the important constants for the numerical algorithm.
%      (If not given, default values ae taken.)
%      In this case, the form is
%
%      Fj(x,j,[lambdaxL,lambdaxH,treshold,lambdatreshold])
%
%      lambdaxL=lower starting value for the Lagrange multiplyier lambdax
%      lambdaxH=higer starting value for the Lagrange multiplyier lambdax
%      treshold=the iteration continues until 
%               (ex(Jx)-ex(Jx)_input)<treshold j
%      lambdatreshold= search ends with error if
%               (lambdaxH-lambdaxL)<lambdatreshold
%      (In this case the lambdax we need is outside of the interval.)
%      For the default value, see the file.
%
%      See example_Fj.m 

% Program by Geza Toth (C) 2013.

function varJz_over_j=Fj(x,j,varargin)


if length(varargin)==1,
    M=varargin{1};
    lambdaxL=M(1);
    lambdaxH=M(2);
    treshold=M(3);
    lambdatreshold=M(4);
    
elseif length(varargin)==0,
    lambdaxL=0.00000001;
    lambdaxH=200000;
    % Treshold in the accuracy of x
    treshold=0.00000000001;
    % Treshold to test whether lamnbaxL amd lambadxH are almost equal
    % This is the case if the point is outside of the interval
    lambdatreshold=1e-12;
else
    error('Only 2 or 3 arguments are allowed for Fj.')
end %if
if x>1
    error('Implemented only for x<=1.');
end %if

% The expectation value of jx from the input
exJx_input=x*j;

% k-producabilty (k must be even!)
k=2*j;
%if ~even(k)
if mod(k,2)>0
    error('Implemented only for integer j.');
end %if

% Dimension of matrices describing Jx, Jy and Jz
dd=k+1;
[Jx,Jy,Jz]=su2(dd);

% Initial lambdax
lambdax=(lambdaxH+lambdaxL)/2;

flag=1;

% Treshold to detect degeneracy
TresholdH=0.0001;

while flag
    
    % Obtain the spin squeezed state as a ground state
    H=Jz^2-lambdax*Jx;
    [v,d]=eig(H);
    indices=find(real(diag(d))<min(real(diag(d)))+TresholdH);
    if length(indices)>1,
        disp('Degenerate eigenvalues!')
    end %if
    
    % Ground state
    groundstate=v(:,indices(1));
    
    % Expectation value of Jx
    exJx=ex(Jx,groundstate);
    
    if  exJx<exJx_input,
        lambdaxL=lambdax;
        lambdax=(lambdaxH+lambdaxL)/2;
    else
        lambdaxH=lambdax;
        lambdax=(lambdaxH+lambdaxL)/2;
    end %if
    
    difference=exJx-exJx_input;
    
    flag=abs(difference)>treshold*j;
    
    if abs(lambdaxL-lambdaxH)<lambdatreshold
       lambdaxH
       lambdaxL
       error('Point outside of the interval. Change lambdaxL and lambdaxH.') 
    end %if
    
end %if

varJz_over_j=va(Jz,groundstate)/j;