(* ::Package:: *)

BeginPackage["ParticleEntanglement`",{"SqueezeRoutines`"}]

(*

Eddy Ardonne and Nicolas Regnault

Squeeze Routines for polarized and spinful quantum Hall states,
bosonic as well as fermionic. Routines for the spinful case
are based on

The structure of spinful quantum Hall states: a squeezing perspective,
E. Ardonne, N. Regnault,
arXiv:1107.2232

See this paper for references for the polarized case, and for reference
about particle and orbital entanglement spectra.

Routines making use of the additional symmetries in case of L=0 and/or
S=0 states are not implemented in this version. A routine to calculate
only the number of solutions to the consistency condition via an
eigenvalue problem is not included either. Please contact E. Ardonne
about these (and possible other) issues.

*)


bosonpes::usage=
 "bosonpes[na] will calculate the particle entanglement spectrum with na particles in
  part A, for a (spinless) bosonic state."

fermionpes::usage=
 "fermionpes[na] will calculate the particle entanglement spectrum with na particles in
  part A, for a (spinless) fermionic state."

spinbosonpes::usage=
 "spinbosonpes[na] will calculate the particle entanglement spectrum with na particles in
  part A, for a spinful bosonic state."

spinfermionpes::usage=
 "spinfermionpes[na] will calculate the particle entanglement spectrum with na particles in
  part A, for a spinful fermionic state."

lmultiplets::usage=
 "lmultiplets[lzstates] gives the number of l multiplets, based on the number of
  lzstates (with lz values non-negagive). Note that no checks are done to see if
  total l is a good symmetry!"

lsmultiplets::usage=
 "lmultiplets[lzszstates] gives the number of l and multiplets, based on the number of
  lz and sz states (with lz and sz values non-negagive). Note that no checks are
  done to see if total l and s are good symmetries!"

lszmultiplets::usage=
 "lszmultiplets[lzszstates] gives the number of l multiplets for each positive sz,
  based on the number of lz ans sz states (with lz and sz values non-negagive). Note
  that no checks are done to see if total l is a good symmetry!"

datalsorted::usage=
 "Reduces the dataset data, such that only the highest lz values of l multiplets
  remain."

datalszsorted::usage=
 "Reduces the dataset data, such that only the highest lz values of l multiplets
  remain. Note that it is assumed that s is not a good quantum number."

datalssorted::usage=
 "Reduces the dataset data, such that only the highest lz and sz values of l and s
  multiplets remain. Degenerate multiplets will appear only once.
  No checks are done if l and s are good quantum numbers!"


Begin["`Private`"]

(** Some functions, needed to generate the root configurations, which
are subsequently used to generate the Hilbertspaces system A and B **)

lz[vec_]:=With[{norb=Length[vec]},Sum[(-(norb-1)/2-1+i)vec[[i]],{i,1,norb}]];

nfluxneeded[0,cellsize_,maxincell_]:=0;
nfluxneeded[Ne_/;Ne>0,cellsize_,maxincell_]:=cellsize( Ceiling[Ne/(maxincell)]-1)

firstconfig[Ne_,nflux_,cellsize_,maxincell_]:=
If[nflux<nfluxneeded[Ne,cellsize,maxincell],{},Table[If[Mod[i,cellsize]==0,Max[Min[maxincell,(Ne-i/cellsize  maxincell)],0],0],{i,0,nflux}]]

hsrootconfigs[Ne_,nflux_,cellsize_,maxincell_]:=
With[{fc=firstconfig[Ne,nflux,cellsize,maxincell]},
If[fc=={},{},
maxp=partition[Reverse[fc]];
res={fc};
While[
lz[res[[-1]]]<Abs[lz[fc]],
tobeincreased=Position[partition[res[[-1]]]-maxp,x_/;x<0,1,1][[1,1]];
res=Append[res,orbocc[ReplacePart[partition[res[[-1]]],tobeincreased->partition[res[[-1]]][[tobeincreased]]+1],nflux]];
];
res
]];

hsbosonrootconfigs[Ne_,nflux_]:=hsrootconfigs[Ne,nflux,1,Ne];
hsfermionrootconfigs[Ne_,nflux_]:=hsrootconfigs[Ne,nflux,1,1];
hsspinbosonrootconfigs[Ne_,nflux_]:=hsrootconfigs[Ne,nflux,1,Ne];
hsspinfermionrootconfigs[Ne_,nflux_]:=hsrootconfigs[Ne,nflux,1,2];

id[num_]:=integerdigits[num,Global`base,Global`norbitals]
fd[rp_]:=fromdigits[rp,Global`base]
sid[num_]:={id[num[[1]]],id[num[[2]]]}
sfd[rp_]:={fd[rp[[1]]],fd[rp[[2]]]}


disknorm[flux_,l_]:=1/Sqrt[2^l l!];
spherenorm[flux_,l_]:=Sqrt[Binomial[flux,l]];
cftnorm[flux_,l_]:=1;

Global`normalization="sphere";

(** Routines to generate the Hilbertspaces of the parts A and B of the system **)

(* Spinless bosons *)

bosonsqueezea[rp_]:=With[{},
norbitalsa=Length[rp];
nparticlesa=Plus@@rp;
base=Global`nparticles+1;
Clear[statea,statepositiona,statepresenta];
hsdima=1;
olddim=1;
currentdim=1;
statea[hsdima]=FromDigits[rp,base];
statepositiona[statea[hsdima]]=hsdima;
statepresenta[statea[hsdima]]=True;
new=True;
While[new,
new=False;
Do[currentpart=integerdigits[statea[i],base,norbitalsa];
If[currentpart[[i1]]>0&&currentpart[[j1]]>0,
newpart=Table[Piecewise[{{currentpart[[jj]]-1,jj==i1},{currentpart[[jj]]-1,jj==j1},{currentpart[[jj]]+2,jj==j1-1&&jj==i1+1},{currentpart[[jj]]+1,jj==i1+1},{currentpart[[jj]]+1,jj==j1-1}},currentpart[[jj]]],{jj,1,norbitalsa}];
newstate=FromDigits[newpart,base];
If[statepresenta[newstate],Null,Null,
new=True;hsdima++;statea[hsdima]=newstate;
statepositiona[newstate]=hsdima;statepresenta[newstate]=True;
];
];
,{i,olddim,currentdim},{j1,3,norbitalsa},{i1,1,j1-2}];
olddim=currentdim+1;currentdim=hsdima;
];
]

bosonsqueezeb[rp_]:=With[{},
norbitalsb=Length[rp];
nparticlesb=Plus@@rp;
base=Global`nparticles+1;
Clear[stateb,statepositionb,statepresentb];
hsdimb=1;
olddim=1;
currentdim=1;
stateb[hsdimb]=FromDigits[rp,base];
statepositionb[stateb[hsdimb]]=hsdimb;
statepresentb[stateb[hsdimb]]=True;
new=True;
While[new,
new=False;
Do[currentpart=integerdigits[stateb[i],base,norbitalsb];
If[currentpart[[i1]]>0&&currentpart[[j1]]>0,
newpart=Table[Piecewise[{{currentpart[[jj]]-1,jj==i1},{currentpart[[jj]]-1,jj==j1},{currentpart[[jj]]+2,jj==j1-1&&jj==i1+1},{currentpart[[jj]]+1,jj==i1+1},{currentpart[[jj]]+1,jj==j1-1}},currentpart[[jj]]],{jj,1,norbitalsb}];
newstate=FromDigits[newpart,base];
If[statepresentb[newstate],Null,Null,
new=True;hsdimb++;stateb[hsdimb]=newstate;
statepositionb[newstate]=hsdimb;statepresentb[newstate]=True;
];
];
,{i,olddim,currentdim},{j1,3,norbitalsb},{i1,1,j1-2}];
olddim=currentdim+1;currentdim=hsdimb;
];
]


(* Spinless fermions *)

fermionsqueezea[rp_]:=With[{},
norbitals=Length[rp];
base=2;
Clear[statea,statepositiona,statepresenta];
hsdima=1;
olddim=1;
currentdim=1;
statea[hsdima]=FromDigits[rp,base];
statepositiona[statea[hsdima]]=hsdima;
statepresenta[statea[hsdima]]=True;
new=True;
While[new,
new=False;
Do[currentpart=integerdigits[statea[i],base,norbitals];
If[currentpart[[i1]]>0&&currentpart[[j1]]>0&&currentpart[[i1+1]]==currentpart[[j1-1]]==0,
newpart=Table[Piecewise[{{currentpart[[jj]]-1,jj==i1},{currentpart[[jj]]-1,jj==j1},{currentpart[[jj]]+2,jj==j1-1&&jj==i1+1},{currentpart[[jj]]+1,jj==i1+1},{currentpart[[jj]]+1,jj==j1-1}},currentpart[[jj]]],{jj,1,norbitals}];
newstate=FromDigits[newpart,base];
If[statepresenta[newstate],Null,Null,
new=True;hsdima++;statea[hsdima]=newstate;
statepositiona[newstate]=hsdima;statepresenta[newstate]=True;
];
];
,{i,olddim,currentdim},{j1,4,norbitals},{i1,1,j1-3}];
olddim=currentdim+1;currentdim=hsdima;
];
]

fermionsqueezeb[rp_]:=With[{},
norbitals=Length[rp];
base=2;
Clear[stateb,statepositionb,statepresentb];
hsdimb=1;
olddim=1;
currentdim=1;
stateb[hsdimb]=FromDigits[rp,base];
statepositionb[stateb[hsdimb]]=hsdimb;
statepresentb[stateb[hsdimb]]=True;
new=True;
While[new,
new=False;
Do[currentpart=integerdigits[stateb[i],base,norbitals];
If[currentpart[[i1]]>0&&currentpart[[j1]]>0&&currentpart[[i1+1]]==currentpart[[j1-1]]==0,
newpart=Table[Piecewise[{{currentpart[[jj]]-1,jj==i1},{currentpart[[jj]]-1,jj==j1},{currentpart[[jj]]+2,jj==j1-1&&jj==i1+1},{currentpart[[jj]]+1,jj==i1+1},{currentpart[[jj]]+1,jj==j1-1}},currentpart[[jj]]],{jj,1,norbitals}];
newstate=FromDigits[newpart,base];
If[statepresentb[newstate],Null,Null,
new=True;hsdimb++;stateb[hsdimb]=newstate;
statepositionb[newstate]=hsdimb;statepresentb[newstate]=True;
];
];
,{i,olddim,currentdim},{j1,4,norbitals},{i1,1,j1-3}];
olddim=currentdim+1;currentdim=hsdimb;
];
]

(* Spinful bosons *)

spinbosonsqueezeintt[rp_]:=With[{},
norbitals=Length[rp];
base=Global`nparticles+1;
Clear[intstate,intstateposition,intstatepresent];
inthsdim=1;
olddim=1;
currentdim=1;
intstate[inthsdim]=FromDigits[rp,base];
intstateposition[intstate[inthsdim]]=inthsdim;
intstatepresent[intstate[inthsdim]]=True;
new=True;
While[new,
new=False;
Do[currentpart=integerdigits[intstate[i],base,norbitals];
If[currentpart[[i1]]>0&&currentpart[[j1]]>0,
newpart=Table[Piecewise[{{currentpart[[jj]]-1,jj==i1},{currentpart[[jj]]-1,jj==j1},{currentpart[[jj]]+2,jj==j1-1&&jj==i1+1},{currentpart[[jj]]+1,jj==i1+1},{currentpart[[jj]]+1,jj==j1-1}},currentpart[[jj]]],{jj,1,norbitals}];
newstate=FromDigits[newpart,base];
If[intstatepresent[newstate],Null,Null,
new=True;inthsdim++;intstate[inthsdim]=newstate;
intstateposition[newstate]=inthsdim;intstatepresent[newstate]=True;
];
];
,{i,olddim,currentdim},{j1,3,norbitals},{i1,1,j1-2}];
olddim=currentdim+1;currentdim=inthsdim;
];
]

addspinbosona[nu_,nd_]:=With[{},
Clear[statea,statepositiona,statepresenta];
nup=nu;
ndown=nd;
sd=spindist[nup,ndown];
hsdima=0;
Do[dist=integerdigits[intstate[i1],base,norbitals];orb=orbtopart[dist];Do[newpart=Catch[Module[{res=Table[0,{s,1,2},{i,1,Length[dist]}]},Do[res[[sd[[i,j]],orb[[j]]]]++,{j,1,nup+ndown}];Throw[res]]];
newstate=spinfromdigits[newpart,base];
If[statepresenta[newstate],Null,Null,hsdima++;statea[hsdima]=newstate;statepositiona[newstate]=hsdima;statepresenta[newstate]=True;]
,{i,1,Length[sd]}]
,{i1,1,inthsdim}]
]


addspinbosonb[nu_,nd_]:=With[{},
Clear[stateb,statepositionb,statepresentb];
nup=nu;
ndown=nd;
sd=spindist[nup,ndown];
hsdimb=0;
Do[dist=integerdigits[intstate[i1],base,norbitals];orb=orbtopart[dist];Do[newpart=Catch[Module[{res=Table[0,{s,1,2},{i,1,Length[dist]}]},Do[res[[sd[[i,j]],orb[[j]]]]++,{j,1,nup+ndown}];Throw[res]]];
newstate=spinfromdigits[newpart,base];
If[statepresentb[newstate],Null,Null,hsdimb++;stateb[hsdimb]=newstate;statepositionb[newstate]=hsdimb;statepresentb[newstate]=True;]
,{i,1,Length[sd]}]
,{i1,1,inthsdim}]
]

spinbosonsqueezea[rp_,nup_,ndown_]:=With[{},spinbosonsqueezeintt[rp];addspinbosona[nup,ndown]]

spinbosonsqueezeb[rp_,nup_,ndown_]:=With[{},spinbosonsqueezeintt[rp];addspinbosonb[nup,ndown]]


(* Spinful fermions *)

spinfermionsqueezeintt[rp_]:=With[{},
norbitals=Length[rp];
base=3;
Clear[intstate,intstateposition,intstatepresent];
inthsdim=1;
olddim=1;
currentdim=1;
intstate[inthsdim]=FromDigits[rp,base];
intstateposition[intstate[inthsdim]]=inthsdim;
intstatepresent[intstate[inthsdim]]=True;
new=True;
While[new,
new=False;
Do[currentpart=integerdigits[intstate[i],base,norbitals];
If[currentpart[[i1]]>0&&currentpart[[j1]]>0&&If[i1==j1-2,currentpart[[i1+1]]==0,currentpart[[i1+1]]<2&&currentpart[[j1-1]]<2],
newpart=Table[Piecewise[{{currentpart[[jj]]-1,jj==i1},{currentpart[[jj]]-1,jj==j1},{currentpart[[jj]]+2,jj==j1-1&&jj==i1+1},{currentpart[[jj]]+1,jj==i1+1},{currentpart[[jj]]+1,jj==j1-1}},currentpart[[jj]]],{jj,1,norbitals}];
newstate=FromDigits[newpart,base];
If[intstatepresent[newstate],Null,Null,
new=True;inthsdim++;intstate[inthsdim]=newstate;
intstateposition[newstate]=inthsdim;intstatepresent[newstate]=True;
];
];
,{i,olddim,currentdim},{j1,3,norbitals},{i1,1,j1-2}];
olddim=currentdim+1;currentdim=inthsdim;
];
]

addspinfermiona[nu_,nd_]:=With[{},
Clear[statea,statepositiona,statepresenta];
base=2;
nup=nu;
ndown=nd;
sd=spindist[nup,ndown];
hsdima=0;
Do[dist=integerdigits[intstate[i1],3,norbitals];orb=orbtopart[dist];Do[newpart=Catch[Module[{res=Table[0,{s,1,2},{i,1,Length[dist]}]},Do[res[[sd[[i,j]],orb[[j]]]]++,{j,1,nup+ndown}];Throw[res]]];
If[FreeQ[newpart,2],
newstate=spinfromdigits[newpart,base];
If[statepresenta[newstate],Null,Null,hsdima++;statea[hsdima]=newstate;statepositiona[newstate]=hsdima;statepresenta[newstate]=True;]]
,{i,1,Length[sd]}]
,{i1,1,inthsdim}]
]

addspinfermionb[nu_,nd_]:=With[{},
Clear[stateb,statepositionb,statepresentb];
base=2;
nup=nu;
ndown=nd;
sd=spindist[nup,ndown];
hsdimb=0;
Do[dist=integerdigits[intstate[i1],3,norbitals];orb=orbtopart[dist];Do[newpart=Catch[Module[{res=Table[0,{s,1,2},{i,1,Length[dist]}]},Do[res[[sd[[i,j]],orb[[j]]]]++,{j,1,nup+ndown}];Throw[res]]];
If[FreeQ[newpart,2],
newstate=spinfromdigits[newpart,base];
If[statepresentb[newstate],Null,Null,hsdimb++;stateb[hsdimb]=newstate;statepositionb[newstate]=hsdimb;statepresentb[newstate]=True;]]
,{i,1,Length[sd]}]
,{i1,1,inthsdim}]
]

spinfermionsqueezea[rp_,nup_,ndown_]:=With[{},spinfermionsqueezeintt[rp];addspinfermiona[nup,ndown]]

spinfermionsqueezeb[rp_,nup_,ndown_]:=With[{},spinfermionsqueezeintt[rp];addspinfermionb[nup,ndown]]





(** Routines to calculate the particle entanglement spectrum **)

(* Spinless bosons *)


bosonpes[na_]:=With[{},

(* Set the normalization *)
Piecewise[{
{onepartnorm[flux_,l_]:=spherenorm[flux,l];,Global`normalization=="sphere"},
{onepartnorm[flux_,l_]:=cftnorm[flux,l];,Global`normalization=="cft"},
{onepartnorm[flux_,l_]:=disknorm[flux,l];,Global`normalization=="disk"}}
,onepartnorm[flux_,l_]:=spherenorm[flux,l];];

nparticles=Global`nparticles;
norbitals=Global`norbitals;
hsdim=Global`hsdim;

(* Normalizing the state *)

hsnorm=Table[With[{vec=id[Global`state[i]]},1/Product[onepartnorm[norbitals-1,j-1]^(vec[[j]])//N,
{j,1,norbitals}]],{i,1,hsdim}];

deg=Table[With[{vec=id[Global`state[i]]},Sqrt[nparticles!/Product[id[Global`state[i]][[j]]!//N,{j,1,norbitals}]]],{i,1,hsdim}];

norm=Sqrt[Sum[(hsnorm[[i]]*deg[[i]]* Global`b[i])^2,{i,1,hsdim}]];

(* Initializing angular momentum varialbles *)
lt=(norbitals-1)/2;
lstate=Sum[(lt-j)id[Global`state[1]][[j+1]],{j,0,norbitals-1}];

(* The number of particles in part B *)
nb=nparticles-na;

rootconfigsa=hsbosonrootconfigs[na,norbitals-1];
rootconfigsb=hsbosonrootconfigs[nb,norbitals-1];
lzvaluesa=lz/@rootconfigsa;
lzvaluesb=lz/@rootconfigsb;

(* Calculating the reduced density matrices, and their singular values *)
Do[
apos=Position[lzvaluesa,lzval][[1,1]];
bpos=Position[lzvaluesb,lstate-lzval][[1,1]];
bosonsqueezea[rootconfigsa[[apos]]];
bosonsqueezeb[rootconfigsb[[bpos]]];
dega=Table[With[{vec=id[statea[i]]},Sqrt[na!/Product[id[statea[i]][[j]]!//N,{j,1,norbitals}]]],{i,1,hsdima}];
degb=Table[With[{vec=id[stateb[i]]},Sqrt[nb!/Product[id[stateb[i]][[j]]!//N,{j,1,norbitals}]]],{i,1,hsdimb}];
cmat[lzval]=Table[With[{hspos=Global`stateposition[statea[i]+stateb[j]]},If[NumberQ[hspos],Global`b[hspos] dega[[i]] degb[[j]]hsnorm[[hspos]]/norm,0]],{i,1,hsdima},{j,1,hsdimb}];
svd[lzval]=SingularValueList[cmat[lzval]//N]//Chop;
Global`svd[lzval]=svd[lzval];
,{lzval,lzvaluesa}];

checknorm=Plus@@Flatten[Table[#^2&/@svd[lzval],{lzval,lzvaluesa}]];
Print["Checking the norm gives: ",checknorm];

ententropy=-Plus@@Flatten[Table[2#^2 Log[#]&/@DeleteCases[Chop[svd[lzval]],0],{lzval,lzvaluesa}]];
Print["The entanglement entropy is: ",ententropy];

(* Collecting the data *)
Clear[data];
Global`lzpresent={};
Do[
If[DeleteCases[svd[lzval],0]=={},Null,
Global`lzpresent=Append[Global`lzpresent,lzval];
data[lzval]=Table[{lzval,-2Log[i]},{i,svd[lzval]}]
];
,{lzval,lzvaluesa}];

Global`data=Flatten[Table[data[i],{i,Global`lzpresent}],1];

lzdatarange=Cases[Global`lzpresent,x_/;x>=0];
Global`lzstates=Table[Count[Global`data,x_/;x[[1]]==l],{l,lzdatarange}];

];





(* Spinless fermions *)

fermionpes[na_]:=With[{},
(* Set the normalization *)
Piecewise[{
{onepartnorm[flux_,l_]:=spherenorm[flux,l];,Global`normalization=="sphere"},
{onepartnorm[flux_,l_]:=cftnorm[flux,l];,Global`normalization=="cft"},
{onepartnorm[flux_,l_]:=disknorm[flux,l];,Global`normalization=="disk"}}
,onepartnorm[flux_,l_]:=spherenorm[flux,l];];

nparticles=Global`nparticles;
norbitals=Global`norbitals;
hsdim=Global`hsdim;

(* Normalizing the state *)

hsnorm=Table[With[{vec=id[Global`state[i]]},1/Product[onepartnorm[norbitals-1,j-1]^(vec[[j]])//N,
{j,1,norbitals}]],{i,1,hsdim}];

norm=Sqrt[nparticles! Sum[(hsnorm[[i]]* Global`b[i])^2,{i,1,hsdim}]];

(* Initializing angular momentum varialbles *)
lt=(norbitals-1)/2;
lstate=Sum[(lt-j)id[Global`state[1]][[j+1]],{j,0,norbitals-1}];

(* The number of particles in part B *)
nb=nparticles-na;

rootconfigsa=hsfermionrootconfigs[na,norbitals-1];
rootconfigsb=hsfermionrootconfigs[nb,norbitals-1];
lzvaluesa=lz/@rootconfigsa;
lzvaluesb=lz/@rootconfigsb;

(* Calculating the reduced density matrices, and their singular values *)
Do[
apos=Position[lzvaluesa,lzval][[1,1]];
bpos=Position[lzvaluesb,lstate-lzval][[1,1]];
fermionsqueezea[rootconfigsa[[apos]]];
fermionsqueezeb[rootconfigsb[[bpos]]];
cmat[lzval]=Table[With[
{hspos=Global`stateposition[statea[i]+stateb[j]],
occapos=Flatten[Position[id[statea[i]],1]],
occbpos=Flatten[Position[id[stateb[j]],1]]},
If[NumberQ[hspos],Signature[Flatten[Append[occapos,occbpos]]]*Global`b[hspos]*Sqrt[na! * nb!]* hsnorm[[hspos]]/norm,0]],{i,1,hsdima},{j,1,hsdimb}];

svd[lzval]=SingularValueList[cmat[lzval]//N]//Chop;
Global`svd[lzval]=svd[lzval];
,{lzval,lzvaluesa}];

checknorm=Plus@@Flatten[Table[#^2&/@svd[lzval],{lzval,lzvaluesa}]];
Print["Checking the norm gives: ",checknorm];

ententropy=-Plus@@Flatten[Table[2#^2 Log[#]&/@DeleteCases[Chop[svd[lzval]],0],{lzval,lzvaluesa}]];
Print["The entanglement entropy is: ",ententropy];

(* Collecting the data *)
Clear[data];
Global`lzpresent={};
Do[
If[DeleteCases[svd[lzval],0]=={},Null,
Global`lzpresent=Append[Global`lzpresent,lzval];
data[lzval]=Table[{lzval,-2Log[i]},{i,svd[lzval]}]
];
,{lzval,lzvaluesa}];

Global`data=Flatten[Table[data[i],{i,Global`lzpresent}],1];

lzdatarange=Cases[Global`lzpresent,x_/;x>=0];
Global`lzstates=Table[Count[Global`data,x_/;x[[1]]==l],{l,lzdatarange}];


];



(* Spinful bosons *)

spinbosonpes[na_]:=With[{},
(* Set the normalization *)
Piecewise[{
{onepartnorm[flux_,l_]:=spherenorm[flux,l];,Global`normalization=="sphere"},
{onepartnorm[flux_,l_]:=cftnorm[flux,l];,Global`normalization=="cft"},
{onepartnorm[flux_,l_]:=disknorm[flux,l];,Global`normalization=="disk"}}
,onepartnorm[flux_,l_]:=spherenorm[flux,l];];

nparticles=Global`nparticles;
norbitals=Global`norbitals;
hsdim=Global`hsdim;

(* Normalizing the state *)

hsnorm=Table[With[{vec=sid[Global`state[i]]},1/Product[onepartnorm[norbitals-1,j-1]^(vec[[1,j]]+vec[[2,j]])//N,
{j,1,norbitals}]],{i,1,hsdim}];

deg=Table[With[{vec=sid[Global`state[i]]},Sqrt[nparticles!/Product[vec[[1,j]]!*vec[[2,j]]!,{j,1,norbitals}]]//N],{i,1,hsdim}];

norm=Sqrt[Sum[(hsnorm[[i]]*deg[[i]]* Global`b[i])^2,{i,1,hsdim}]];

(* Initializing angular momentum varialbles *)
lt=(norbitals-1)/2;
lstate=Sum[(lt-j)*(sid[Global`state[1]][[1,j+1]]+sid[Global`state[1]][[2,j+1]]),{j,0,norbitals-1}];
sstate=(Plus@@(sid[Global`state[1]][[1]]))-(Plus@@(sid[Global`state[1]][[2]]));


(* The number of particles in part B *)
nb=nparticles-na;

rootconfigsa=hsspinbosonrootconfigs[na,norbitals-1];
rootconfigsb=hsspinbosonrootconfigs[nb,norbitals-1];
lzvaluesa=lz/@rootconfigsa;
lzvaluesb=lz/@rootconfigsb;
szvaluesa=Table[i,{i,-na/2,na/2}];
szvaluesb=Table[i,{i,-nb/2,nb/2}];

szlzpresent={};
szpresent={};
lzpresent={};


(* Calculating the reduced density matrices, and their singular values *)
Do[
apos=Position[lzvaluesa,lzval][[1,1]];
bpos=Position[lzvaluesb,lstate-lzval][[1,1]];


spinbosonsqueezea[rootconfigsa[[apos]],na/2-szval,na/2+szval];
spinbosonsqueezeb[rootconfigsb[[bpos]],nb/2+szval+sstate/2,nb/2-szval-sstate/2];

If[
hsdima>0&&hsdimb>0,
dega=Table[With[{vec=sid[statea[i]]},Sqrt[na!/Product[vec[[1,j]]!vec[[2,j]]!//N,{j,1,norbitals}]]],{i,1,hsdima}];
degb=Table[With[{vec=sid[stateb[i]]},Sqrt[nb!/Product[vec[[1,j]]!vec[[2,j]]!//N,{j,1,norbitals}]]],{i,1,hsdimb}];
szlzpresent=Append[szlzpresent,{szval,lzval}];
szpresent=Union[Append[szpresent,szval]];
lzpresent=Union[Append[lzpresent,lzval]];
cmat[{szval,lzval}]=Table[With[{hspos=Global`stateposition[statea[i]+stateb[j]]},If[NumberQ[hspos],Global`b[hspos]* dega[[i]]* degb[[j]]*hsnorm[[hspos]]/norm,0]],{i,1,hsdima},{j,1,hsdimb}];
svd[{szval,lzval}]=SingularValueList[cmat[{szval,lzval}]//N]//Chop;
Global`svd[{szval,lzval}]=svd[{szval,lzval}];
]
,{lzval,lzvaluesa},{szval,szvaluesa}];

checknorm=Plus@@Flatten[Table[#^2&/@svd[szlzval],{szlzval,szlzpresent}]];
Print["Checking the norm gives: ",checknorm];

ententropy=-Plus@@Flatten[Table[2 * #^2 Log[#]& /@DeleteCases[Chop[svd[szlzval]],0],{szlzval,szlzpresent}]];
Print["The entanglement entropy is: ",ententropy];

(* Collecting the data *)
Clear[data];
Global`szlzpresent={};
Do[
If[DeleteCases[svd[{szval,lzval}],0]=={},Null,
Global`szlzpresent=Append[Global`szlzpresent,{szval,lzval}];
data[{szval,lzval}]=Table[{lzval,-2Log[i]},{i,svd[{szval,lzval}]}]
]
,{szval,szpresent},{lzval,Cases[szlzpresent,x_/;x[[1]]==szval][[All,2]]}];

szrange=Cases[szpresent,x_/;x>=0];
Global`data=Table[Flatten[Table[data[{szval,i}],{i,Cases[Global`szlzpresent,x_/;x[[1]]==szval][[All,2]]}],1],{szval,szrange}];

szdatarange=Cases[Union[Global`szlzpresent[[All,1]]],x_/;x>=0];
lzdatarange=Cases[Union[Global`szlzpresent[[All,2]]],x_/;x>=0];

Global`lzszstates=Table[Count[Global`data[[s]],x_/;x[[1]]==l],{s,1,Length[Global`data]},{l,lzdatarange}];

];

(* Spinful fermions *)

spinfermionpes[na_]:=With[{},
(* Set the normalization *)
Piecewise[{
{onepartnorm[flux_,l_]:=spherenorm[flux,l];,Global`normalization=="sphere"},
{onepartnorm[flux_,l_]:=cftnorm[flux,l];,Global`normalization=="cft"},
{onepartnorm[flux_,l_]:=disknorm[flux,l];,Global`normalization=="disk"}}
,onepartnorm[flux_,l_]:=spherenorm[flux,l];];

nparticles=Global`nparticles;
norbitals=Global`norbitals;
hsdim=Global`hsdim;

(* Normalizing the state *)

hsnorm=Table[With[{vec=sid[Global`state[i]]},1/Product[onepartnorm[norbitals-1,j-1]^(vec[[1,j]]+vec[[2,j]])//N,
{j,1,norbitals}]],{i,1,hsdim}];


norm=Sqrt[nparticles! * Sum[(hsnorm[[i]]* Global`b[i])^2,{i,1,hsdim}]];

(* Initializing angular momentum varialbles *)
lt=(norbitals-1)/2;
lstate=Sum[(lt-j)*(sid[Global`state[1]][[1,j+1]]+sid[Global`state[1]][[2,j+1]]),{j,0,norbitals-1}];
sstate=(Plus@@(sid[Global`state[1]][[1]]))-(Plus@@(sid[Global`state[1]][[2]]));


(* The number of particles in part B *)
nb=nparticles-na;

rootconfigsa=hsspinfermionrootconfigs[na,norbitals-1];
rootconfigsb=hsspinfermionrootconfigs[nb,norbitals-1];
lzvaluesa=lz/@rootconfigsa;
lzvaluesb=lz/@rootconfigsb;
szvaluesa=Table[i,{i,-na/2,na/2}];
szvaluesb=Table[i,{i,-nb/2,nb/2}];

szlzpresent={};
szpresent={};
lzpresent={};


(* Calculating the reduced density matrices, and their singular values *)
Do[
apos=Position[lzvaluesa,lzval][[1,1]];
bpos=Position[lzvaluesb,lstate-lzval][[1,1]];


spinfermionsqueezea[rootconfigsa[[apos]],na/2-szval,na/2+szval];
spinfermionsqueezeb[rootconfigsb[[bpos]],nb/2+szval+sstate/2,nb/2-szval-sstate/2];

If[
hsdima>0&&hsdimb>0,
szlzpresent=Append[szlzpresent,{szval,lzval}];
szpresent=Union[Append[szpresent,szval]];
lzpresent=Union[Append[lzpresent,lzval]];
cmat[{szval,lzval}]=Table[With[
{hspos=Global`stateposition[statea[i]+stateb[j]],
occaposup=Flatten[Position[sid[statea[i]][[1]],1]],
occbposup=Flatten[Position[sid[stateb[j]][[1]],1]],
occaposdown=Flatten[Position[sid[statea[i]][[2]],1]],
occbposdown=Flatten[Position[sid[stateb[j]][[2]],1]]}
,If[NumberQ[hspos],Global`b[hspos]*Signature[Flatten[Append[occaposup,occbposup]]]*Signature[Flatten[Append[occaposdown,occbposdown]]]*Sqrt[na! * nb!]*hsnorm[[hspos]]/norm,0]],{i,1,hsdima},{j,1,hsdimb}];
svd[{szval,lzval}]=SingularValueList[cmat[{szval,lzval}]//N]//Chop;
Global`svd[{szval,lzval}]=svd[{szval,lzval}];
]
,{lzval,lzvaluesa},{szval,szvaluesa}];

checknorm=Plus@@Flatten[Table[#^2&/@svd[szlzval],{szlzval,szlzpresent}]];
Print["Checking the norm gives: ",checknorm];

ententropy=-Plus@@Flatten[Table[2 * #^2 Log[#]& /@DeleteCases[Chop[svd[szlzval]],0],{szlzval,szlzpresent}]];
Print["The entanglement entropy is: ",ententropy];

(* Collecting the data *)
Clear[data];
Global`szlzpresent={};
Do[
If[DeleteCases[svd[{szval,lzval}],0]=={},Null,
Global`szlzpresent=Append[Global`szlzpresent,{szval,lzval}];
data[{szval,lzval}]=Table[{lzval,-2Log[i]},{i,svd[{szval,lzval}]}]
]
,{szval,szpresent},{lzval,Cases[szlzpresent,x_/;x[[1]]==szval][[All,2]]}];

szrange=Cases[szpresent,x_/;x>=0];
Global`data=Table[Flatten[Table[data[{szval,i}],{i,Cases[Global`szlzpresent,x_/;x[[1]]==szval][[All,2]]}],1],{szval,szrange}];


szdatarange=Cases[Union[Global`szlzpresent[[All,1]]],x_/;x>=0];
lzdatarange=Cases[Union[Global`szlzpresent[[All,2]]],x_/;x>=0];

Global`lzszstates=Table[Count[Global`data[[s]],x_/;x[[1]]==l],{s,1,Length[Global`data]},{l,lzdatarange}];


];


(**  Some routines, to extract the l and s (or sz) multiplet from the
number of lz and sz multiplets, as stored in the matrix lzszstates, which only
contains the number of lz & sz states, with non-negative lz and sz !! **)

lmultiplets[vec_]:=With[{dim=Length[vec]},
Table[If[l<dim,vec[[l]]-vec[[l+1]],vec[[l]]],{l,1,dim}]];


lszmultiplets[mat_]:=With[{dims=Dimensions[mat]},
Table[If[l<dims[[2]],mat[[s,l]]-mat[[s,l+1]],mat[[s,l]]],{s,1,dims[[1]]},{l,1,dims[[2]]}]];

lsmultiplets[mat_]:=With[{dims=Dimensions[mat]},
Table[
Piecewise[{
{mat[[s,l]]-mat[[s,l+1]]-(mat[[s+1,l]]-mat[[s+1,l+1]]),l<dims[[2]]&&s<dims[[1]]},
{mat[[s,l]]-mat[[s,l+1]],l<dims[[2]]&&s==dims[[1]]},
{mat[[s,l]]-mat[[s+1,l]],l==dims[[2]]&&s<dims[[1]]},
{mat[[s,l]],l==dims[[2]]&&s==dims[[1]]}}]
,{s,1,dims[[1]]},{l,1,dims[[2]]}]];


(* Some routines to sort the data according to the highest lz (and sz) values of
l (and s) multiplets *)

datalsorted:=With[{rdata=Table[Round[Global`data[[i,2]],10^(-10)],{i,1,Length[Global`data]}],rdatau=Union[Table[Round[Global`data[[i,2]],10^(-10)],{i,1,Length[Global`data]}]]},
Table[{Max[Global`data[[Flatten[Position[rdata,rdatau[[j]]]],1]]],rdatau[[j]]//N},{j,1,Length[rdatau]}]
];

datalszsorted:=Table[With[{rdata=Table[Round[Global`data[[sz,i,2]],10^(-10)],{i,1,Length[Global`data[[sz]]]}],rdatau=Union[Table[Round[Global`data[[sz,i,2]],10^(-10)],{i,1,Length[Global`data[[sz]]]}]]},
Table[{Max[Global`data[[sz,Flatten[Position[rdata,rdatau[[j]]]],1]]],rdatau[[j]]//N},{j,1,Length[rdatau]}]
],{sz,1,Length[szdatarange]}];

datalssorted:=Module[{datalist,rdatalist,rdata,rdatau,datalistls},
datalist=Flatten[Table[{Global`data[[si,j,1]],si,Global`data[[si,j,2]]},{si,1,Length[szdatarange]},{j,1,Length[Global`data[[si]]]}],1];
rdatalist=Table[{datalist[[i,1]],datalist[[i,2]],Round[datalist[[i,3]],10^(-10)]},{i,1,Length[datalist]}];
rdata=rdatalist[[All,3]];
rdatau=Union[rdata];
datalistls=Table[{Max[rdatalist[[Flatten[Position[rdatalist,rdatau[[i]]]],1]]],Max[rdatalist[[Flatten[Position[rdatalist,rdatau[[i]]]],2]]],rdatau[[i]]//N},{i,1,Length[rdatau]}];
Table[Cases[datalistls,x_/;x[[2]]==si][[All,{1,3}]],{si,1,Length[szdatarange]}]
]



End[ ]

EndPackage[ ]

