      program simf_b5  EPR CORRELATION
c=======================================================================
c
c    Simulates the EPR correlations from coincidence experiments
c 
c    Shows violation of Bell's Inequalities using Quaternion Spin
c
c  The model gives a superpostion of two phases for spin in free flight. 
c  For polarizations only one phase is present while for coherence there 
c  are two, and we must choose the larger. Then the sign of the presession
c  (spin)is determined by varying the parameters and the LHV in the click function.
c
c   thA     Alice's filter setting deg
c   thARad  in radians
c   thB     Bob's filter setting deg
c   thBRad  in radians
c   thAB    the difference deg  degAB = degA - degB 
c   thABRad in radians
c 
c   the integer LHV:
c   the states are  |k><kp|, +/- each:  
c      k .EQ. kp  is polarization; 
c      k .NE. kp  is coherence.
c
c   The program chooses the difference thAB from 0 to 360 in intervals of
c   5 degree and calculates the correlation EAB = EABP + EABC as the sum 
c   of the correlation from polarizations EABP and that for coherence EABC.
c   For each thAB the program chooses cobinations of thA and thB to include 
c   different Alice and Bob settings. 
c             
c   The Phase depends upon the signs of the trig functions of the angles from
c   the theory. A Left spinning spin is +1 and -1 for a Rght spinning spin. 
c   From these parameters a click is produced in the function click. 
c 
c       superposition is derived in Appendix 1
c 
c         A_(axis(1,theta))*exp(i(pi/2 - thA)Y) + A_(axis(3,theta))*exp(-i(thA)Y) 
c 
c   For polarizations: choose either axis and find if the imaginary part spins L or R
c   for click +1 or -1,
c 
c   For coherences, first choose the larger of axis 1 or 3 from the real part, then find the 
c   clicks in the same way as for polarization. 
c
c   th (thRad) is the LHV theta which is the orientation of the 2D spin in the 
c   LAB frame.  
c 
C   Bob's ClickB is found by passing again but with the following change:  
c 
c       th --> th + Pi  (or Alice has th+Pi/2 and Bob th-Pi/2, th must differ by 
c       Pi from Alice to Bob for conservation of spin ang. mom.)
c
c   The correlation is calcuated separately for polarization EABP, and coherence EABC, 
c   then they are added since they are independent and complementary attributes of spin. 
c   The clicks from Alice and Bob are multiplied to give clickA*clickB which is even or odd. 
c   If this is +1 then they accumulate in EQ, and if -1 they accumulat in NEQ, (for both P and C)
c   At the end of the calculation the correlation is calculated from
c      
c        EABP(thAB) = (EQP - NEQP)/(EQP + NEQP)
c        EABC(thAB) = (EQC - NEQC)/(EQC + NEQC)
c        EAB = EABP + EABC  is the total correlation which is stored in correl(i)  
c
c   To compile:  g77 sim_b5.f
c   To plot, use say gnu plot 
c     
c   The simulated data is written to "points.txt" Unit 10;
c   The raw data of clicks is written to "check.txt" Unit 11. 
c   The CHSH value is written to unit 12.
c
c      Bryan Sanctuary
c      bryan.sanctuary@mcgill.ca
c      December 2022
c
c=======================================================================
c                                                                          
        implicit none
		integer          k,kk,kp,kpp,i,count1           
        integer          thAB,thABn,thA,thAn,thAan,thBan,thB,thBn,thn           
		integer          click,clickA,clickB 
		integer          angles(361)
	    double precision correl(361),correEQP(361),correNEQP(361)
		double precision correEQC(361),correNEQC(361)
	    double precision Phase,Pi,EQ,NEQ,EAB,EQC,NEQC,EABC,EQP,NEQP,EABP	
		double precision thABRad,thARad,thBRad,thRad,th,CHSH            

	   
	    open (10, file ='points.txt', status ='replace')
        open (11, file ='check.txt',  status ='replace')
        open (12, file ='CHSH.txt',   status ='replace')	   
c
        Pi = 4.0*atan(1.0)
c 
c   Initialize the storage of the results
c 
        do 1031 i    = 1,361
        correl(i)    = 0
	    angles(i)    = 0
        correEQP(i)  = 1  
		correNEQP(i) = 1
        correEQC(i)  = 1
		correNEQC(i) = 1		
1031    continue
c
c   set the counter 
c
	    count1 = 1
c 
c   set the angle difference (the plot parameter) of thAB = (thA - thB)
c
	    do 1 thABn   = 5,365,5
	         thAB    = thABn - 5
	         thABRad = thAB*Pi/180.
c	
c   initialize these to  zero: integer counts are large and set to 1 to 
c   avoid divide errors. 
c
	    EAB  = 0
	    EABP = 0
	    EABC = 0
	    EQ   = 1
	    EQP  = 1
	    NEQP = 1
	    NEQC = 1
	    EQC  = 1
c 
c  obtain Alice's filter angle thArad and Bob's as thBrad
c 
        do 2   thBan  = 5,365,5
	         thB    = thBan - 5
	         thA    = thAB+thB  
	         thARad = thA*Pi/180.
	         thBRad = thB*Pi/180.
c 
c   set the LHV: the k and kp values to distinguish polarization (k = kp)  
c   from coherence (k = -kp)
c 
           do 3 kk  = 1,3,2
	            k   = kk-2
	       do 4 kpp = 1,3,2
                kp  = kpp-2
c               kp  = -k 
c               kp  = +k
c
c  this calculates both polarization and coherence and their sum. To calculate only one,
c  first comment out the do loop 4 above.  Set kp = k for plarizations
c  and kp = -k for coherenes.
c
c  set the LHV theta (thRad) which differs between Alice and Bob by pi. 
c  for faster running, change the increment to say 5.
c
 	    do 7 thn    = 1,361,1
	         th     = thn-1
             thRad  = th*Pi/180.
c 
c   determine the clicks from the function click
c
        clickA = click(k,kp,nX,nZ,thrad-pi/2.,thArad)
        clickB = click(k,kp,nX,nZ,thrad+pi/2.,thBrad) 
c 
c  Accumulate the clicks for polarization separately from the coherence
c
c 	    accumulate polarization clicks
c
        if (k.eq.kp) goto 90
	    goto 91
90	    if (clickA*clickB.eq.1) then
	    EQP = EQP + 1
	    elseif (clickA*clickB.eq.-1) then
	    NEQP = NEQP + 1
	    endif
        goto 93	 
c 
c       accumulate coherence clicks   
91	    if (k.ne.kp) goto 92
        goto 93
92	    if (clickA*clickB.EQ.1) then
	    EQC = EQC + 1
	    else if (clickA*clickB.EQ.-1) then
	    NEQC = NEQC + 1
	    endif
c
93      continue	 
7       continue 
4       continue
3       continue
2       continue
c
c  calculate the correlation for the polarizations and the coherences and then add the two
c 
	    EABP = (EQP - NEQP)/(EQP + NEQP)
	    EABC = (EQC - NEQC)/(EQC + NEQC)
	    EAB  = (EABP + EABC)	
c 
c  store the points and increment count1 for the next angle thAB.

		 correEQP (count1)    =  EQP
		 correNEQP(count1)    =  NEQP
         correEQC (count1)    =  EQC
		 correNEQC(count1)    =  NEQC

		 
	     angles   (count1)     = thAB
	     correl   (count1)     = EAB
c 
	    count1                = count1 + 1
c 
1       continue
c 	
c    calculate CHSH
c 
        CHSH = abs(3.*correl(10)-correl(28)) 
c
c       write results to a file to plot
c   
 	    do 1000 i = 1,73
	    write(10,997) angles(i), correl(i) 
997     format(I3,f12.8)
c
c    write the raw data clicks
c 
		write(11,700)angles(i),correEQP(i),correNEQP(i),correEQC(i),
     +  correNEQC(i)	
700     format(1x,I4,4f20.0)	
1000    continue 
	    write(12,1002)CHSH,correl(10),correl(28)
1002    format("CHSH = ",f8.3," at 45 ",f10.5," at 135 degrees ",f10.5)
	    stop
        end 
c=======================================================================   
	    integer function click(k,kp,thrad,thArad)
        integer          k,kp,nZ,nX       
        double precision thRad,thARad,Pi,PhaseM,PhaseMO,PhaseME,PhaseS2
	    double precision Phase,PhaseCC,PhaseCS,PhaseSC,PhaseSS,PhaseS1
c 
        pi = 4.0*atan(1.0)  
c
c    Toggle between polarizations and coherences
c
c       polarizations--only give inverted triangle: CHSH = 2.  Can use either sin or cos. 
c 
	    if (k.eq.kp) goto 130
 	    goto 131
130	    continue  	
         Phase =  cos(thArad )*(nZ*cos(thRad)-nX*sin(nX*thRad))
 	    if (Phase.GT.0.) then
 	    click = +1
 	    elseif (Phase.LE.0.) then
 	    click = -1
	    endif 
	    goto 132
131     continue
c 
c       The model below can be changed to test others. 
c      
        if (k.eq.kp) goto 132
c 
c   choose the larger spin axis	
c 	
	    PhaseCC = cos(2*thArad )*(nZ*cos(thRad)-nX*sin(thRad))
	    PhaseSS = sin(2*thArad )*(nZ*sin(thRad)+nX*cos(thRad))
	    PhaseM  = max(PhaseCC,PhaseSS)   
		if (phaseM.EQ.PhaseCC) then
	    Phase = sin(thArad)*(nZ*cos(thRad)-nX*sin(thRad))
	    elseif (phaseM.EQ.PhaseSS) then
	    Phase = cos(thArad)*(nZ*sin(thRad)+nX*cos(thRad))
	    endif
        if (Phase.GT.0.) then
        click   = +1
 	    elseif (Phase.LT.0.) then
  	    click   = -1
        endif
132     continue
123     return
 	    end
c=======================================================================