(* ::Package:: *)

(* ::Title:: *)
(*Different photon statistics for entangled and separable states*)


(* ::Subtitle:: *)
(*Video that demonstrates measurements of photon pairs in different bases of polarization*)


(* ::Item:: *)
(*Created in version: "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)"*)


(* ::Chapter:: *)
(*Photon statistics*)


(* ::Subchapter::Closed:: *)
(*Theoretical background (not part of the rest of the code)*)


(* ::Input:: *)
(*u = RotationMatrix[\[Theta]];*)
(*uu = KroneckerProduct[u, u];*)


(* ::Input:: *)
(*stateEnt={{0,0,0,0},{0,1/2,-(1/2),0},{0,-(1/2),1/2,0},{0,0,0,0}};*)
(*stateSep={{0,0,0,0},{0,1/2,0,0},{0,0,1/2,0},{0,0,0,0}};*)


(* ::Input:: *)
(*stateEntRot = ComplexExpand[uu.stateEnt.uu\[ConjugateTranspose]]//Simplify;*)
(*stateSepRot = ComplexExpand[uu.stateSep.uu\[ConjugateTranspose]]//Simplify;*)


(* ::Input:: *)
(*Diagonal/@{stateEntRot,stateSepRot}*)


(* ::Input:: *)
(*ClearAll[plotProbs]*)
(*plotProbs[probFun_,title_]:=Plot[Evaluate@probFun[\[Theta]],{\[Theta],0,2\[Pi]},PlotLabels->(Subscript["p",Row[{#1,#2}/.{0->"H",1->"V"}]]&@@@{{0,1},{0,0},{1,1},{1,0}}),PlotRange->{All,{0,1}},Ticks->{\[Pi]/2 Range[0,4],All},PlotLabel->title]*)


(* ::Input:: *)
(*plotProbs[probsEnt,"Probabilities for an entangled state"]*)


(* ::Input:: *)
(*plotProbs[probsSep,"Probabilities for a separable state"]*)


(* ::Subchapter::Closed:: *)
(*Measurement probabilities*)


(* ::Input::Initialization:: *)
ClearAll[probsEnt]
(*probability of detection of am entangled photon pair in one of four outputs, when detectors are rotated through angle \[Theta]*)
probsEnt[\[Theta]_]:={0.5,0,0,0.5}


(* ::Input::Initialization:: *)
ClearAll[probsSep]
(*probability of detection of a separable photon pair in one of four outputs, when detectors are rotated through angle \[Theta]*)
probsSep[\[Theta]_]:={1/8. (3+Cos[4 \[Theta]]),Cos[\[Theta]]^2 Sin[\[Theta]]^2,Cos[\[Theta]]^2 Sin[\[Theta]]^2,1/8. (3+Cos[4 \[Theta]])}


(* ::Subchapter::Closed:: *)
(*Photon sequences*)


(* ::Input::Initialization:: *)
ClearAll[generateSinglePhotonSequence]
generateSinglePhotonSequence[probs_,numOfPairs_,sampleGenFun_:sampleGenerationCustom]:=Module[{samples,histlist,seqPh},

(*generate a train of photons according to probabilities probs; the detection events are generated by function 'sampleGenFun'*)
(*because in the video only a moderate number of photons is used, the collected statistics given by sampleGenFun=sampleGenerationMathem differ quite significantly from the expected large-number averages; to counter this artefact, sampleGenFun=sampleGenerationCustom is chosen such that the resulting statistics follow more closely the expected averages at the cost of being not random *)
samples=sampleGenFun[probs,numOfPairs];
histlist=FoldList[Plus,samples];
seqPh=Rest[samples]/.{{0,0,0,1}->{True,False},{0,0,1,0}->{True,True},{0,1,0,0}->{False,False},{1,0,0,0}->{False,True}};
{AppendTo[seqPh,{False,False}],histlist}
]


(* ::Input::Initialization:: *)
ClearAll[sampleGenerationMathem]
(*random generation given by function RandomChoice*)
sampleGenerationMathem[probs_,numOfPairs_]:=Prepend[RandomChoice[probs->{{0,0,0,1},{0,0,1,0},{0,1,0,0},{1,0,0,0}},numOfPairs],{0,0,0,0}];


(* ::Input::Initialization:: *)
ClearAll[sampleGenerationCustom]
(*"random" generation that produces well-behaved statistics*)
(*detection events are built consecutively by looking at previous events and excluding those that differ too much from the expected values, see customRandomChoiceSingleRun*)
sampleGenerationCustom[probs_,numOfPairs_]:=NestList[customRandomChoiceSingleRun[probs,numOfPairs,#]&,{0,0,0,0},numOfPairs];


(* ::Input::Initialization:: *)
customRandomChoiceSingleRun[probs_,numOfPairs_,accum_,dev_:.8]:=Module[{samples,dists,entrs,argmin,randomness,batchSize=5},

(*accum are accumulated detections from previous events; this function generates a new event that closely follows the expected large-number averages*)
(*at first a batch of batchSize events is generated and only one event is chosen in the end according to two criteria*)
samples=RandomChoice[probs->{{0,0,0,1},{0,0,1,0},{0,1,0,0},{1,0,0,0}},batchSize];

(*to introduce "outliers", sometimes we use the standard approach*)
randomness=RandomChoice[{1-dev,1+dev}->{True,False}];
If[randomness,Return[samples[[1]]]];

(*otherwise we use the batch and find the event that is close to what we expect (1st criterion) and is also uniform enough (2nd criterion)*)
(*1st criterion calculates the distance between what we want and what we got*)
dists=Norm[numOfPairs probs-(accum+#)]&/@samples;
(*2nd criterion measures uniformity by calculating corresponding entropy*)
entrs=sampleEntropy[accum+#]&/@samples;
(*we want the distance to be small and entropy to be large*)
argmin=First@Ordering[dists-entrs,1];

(*return the "best" event*)
samples[[argmin]]
]


(* ::Input::Initialization:: *)
sampleEntropy[sample_]:=Module[{aux=sample},

(*Mathematica's built-in Entropy does not help here*)
aux=N[aux/.{0->Nothing}];
If[aux!={0,0,0,0},aux/=Total[aux]];
-aux.Log2[aux]
]


(* ::Chapter:: *)
(*Scene*)


(* ::Subchapter:: *)
(*Constants*)


(* ::Input::Initialization:: *)
fontFamily="Adobe Devanagari"(*"Arial"*)(*"Times New Roman"*);
fontSize=20;
grayColor=GrayLevel[0.41];
reCol=RGBColor[1,0.77,0](*Red*)
grCol=Magenta(*Green*)


(* ::Input::Initialization:: *)
With[{lab0="H",lab1="V"},
labelEnt=Text[Style[ToString[Ket["\[Psi]"],TraditionalForm]<>" = "<>ToString[HoldForm[1/Sqrt[2]],TraditionalForm]<>ToString[HoldForm[""(Ket[lab0,lab1]-Ket[lab1,lab0])],TraditionalForm],fontSize,FontFamily->fontFamily],Scaled@{.5,.88},{0,0}];
labelSep=Text[Style["\[Rho]"<>" = "<>ToString[HoldForm[1/2],TraditionalForm]<>ToString[HoldForm[""(Ket[lab0,lab1]Bra[lab0,lab1]+Ket[lab1,lab0]Bra[lab1,lab0])],TraditionalForm],fontSize,FontFamily->fontFamily],Scaled@{.5,.88},{0,0}];
]
(*{labelEnt,labelSep}*)


(* ::Subchapter:: *)
(*Source*)


(* ::Input::Initialization:: *)
(*credit to "J.M.'s discontentment"; https://mathematica.stackexchange.com/questions/49313/drawing-a-cuboid-with-rounded-corners*)
ClearAll[roundedCuboid]
roundedCuboid[p1_?VectorQ, p2_?VectorQ, r_?NumericQ]:=Module[{csk, csw, cv, ei, fi, ocp, osk, owt},
cv=Tuples[Transpose[{p1 + r, p2 - r}]];
ocp={{{1, 0, 0}, {1, 1, 0}, {0, 1, 0}}, {{1, 0, 1}, {1, 1, 1}, {0,1, 1}}, {{0, 0, 1}, {0, 0, 1}, {0, 0, 1}}};
osk={{0, 0, 0, 1, 1, 1}, {0, 0, 0, 1, 1, 1}};
owt={{1, 1/Sqrt[2], 1}, {1/Sqrt[2], 1/2, 1/Sqrt[2]}, {1,1/Sqrt[2], 1}};
ei={{{4, 8}, {2, 6}, {1, 5}, {3, 7}}, {{6, 8}, {2, 4}, {1, 3}, {5,7}}, {{7, 8}, {3, 4}, {1, 2}, {5, 6}}};
csk={{0, 0, 1, 1}, {0, 0, 0, 1, 1, 1}};
csw={{1, 1/Sqrt[2], 1}, {1, 1/Sqrt[2], 1}};
fi={{8, 6, 5, 7}, {8, 7, 3, 4}, {8, 4, 2, 6}, {4, 3, 1, 2}, {2, 1,5, 6}, {1, 3, 7, 5}};

Flatten[{EdgeForm[],BSplineSurface3DBoxOptions->{Method->{"SplinePoints" -> 35}},
MapIndexed[
               BSplineSurface[Map[AffineTransform[{RotationMatrix[\[Pi] Mod[#2[[1]] - 1, 4]/2, {0, 0, 1}], #1}],ocp.DiagonalMatrix[r {1,1,If[Mod[#2[[1]] - 1, 8] < 4, 1, -1]}],{2}
],SplineDegree->2,SplineKnots->osk,SplineWeights->owt]&
,cv[[{8, 4, 2, 6, 7, 3, 1, 5}]]
]
,
MapIndexed[
Function[{idx, pos},BSplineSurface[Outer[Plus, cv[[idx]],Composition[Insert[#,0,pos[[1]]]&,RotationTransform[\[Pi] (pos[[2]] - 1)/2]]/@(r {{1,0}, {1, 1}, {0, 1}}), 1]
,SplineDegree->{1, 2},SplineKnots-> csk,SplineWeights->csw]]
,ei,{2}
]
,
Polygon[MapThread[
Map[TranslationTransform[r #2],cv[[#1]]]&,{fi,Join[#,-#]&[IdentityMatrix[3]]}
]]}
]
]


(* ::Input::Initialization:: *)
ClearAll[sourceCuboid]
sourceCuboid[fine:(True|False):True,scale_:0.8]:=
sourceCuboid[fine,scale]=
 Module[{pt,cyl,outlet}, 
   pt = scale {1, 1, 1};
   cyl = {Black, Cylinder[{{-.2, 0, 0}, {0.1, 0, 0}}, 0.2]};
   
   {GrayLevel[.8], EdgeForm[None],
If[fine,
{
roundedCuboid[-pt, pt, .1],

outlet = First@Show@Region[
RegionProduct[BoundaryDiscretizeRegion@Annulus[{0, 0},{0.5,1}],Line[{{-.5}, {0.5}}]]
]/.x_Directive -> Directive[EdgeForm[None]];
outlet=Delete[outlet,{2,2,-1}];
Translate[#, scale {0, -1, 0}] &@Rotate[Scale[outlet, 0.3], \[Pi]/2, {1, 0, 0}],
Translate[#, scale {1, 0, 0}] &@Rotate[Scale[outlet, 0.3], \[Pi]/2, {0, 1, 0}]
}
,
Cuboid[-pt, pt]
],
Translate[cyl, scale {.9, 0, 0}],
Translate[Rotate[cyl, -(\[Pi]/2), {0, 0, 1}], scale {0, -0.9, 0}],
}];


(* ::Input:: *)
(*(*{Graphics3D[{sourceCuboid[True]}, Boxed -> False, Lighting -> "Neutral"],Graphics3D[{sourceCuboid[False]}, Boxed -> False, Lighting -> "Neutral"]}*)*)


(* ::Subchapter:: *)
(*Photon*)


(* ::Input::Initialization:: *)
photon={Orange,Ball[{0, 0, 0}, .1]};


(* ::Input::Initialization:: *)
travelFunction[gr_, pt1_, pt2_, rat_]:=Translate[gr, (1 - rat) pt1 + rat pt2]


(* ::Input::Initialization:: *)
ClearAll[photonTravelAll]
photonTravelAll[refl:(True|False),ptCr_, ptBS_, lenOut_, rat_,ang_]:=Module[{ratLoc, ptOut, distCrBS = Norm[ptCr - ptBS], distBSOut, distRatio,distTotal, incr},

If[refl,
incr=RotationTransform[ang, ptCr - ptBS][lenOut Cross[Normalize[ptCr - ptBS], {0, 0, 1}]];
,
incr=-lenOut Normalize[ptCr - ptBS];
];
ptOut=ptBS+incr;

distBSOut = Norm[ptBS - ptOut];
distTotal = distCrBS + distBSOut;
distRatio = distCrBS/distTotal;

If[rat <= distRatio,
ratLoc = rat distTotal/distCrBS;
travelFunction[photon, ptCr, ptBS, ratLoc]
,
ratLoc = (rat distTotal - distCrBS)/distBSOut;
travelFunction[photon, ptBS, ptOut, ratLoc]
]
]


(* ::Subchapter:: *)
(*PBS setup*)


(* ::Input::Initialization:: *)
pbs=Module[{p1={0, 0, 0},p2={1, 0, 0},p3={0, 1, 0},p4={0, 0, 1},p5={1, 0, 1},p6={0, 1, 1},prism},
prism=Translate[Prism[{p1, p2, p3, p4, p5, p6}], {-.505, -.505, -.5}];
{EdgeForm[None],
{Opacity[.8, Lighter[Blue, .7]],FaceForm[Opacity[.95, Lighter[Blend[{Cyan, Blue}, .2], .5]]],
prism
},
{Opacity[.9, Lighter[Blue, .7]],FaceForm[Opacity[.9, Lighter[Blend[{Cyan, Blue}, .4], .5]]],
Rotate[prism, \[Pi], {0, 0, 1}]
}
}
];


(* ::Input::Initialization:: *)
ClearAll[arrowStrap3D]
arrowStrap3D[pltstyle_:{},arrowlen_:0.7,arrwid_:0.2,strokegap_:.1,strokewid_:.1]:=Module[{maxang=2\[Pi]-strokegap,strapStroke,strapArrow,opts},

opts={Mesh->None,PlotStyle->pltstyle,Lighting->"Neutral"};

strapStroke=ParametricPlot3D[{Cos[ang],Sin[ang],u},{ang,0,maxang-arrowlen},{u,-strokewid/2.,strokewid/2.},PlotPoints->8,Evaluate[Sequence@@opts]];strapStroke=First@Cases[strapStroke,_GraphicsComplex,Infinity];

strapArrow=ParametricPlot3D[{Cos[ang],Sin[ang],u arrwid(maxang-ang)},{ang,maxang-arrowlen,maxang},{u,-1,1},PlotPoints->5,Evaluate[Sequence@@opts]];strapArrow=First@Cases[strapArrow,_GraphicsComplex,Infinity];

{strapStroke,strapArrow}
]


(* ::Input::Initialization:: *)
strap=arrowStrap3D[{Black},.7];
circle=ParametricPlot3D[{Cos[ang],Sin[ang],0},{ang,0,2\[Pi]}];
circle=First@Cases[InputForm[circle],_Line,Infinity];


(* ::Input::Initialization:: *)
ClearAll[pbsSetup]
pbsSetup[fine:(True|False):True,lengthOut_,color1_:grayColor,color2_:grayColor,arrow_:True]:=
pbsSetup[fine,lengthOut,color1,color2,arrow]=Module[{pt, plate, platescale = 0.5},
pt = platescale {0.05, 1, 1};
plate = Cuboid[-pt, pt];
{
pbs,
{Opacity[.8], EdgeForm[None],
{color1,Translate[plate, lengthOut {-1, 0, 0}]},
{color2,Translate[Rotate[plate, \[Pi]/2, {0, 0, 1}], lengthOut {0, 1, 0}]}
},
If[arrow,
{Dashed,Black,
If[fine,
Rotate[#,-\[Pi]/2,{0,1,0}]&@Rotate[Scale[strap,1.1lengthOut],\[Pi]/2,{1,0,0}],
Rotate[Scale[Arrow@circle,1.1lengthOut],\[Pi]/2,{1,0,0}]
]
},
{}
]
}
]


(* ::Input:: *)
(*(*{Graphics3D[pbsSetup[True,2,grayColor,grayColor,True] ],Graphics3D[pbsSetup[False,2,grayColor,grayColor,True] ]}*)*)


(* ::Subchapter:: *)
(*Histograms*)


(* ::Input::Initialization:: *)
ClearAll[histogramPhoton]
histogramPhoton[list_,ranmax_:5,problist_:{0,0,0,0},showTicks_:True]/;Length[list]==Length[problist]==4:=
Module[{ticks,minlim=-.2,maxlim=1.2,data,probs,colorTab,col1=grCol(*Green*),col2=reCol(*Red*)},

(*generate ticks*)
colorTab={{col1,col1},{col1,col2},{col2,col1},{col2,col2}};
ticks={#,Row[colorTab[[#]]],0}&/@Range[4];

(*plot of data*)
data=ListStepPlot[Transpose[{Range[0,5],Join[{0},list,{0}]}],Center,
LabelingFunction->If[showTicks,(Placed[Style[Round[#1[[2]]],FontSize->fontSize,FontFamily->fontFamily],Above]&),None],
PlotStyle->Directive[Orange,EdgeForm[None]],
Filling->Axis,FillingStyle->Directive[Orange,Opacity[.6]],
PlotRange->{{0.5,4.5},{minlim ranmax,maxlim ranmax}}
(*it is necessary to have PlotRange option present in both plots, otherwise the plots "wiggle" from one frame to the other, for some reason*)
];

(*plot of rescaled probabilities*)
probs=ListStepPlot[Transpose[{Range[0,5],Join[{0},problist,{0}]}],Center,
PlotStyle->Directive[Lighter[Blue,.5],EdgeForm[None]],
Filling->Axis,FillingStyle->Directive[Lighter[Blue,.7],Opacity[.7]],
PlotRange->{{0.5,4.5},{minlim ranmax,maxlim ranmax}}
(*it is necessary to have PlotRange option present in both plots, otherwise the plots "wiggle" from one frame to the other, for some reason*)
];

(*both plots together*)
Show[probs,data,Ticks->{ticks,None},ImageSize->250,Axes->{True,False}]
]


(* ::Input::Initialization:: *)
ClearAll[histogramPhotonBlank]
histogramPhotonBlank=histogramPhoton[{0,0,0,0},1,{0,0,0,0},False];


(* ::Input::Initialization:: *)
ClearAll[histogramPhotonProb]
histogramPhotonProb[problist_,ranmax_:1]/;Length[problist]==4:=Module[{ticks,minlim=-.2,maxlim=1.2,probs,colorTab,probsAux,col1=grCol(*Green*),col2=reCol(*Red*)},

(*generate ticks*)
probsAux=If[Total[problist]==0,{0,0,0,0},Round[N[problist/Total[problist]],.01]];
probsAux=Join[{0},probsAux,{0}];
(*colorTab={{Green,Green},{Green,Red},{Red,Green},{Red,Red}};*)
colorTab={{col1,col1},{col1,col2},{col2,col1},{col2,col2}};
ticks={#,Row[colorTab[[#]]],0}&/@Range[4];

(*plot of rescaled probabilities*)
probs=ListStepPlot[Transpose[{Range[0,5],Join[{0},problist,{0}]}],Center,
LabelingFunction->(Placed[Style[probsAux[[#2[[2]]]],FontSize->fontSize,FontFamily->fontFamily],Above]&),
PlotStyle->Directive[Lighter[Blue,.5],EdgeForm[None]],
Filling->Axis,FillingStyle->Directive[Lighter[Blue,.7],Opacity[.7]],
PlotRange->{{0.5,4.5},{minlim ranmax,maxlim ranmax}},
Ticks->{ticks,None},ImageSize->250,Axes->{True,False}
];

probs
]


(* ::Input:: *)
(*(*histogramPhotonProb[{0,0,0,0}(*{1,3,5,5}*),5]*)*)
(*(*histogramPhoton[{1,3,5,5},5,{3,2,2.5,3},True]*)*)


(* ::Subchapter:: *)
(*Scene*)


(* ::Input::Initialization:: *)
ClearAll[scene]
scene[fine_,ang_, ratIn_, refla_: True,reflb_: True,angViewIn_:0.1,arrows_:True,angleLabel_:True,label_:Text["",{0,0}],imgSize_:Automatic,sphRad_:Automatic] :=
 Module[{scale=3,ptCr,ptBS1,ptBS2,reflColor=reCol,transColor=grCol,lenghtOut=2,color1a=grayColor,color2a=grayColor,color1b=grayColor,color2b=grayColor,rat,angView,angleLab,aliceLab,bobLab,detFireLim=0.8},

{ptCr,ptBS1,ptBS2} = scale{ {0, 1, 0},{-1, 0, 0},{1, 0, 0}};
angleLab=If[angleLabel,Text[Style["\[Theta] = "<>ToString[Round[Mod[ang ,2.\[Pi]]/Degree,.1]]<>"\[Degree]",fontSize,FontFamily->fontFamily],Scaled@{.95,.73},{-1,0}],{}];
aliceLab=Text[Style["A",Bold,1.5fontSize,FontFamily->fontFamily],{0.1,0.03}];
bobLab=Text[Style["B",Bold,1.5fontSize,FontFamily->fontFamily],{0.9,0.03}];

rat=Clip[ratIn,{0.,1}];
angView=Clip[angViewIn,{0.1,\[Pi]}];

If[rat>detFireLim,
If[refla,color1a=reflColor,color2a=transColor];
If[reflb,color1b=transColor,color2b=reflColor]
];

Graphics3D[{
Translate[Rotate[Rotate[pbsSetup[fine,lenghtOut + 0.1,color1a,color2a,arrows], 3 \[Pi]/4., {0, 0, 1}], ang,ptCr - ptBS1], ptBS1],
Translate[Rotate[Rotate[pbsSetup[fine,lenghtOut + 0.1,color1b,color2b,arrows], - 3\[Pi]/4., {0, 0, 1}], ang+\[Pi],ptCr - ptBS2], ptBS2],
Translate[Rotate[sourceCuboid[fine], -\[Pi]/4., {0, 0, 1}], ptCr],
photonTravelAll[refla,ptCr,ptBS1,lenghtOut,rat,ang],
photonTravelAll[reflb,ptCr,ptBS2,lenghtOut,rat,ang+\[Pi]]
},
Boxed -> False, Lighting -> "Neutral",ViewCenter -> {0.6, 0.9, 0.5},ViewPoint -> FromSphericalCoordinates[{1,angView,-\[Pi]/2}], ViewVertical -> {0, 1, 0},
ImageSize -> If[imgSize===Automatic,250{1.6,1},imgSize],
SphericalRegion->Sphere[{0,1.8,0},If[sphRad===Automatic,4.258,sphRad]],
Epilog->{label,angleLab,aliceLab,bobLab}
]
]


(* ::Input::Initialization:: *)
ClearAll[sceneHist]
sceneHist[sceneEnt_,sceneSep_,histEnt_,histSep_,histSize_:220]:=
Grid[{
{sceneEnt,sceneSep},
{Show[histEnt,ImageSize->histSize],Show[histSep,ImageSize->histSize]}
},Alignment->Center
]


(* ::Input:: *)
(*(*Manipulate[*)
(* scene[fine,ang, rat, refla,reflb,av,arrows,True,labelEnt,250{1.6(*1.328125`*),1},5(*4.258*)], {ang, 0, 2 \[Pi]}, {rat, 0, 1.1}, {refla, {True, False}}, {reflb, {True, False}},{{av,0.1(*0.93*)},0.1,0.938},{{arrows,True},{True,False}},{{fine,True},{True,False}},Deployed\[Rule]True]*)*)


(* ::Chapter:: *)
(*Video*)


(* ::Input::Initialization:: *)
modDivRatio[rat_,num_]:={num Mod[rat,1/num],Floor[rat (num)]+1}


(* ::Subchapter:: *)
(*Smooth rotation of detectors*)


(* ::Input::Initialization:: *)
ClearAll[probabsVideoSegment]
probabsVideoSegment[fine_,max_,numOfPairs_,angInit_,angFinal_,angView_:0.1,arrows_:False,angleLabel_:True,ratRotLimit_:0.9]:=Module[{playSegment,paddAng=0.2(angFinal-angInit)},

(*generate function that governs the stage where detectors smoothly rotate*)
playSegment[ratIn_]:=Module[{rat,gr,histEnt,histSep,sceneEnt,sceneSep,ang,cond},

rat=Clip[ratIn,{0,1}];
cond=rat<=ratRotLimit;
ang=Rescale[rat,{0,ratRotLimit},{angInit-paddAng,angFinal+paddAng}];
ang=Clip[ang,{angInit,angFinal}];

sceneEnt=scene[fine,ang,0,False,False,angView,If[cond,arrows,False],If[cond,angleLabel,False],labelEnt];
sceneSep=scene[fine,ang,0,False,False,angView,If[cond,arrows,False],If[cond,angleLabel,False],labelSep];
histEnt=histogramPhotonProb[numOfPairs probsEnt[ang],max];
histSep=histogramPhotonProb[numOfPairs probsSep[ang],max];

gr=sceneHist[sceneEnt,sceneSep,histEnt,histSep];
gr
];

playSegment
]


(* ::Subchapter:: *)
(*Photon emission and detection*)


(* ::Input::Initialization:: *)
ClearAll[photonsVideoSegment]
photonsVideoSegment[fine_,max_,numOfPairs_,seqEnt_,seqSep_,angBounds_,angView_:0.1,arrows_:False,angleLabel_:True,ratPhotonsStart_:0.2,ratPhotonsEnd_:0.9]:=Module[{playSegment,seqPhEnt,seqPhSep,histListEnt,histListSep,histogramPhotonListEnt,histogramPhotonListSep,paddAng,angInit=angBounds[[1]],angFinal=angBounds[[2]]},

(*sequences of photons are generated in a separate function and piped to this function*)
{seqPhEnt,histListEnt}=seqEnt;
{seqPhSep,histListSep}=seqSep;
paddAng=0.3(angFinal-angInit);

histogramPhotonListEnt=histogramPhoton[#,max,numOfPairs probsEnt[angFinal]]&/@histListEnt;
histogramPhotonListSep=histogramPhoton[#,max,numOfPairs probsSep[angFinal]]&/@histListSep;

(*generate function that governs the stage where photons are emitted by the source and then detected by rotated detectors*)
playSegment[ratIn_]:=Module[{rat,idx,ratLoc,gr,histEnt,histSep,sceneEnt,sceneSep,ang},

rat=Clip[ratIn,{0,1}];
ang=Rescale[rat,{0,ratPhotonsStart},{angInit-paddAng,angFinal+paddAng}];
ang=Clip[ang,{angInit,angFinal}];

Which[
rat<=ratPhotonsStart,
(*at first, detectors are rotated*)
sceneEnt=scene[fine,ang,0,False,False,angView,arrows,angleLabel,labelEnt];
sceneSep=scene[fine,ang,0,False,False,angView,arrows,angleLabel,labelSep];
histEnt=histogramPhotonBlank;
histSep=histogramPhotonBlank;
,
rat<=ratPhotonsEnd,
(*second, a train of photons is emitted and detected*)
ratLoc=Rescale[rat,{ratPhotonsStart,ratPhotonsEnd},{0,1}];
{ratLoc,idx}=modDivRatio[ratLoc,numOfPairs];
sceneEnt=scene[fine,ang,ratLoc, seqPhEnt[[idx,1]],seqPhEnt[[idx,2]],angView,False,angleLabel,labelEnt];
sceneSep=scene[fine,ang,ratLoc, seqPhSep[[idx,1]],seqPhSep[[idx,2]],angView,False,angleLabel,labelSep];
histEnt=histogramPhotonListEnt[[idx]];
histSep=histogramPhotonListSep[[idx]];
,
True,
(*at last, some time is left after the last photon gets detected*)
idx=Length[histogramPhotonListEnt];
sceneEnt=scene[fine,ang,0, seqPhEnt[[idx,1]],seqPhEnt[[idx,2]],angView,False,angleLabel,labelEnt];
sceneSep=scene[fine,ang,0, seqPhSep[[idx,1]],seqPhSep[[idx,2]],angView,False,angleLabel,labelSep];
histEnt=histogramPhotonListEnt[[idx]];
histSep=histogramPhotonListSep[[idx]];
];

gr=sceneHist[sceneEnt,sceneSep,histEnt,histSep];
gr
];

playSegment
]


(* ::Subchapter:: *)
(*Short video*)


(* ::Input::Initialization:: *)
ClearAll[generateVideoShort]
generateVideoShort[numOfPairs_,angDelta_,initAngle_,angNum_,angView_,fine:(True|False):True]:=Module[{max,stages,videoFun,angList,angFin,funList,photonSeriesInit,rotationViewInit,photonSeriesMid,probabRotatSetup,rotationViewFinal,photonSeriesList,arg,seqsEnt,seqsSep,argList,compareList,rescaleList},

(*generate the list of angles for which measurements are performed and corresponding photon statistics*)
angList=initAngle+angDelta Range[0,angNum];
angFin=Last[angList];
seqsEnt=generateSinglePhotonSequence[probsEnt[#],numOfPairs]&/@angList;
seqsSep=generateSinglePhotonSequence[probsSep[#],numOfPairs]&/@angList;
max=Max[seqsEnt/.True|False->0,seqsSep/.True|False->0];

(*the middle stage with photon detection in different bases*)
angList=Partition[angList,2,1];
argList=Transpose[{Rest@seqsEnt,Rest@seqsSep,angList}];
photonSeriesList=photonsVideoSegment[fine,max,numOfPairs,#1,#2,#3,angView,True,True]&@@@argList;
photonSeriesMid[x_]:=Module[{ratLoc,idx},
{ratLoc,idx}=modDivRatio[x,angNum];
photonSeriesList[[idx]][ratLoc]
];

(*all the other stages of the video*)
photonSeriesInit=photonsVideoSegment[fine,max,numOfPairs,First[seqsEnt],First[seqsSep],{initAngle,initAngle},angView,False,False,0.01];
probabRotatSetup=probabsVideoSegment[fine,max,numOfPairs,angFin,angFin+(*2.\[Pi]+*)Mod[initAngle-angFin,2\[Pi]],angView,True,True];

(*list of all functions and time delimiters*)
funList={photonSeriesInit,photonSeriesMid,probabRotatSetup};
stages={0.2,0.8}; (*time instants when one stage should change into another*)

(*from here on a general code...*)
stages=Prepend[stages,0];
compareList=Table[stages[[j]]<=arg<stages[[j+1]],{j,Length[funList]-1}];
rescaleList=Table[Rescale[arg,{stages[[j]],stages[[j+1]]}],{j,Length[funList]-1}];

(*generate function that governs the flow of the video*)
videoFun[ratIn_]:=Module[{rat,fun},

rat=Clip[ratIn,{0,1}];

(*choose a correct function from the list*)
fun=Piecewise[Transpose[{Most@funList,compareList}]/.arg->rat,Last[funList]];

(*choose a correct rescaling for the input parameter*)
rat=Piecewise[Transpose[{rescaleList,compareList}]/.arg->rat,Rescale[rat,{Last[stages],1}]];

(*return value*)
fun[rat]
];

videoFun
]


(* ::Chapter:: *)
(*Rasterization*)


(* ::Input::Initialization:: *)
ClearAll[rasterizeFrameSequence]
rasterizeFrameSequence[fun_,numOfFrames_:10,imgResolution_:70]:=Module[{time,frames},
{time,frames} =AbsoluteTiming[
ParallelMap[
Rasterize[fun[#],Background -> None,ImageResolution ->imgResolution]&,
Subdivide[0, 1.,numOfFrames-1]
]
];
Print["execution time: ",DateString[time,{"Minute"," m ","Second"," s"}]];
Print["size: ",ByteCount[frames]/1024/1024.," MB"];

frames
]


(* ::Chapter:: *)
(*Export*)


(* ::Input:: *)
filename="movieShort.gif";


(* ::Input:: *)
funAllshort=generateVideoShort[8,3\[Pi]/8.,0.,2,0.9,True];


(* ::Input:: *)
framesShort=rasterizeFrameSequence[funAllshort,210,60];


(* ::Input:: *)
(*(*ListAnimate[framesShort,AnimationRate->3.]*)*)


(* ::Input:: *)
SetDirectory[NotebookDirectory[]]
Export[filename, framesShort,"DisplayDurations"->0.2,AnimationRepetitions->Infinity,"ColorMapLength"->256,Dithering->None]
FileSize[filename]