% Fjinv   Inverse of 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/2.
%         Fj_inv(x,j) gives F_j^{-1}(x). The algorithm looks for the 
%         ground state 
%
%         H=Jz^2-lambdax*Jx;
%
%         It searches for the lambdax corresponding to the state 
%         with var(Jz)=xj. Then, it gives back the expectation value of Jx
%         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_inv(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 
%               (var(Jz)-var(Jz)_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_inv.m 

% Program by Geza Toth (C) 2013.

function exJz_over_j=Fj_inv(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=20000;
    % Treshold in the accuracy of x
    treshold=0.0000000001;
    % 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/2
    error('Implemented only for x<=1/2.');
end %if

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

% The variance of jz from the input
vaJz_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;

while flag
    
    % Obtain the spin squeezed state as a ground state
    H=Jz^2-lambdax*Jx;
    [v,d]=eig(H);
    DeltaH=0.00001;
    indices=find(real(diag(d))<min(real(diag(d)))+DeltaH);
    if length(indices)>1,
        disp('Degenerate eigenvalues!')
    end %if
    
    % Ground state
    groundstate=v(:,indices(1));
    
    % Expectation value of Jx
    %exJx=ex(Jx,groundstate);
    
    % Variance of Jz
    vaJz=va(Jz,groundstate);
    
    if  vaJz<vaJz_input,
        lambdaxL=lambdax;
        lambdax=(lambdaxH+lambdaxL)/2;
    else
        lambdaxH=lambdax;
        lambdax=(lambdaxH+lambdaxL)/2;
    end %if
    
    difference=vaJz-vaJz_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;
exJz_over_j=ex(Jx,groundstate)/j;