View source for Simulated optics experiment/Data analysis

{{draft task}}
{{alertbox|#E9EC6B|There was a typographical error in the formula for CHSH contrast, in which ''kappaL2R2'' appeared where ''kappaL2R1'' belonged. Both task description and Python code contained the error, which was discovered when Object Icon gave different results on the same raw data. Please review your examples to see if this error was copied. We apologize for the inconvenience.}}

===Introduction===

In this task, you will write a program to analyze the "raw data" outputted by one of the ''[[Simulated optics experiment/Simulator]]'' implementations. The analysis will estimate the following numbers:
* Four ''[[wp:Correlation_coefficient|correlation coefficients]]'': one for each of the four possible combinations of the two polarizing beam splitter angles. These measure the interrelatedness of light detection events on the two sides of the experiment.
* A ''[[wp:CHSH_inequality|CHSH contrast]]'', which is derived from the correlation coefficients.
The estimated CHSH contrast is what we are really interested in, because what we intend to show is that it contradicts what many physicists believe.

Quantum mechanics predicts a CHSH contrast of approximately <math>2\sqrt{2}\approx{2.828427}</math>. The ''[[Simulated optics experiment/Simulator]]'' simulation will be shown (by our analysis) to have a CHSH contrast very close to that. Note that our simulation does not employ any quantum mechanics, but instead simulates light [[wp:Classical_physics|"classically"]]. Also our simulation clearly involves only causation that in each instance takes a finite amount of forward-moving time. That is, there is no "[[wp:Action_at_a_distance|instantaneous action at a distance]]".

The prevailing view in physics, however, is that a CHSH contrast greater than 2 (which we certainly achieve) is possible only with quantum mechanics, and also that taking measurements on one side of the apparatus "instantaneously" determines what the measurements will be on the other side. Although quantum mechanics actually has no concept of "causation", this "determination" is often described as being "instantaneous action at a distance".

Thus, we are left with a conundrum: we clearly get a CHSH contrast not only greater than 2, but in fact very close to the theoretical value that quantum mechanics predicts. This is, in fact, a validation of the quantum mechanical mathematics: the mathematics is producing the correct result for this kind of experiment. On the other hand, the ''strong'' prevailing view in physics is that ''this cannot be happening'' if the physics involved is "classical"--which our simulated physics very much is. The prevailing view is that experiments similar to ours ''prove'' that "instantaneous action at a distance" is a real thing that exists. One sometimes sees such "proof" touted in the popular media.

I leave it to each individual to think about that conundrum and come to their own opinions. I am not a physicist and probably neither are you, so we can think whatever we want, and no one cares. But we likely ''are'' both computer programmers. What we mainly want to do here is demonstrate the use of computers to simulate physical events and analyze the data generated.

===Task description===

You will write a program to perform the mathematical data analysis described in the [[wp:Open_access|open access]] paper
    A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
    J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
    https://doi.org/10.2991/jnmp.2004.11.s1.13
The program will read the raw data that was outputted by any of the ''[[Simulated optics experiment/Simulator]]'' implementations. The [[#ATS|ATS]], [[#ObjectIcon|Object Icon]], [[#Python|Python]], and [[#RATFOR|RATFOR]] programs below can serve as reference implementations. The ATS program stores the raw data in a big array. The RATFOR does not store the raw data at all, but instead calculates as it inputs. The other programs use ''set'' data structures.

The analysis proceeds as follows.
# Read in the raw data. You do ''not'' have to retain the order of the data lines. In other words, you can, if you wish, store the data in a ''set'' structure.
# For each of the data lines, if its first entry is a "1", then swap the last two entries with each other, and do the same with the two entries before them. This is to account for the difference in geometry: the original light pulses are rotated 90° with respect to the apparatus, relative to if the first entry were a "0".
# Split the data into four groups, based on the values of the second and third line entries. In other words, one group for lines that look like <math>\circ~0~0~\circ~\circ~\circ~\circ</math>, one for lines that look like <math>\circ~0~1~\circ~\circ~\circ~\circ</math>, one for lines that look like <math>\circ~1~0~\circ~\circ~\circ~\circ</math>, and one for lines that look like <math>\circ~1~1~\circ~\circ~\circ~\circ</math>. Thus each group corresponds to a particular setting of the two polarizing beam splitters.
# Separately compute a correlation coefficient for each group. This is done as follows.
## Let <math>N</math> equal the total number of lines in a group.
## Let <math>N_{L1}</math> equal the number of lines that have a "1" in the fourth entry. This is the number of detections by one of the light detectors on the left side.
## Let <math>N_{L2}</math> equal the number of lines that have a "1" in the fifth entry. This is the number of detections by the other light detector on the left side.
## Add up similar numbers for the right side: <math>N_{R1}</math> for the sixth entry in a line, and <math>N_{R2}</math> for the seventh entry.
## Compute the following numbers: <math>sin_L = \sqrt{N_{L1}/N}</math>, <math>cos_L = \sqrt{N_{L2}/N}</math>, <math>sin_R = \sqrt{N_{R1}/N}</math>, and <math>cos_R = \sqrt{N_{R2}/N}</math>. These numbers represent an estimate of the orientations of the two beam splitters, by inference from detection events in the light detectors.
## Compute <math>cos_{LR} = (cos_R\cdot cos_L) + (sin_R\cdot sin_L)</math> and <math>sin_{LR} = (sin_R\cdot cos_L) - (cos_R\cdot sin_L)</math>. These are angle difference formulas from trigonometry. We are "eliminating the [[wp:Coordinate_system|coordinate system]]".
## Compute <math>\kappa = cos_{LR}^2 - sin_{LR}^2</math>. This is the correlation coefficient for the group, estimated by inference from detection events.
# Designate each of the four correlation coefficients according to which group it is for. The [[#Python|Python example]], for instance, refers to them as <code>kappaL1R1</code>, <code>kappaL1R2</code>, <code>kappaL2R1</code>, and <code>kappaL2R2</code>.
# Compute the estimated CHSH contrast: <code>-kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2</code>.
# Display the number of "pulse pair events" that were simulated, the four correlation coefficients, and the estimated CHSH contrast. Also display the theoretical value of the CHSH contrast (<math>2\sqrt{2}\approx{2.828427}</math>) and the difference between the estimated CHSH contrast and the theoretical value. Also display what we will call the ''CHSH violation'': the difference between the estimated CHSH contrast and the number 2. Any value significantly greater than zero seems to be ''supposed to'' mean that the simulation is inherently "quantum mechanical" and involves "instantaneous action at a distance". But you can decide for yourself whether it really means that. Just do the calculations and print them. One should get an estimated CHSH contrast that is close to the theoretical value, and thus a relatively "huge" CHSH violation of about +0.8.

===Extra "credit"===

Here is a snapshot of what contains a [[wp:Thought_experiment|thought experiment]] I have been working on: http://web.archive.org/web/20230529223105/https://crudfactory.com/Bell-assumes-his-conclusion.pdf The thought experiment concerns the same subject as this task, but more directly addresses the conundrum. It aims to illustrate the fallaciousness (despite their widespread acceptance by people regarded as "experts") of certain arguments underlying the conundrum, and thus to resolve the technical aspects of the conundrum. Whether it succeeds is up to you to decide. I am just a computer programmer, after all. And, even if I succeed in resolving the technical aspects of the conundrum, I have barely touched on the social aspects of it. Nevertheless, the thought experiment begs to be computer-animated, and I am not ''that'' kind of a programmer.

For extra "credit" (in the form of contribution to the community), write an animation or animations of the thought experiment.

--[[User:Chemoelectric|Chemoelectric]] ([[User talk:Chemoelectric|talk]]) 22:50, 29 May 2023 (UTC)

===See also===
* ''[[Simulated optics experiment/Simulator]]''
* E. T. Jaynes, [https://web.archive.org/web/20200915075028/https://bayes.wustl.edu/etj/articles/cmystery.pdf "Clearing up mysteries -- the original goal"] Opening lecture of the eighth International MAXENT Workshop. 1988. This is a classic from the same physicist as the [[wp:Jaynes–Cummings_model|Jaynes-Cummings model]] is named for.

=={{header|ATS}}==
This implementation stores the data in an n-by-7 runtime-allocated array, and does not use ''set'' data structures. It should prove translatable to C, Ada, modern Fortran, etc., without much trouble.

I know ATS can be very difficult to comprehend if you do not know the language, but there is not much "ATSism" in this code. It likely is readable by many. The following peculiarities do occur: some of the loops are written as tail recursions (which Schemers, ML-users, etc., probably will have no difficulty with), and there is C code embedded in the program in more than one way. I embed C code in part to avoid using libraries other than the ATS prelude, which is minimal.

The program will have to be linked with the C math library and with an allocator, which can be malloc(3) (despite free(3) never being called), because the only allocation is for the array:
<pre>patscc -O2 -std=gnu2x -DATS_MEMALLOC_LIBC optics_simulation_analysis_task.dats -lm</pre>
(It is my habit to tell patscc to use C beyond -std=c99)

<syntaxhighlight lang="ats">
// Reference:
//
//  A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
//  J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
//  https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)

#include "share/atspre_staload.hats"

typedef zero_or_one = intBtwe (0, 1)

val angleL1 = 0.0
and angleL2 = 45.0
and angleR1 = 22.5
and angleR2 = 67.5

(*------------------------------------------------------------------*)

%{^

#include <math.h>
#include <stdio.h>

//
// fscanf(3) will be a most convenient way to input the necessary
// data. It will ignore line breaks, but the line breaks are in the
// raw data for the convenience of line-by-line input, not because
// they matter.
//
static atstype_uint
_read_uint (FILE *f)
{
  unsigned int u = 0U - 1U;     // Yes, this is legal in standard C.
  int bitbucket = fscanf (f, "%u", &u);
  return u;
}

%}

fn
read_uint (inpf : FILEref) : uint =
  let
    typedef FILEstar = $extype"FILE *"
    extern castfn FILEref2star : FILEref -<> FILEstar
    extern fn _read_uint : FILEstar -> uint = "mac#_read_uint"
  in
    _read_uint (FILEref2star inpf)
  end

fn
read_intGte0 (inpf : FILEref) : intGte 0 =
  let
    val u = g1ofg0 (read_uint inpf)
    prval () = lemma_g1uint_param u (* Proves u >= 0. *)
  in
    u2i u
  end

fn
read_zero_or_one (inpf : FILEref) : zero_or_one =
  let
    val i = read_intGte0 inpf
  in
    if i <= 1 then
      i
    else
      begin
        fprintln! (stderr_ref, "bad data");
        exit (1)
      end
  end

(*------------------------------------------------------------------*)

fn
read_raw_data (inpf : FILEref) : mtrxszref zero_or_one =
  let
    val n = read_intGte0 inpf
    val data = mtrxszref_make_elt<zero_or_one> (i2sz n, i2sz 7, 0)
    var i : intGte 0
  in
    for (i := 0; i2sz i <> data.nrow(); i := succ i)
      let
        var j : intGte 0
      in
        for (j := 0; j <> 7; j := succ j)
          data[i, j] := read_zero_or_one inpf
      end;
    data
  end

fn {a : t@ype}
mtrxszref_swap (M  : mtrxszref a,
                i1 : int,
                j1 : int,
                i2 : int,
                j2 : int) : void =
  let
    val m1 = M[i1, j1] and m2 = M[i2, j2]
  in
    M[i1, j1] := m2;  M[i2, j2] := m1
  end

fn
adjust_for_geometry (data : mtrxszref zero_or_one) : void =
  let
    var i : intGte 0
  in
    for (i := 0; i2sz i <> data.nrow(); i := succ i)
      if data[i, 0] = 1 then
        begin
          mtrxszref_swap (data, i, 3, i, 4);
          mtrxszref_swap (data, i, 5, i, 6)
        end
  end

fn
group_size (data   : mtrxszref zero_or_one,
            group  : @(zero_or_one, zero_or_one)) : intGte 0 =
  let
    fun
    loop (i     : intGte 0,
          count : intGte 0) : intGte 0 =
      if i2sz i = data.nrow() then
        count
      else if data[i, 1] = group.0 && data[i, 2] = group.1 then
        loop (succ i, succ count)
      else
        loop (succ i, count)
  in
    loop (0, 0)
  end

fn
count_detections (data   : mtrxszref zero_or_one,
                  group  : @(zero_or_one, zero_or_one),
                  column : intBtwe (3, 6)) : intGte 0 =
  let
    fun
    loop (i     : intGte 0,
          count : intGte 0) : intGte 0 =
      if i2sz i = data.nrow() then
        count
      else if data[i, 1] = group.0 && data[i, 2] = group.1 then
        loop (succ i, count + data[i, column])
      else
        loop (succ i, count)
  in
    loop (0, 0)
  end

fn
correlation_coefficient (data  : mtrxszref zero_or_one,
                         group : @(zero_or_one, zero_or_one))
    : double =
  let
    extern fn sqrt : double -<> double = "mac#sqrt"

    val N = group_size (data, group)
    and NL1 = count_detections (data, group, 3)
    and NL2 = count_detections (data, group, 4)
    and NR1 = count_detections (data, group, 5)
    and NR2 = count_detections (data, group, 6)

    (* Equation (2.4) in the reference. *)
    val sinL = sqrt (g0i2f NL1 / g0i2f N)
    and cosL = sqrt (g0i2f NL2 / g0i2f N)
    and sinR = sqrt (g0i2f NR1 / g0i2f N)
    and cosR = sqrt (g0i2f NR2 / g0i2f N)

    (* Equation (2.3). *)
    val cosLR = (cosR * cosL) + (sinR * sinL)
    and sinLR = (sinR * cosL) - (cosR * sinL)

    (* Equation (2.5). *)
    val kappa = (cosLR * cosLR) - (sinLR * sinLR)
  in
    kappa
  end

fn
analyze_data (data : mtrxszref zero_or_one) : void =
  let
    extern fn sqrt : double -<> double = "mac#sqrt"

    val kappa00 = correlation_coefficient (data, @(0, 0))
    and kappa01 = correlation_coefficient (data, @(0, 1))
    and kappa10 = correlation_coefficient (data, @(1, 0))
    and kappa11 = correlation_coefficient (data, @(1, 1))

    val chsh_contrast = ~kappa00 + kappa01 + kappa10 + kappa11
    and chsh_contrast_nominal = 2.0 * sqrt 2.0

    macdef X (x) = ignoret (,(x)) (* "ignore return value" *)
  in
    X ($extfcall (int, "printf", "\n   light pulse events   %9zu\n",
                  data.nrow()));
    X ($extfcall (int, "printf", "\n    correlation coefs\n"));

    X ($extfcall (int, "printf", "          %4.1f° %4.1f°   %+9.6f\n",
                  angleL1, angleR1, kappa00));
    X ($extfcall (int, "printf", "          %4.1f° %4.1f°   %+9.6f\n",
                  angleL1, angleR2, kappa01));
    X ($extfcall (int, "printf", "          %4.1f° %4.1f°   %+9.6f\n",
                  angleL2, angleR1, kappa10));
    X ($extfcall (int, "printf", "          %4.1f° %4.1f°   %+9.6f\n",
                  angleL2, angleR2, kappa11));

    X ($extfcall (int, "printf", "\n        CHSH contrast   %+9.6f\n",
                  chsh_contrast));
    X ($extfcall (int, "printf", "  2*sqrt(2) = nominal   %+9.6f\n",
                  chsh_contrast_nominal));
    X ($extfcall (int, "printf", "           difference   %+9.6f\n",
                  chsh_contrast - chsh_contrast_nominal));

    X ($extfcall (int, "printf",
                  "\n       CHSH violation   %+9.6f\n\n",
                  chsh_contrast - 2.0))
  end

implement main0 () =
  let
    val data = read_raw_data stdin_ref
  in
    adjust_for_geometry data;
    analyze_data data
  end
</syntaxhighlight>

{{out}}
An example run:
<pre>

   light pulse events      100000

    correlation coefs
           0.0° 22.5°   -0.707806
           0.0° 67.5°   +0.705415
          45.0° 22.5°   +0.712377
          45.0° 67.5°   +0.718882

        CHSH contrast   +2.844480
  2*sqrt(2) = nominal   +2.828427
           difference   +0.016053

       CHSH violation   +0.844480

</pre>

=={{header|Julia}}==
Note that, statistically, the extra correlation in the analysis is inserted via the call to swap_LR_channels(), which 
is statistically dependent on whether the initial photon pair is [0, 90] or [90, 0].
{{trans|Python}}
<syntaxhighlight lang="julia">using Printf

""" The data for one light pulse. """
mutable struct PulseData
    logS::Int
    logL::Int
    logR::Int
    detectedL1::Int
    detectedL2::Int
    detectedR1::Int
    detectedR2::Int
end

""" Swap detectedL1 and detectedL2. """
swap_L_channels(p) = begin p.detectedL1, p.detectedL2 = p.detectedL2, p.detectedL1; p end


""" Swap detectedR1 and detectedR2. """
swap_R_channels(p) = begin p.detectedR1, p.detectedR2 = p.detectedR2, p.detectedR1; p end

"""
Swap channels on both right and left. This is done if the light source was (1,90°)
instead of (1,0°). For in that case the orientations of the polarizing beam splitters,
relative to the light source, is different by 90°.
"""
swap_LR_channels(p::PulseData) = (swap_L_channels(p); swap_R_channels(p))


"""
Split the data into two subsets, according to whether a set item satisfies a predicate.
The return value is a tuple, with those items satisfying the predicate in the first
tuple entry, the other items in the second entry.
"""
split_data(predicate, data) = (filter(predicate, data), filter(!predicate, data))

"""
Some data items are for a (1,0°) light pulse. The others are for a (1,90°) light pulse.
Thus the light pulses are oriented differently with respect to the polarizing beam splitters.
We adjust for that distinction here.
"""
function adjust_data_for_light_pulse_orientation(data)
    data0, data90 = split_data(x -> x.logS == 0, data)
    return vcat(data0, swap_LR_channels.(data90))
end

"""
Split the data into subsets per each arrangement of polarizing beam splitters.
"""
function split_data_according_to_PBS_setting(data)
    dataL1, dataL2 = split_data(x -> x.logL == 0, data)
    dataL1R1, dataL1R2 = split_data(x -> x.logR == 0, dataL1)
    dataL2R1, dataL2R2 = split_data(x -> x.logR == 0, dataL2)
    return (dataL1R1, dataL1R2, dataL2R1, dataL2R2)
end

"""
Compute the correlation coefficient for the subset of the data that corresponding
to a particular setting of the polarizing beam splitters.
"""
function compute_correlation_coefficient(angleL, angleR, data)

    # We make certain the orientations of beam splitters are
    # represented by perpendicular angles in the first and fourth
    # quadrant. This restriction causes no loss of generality, because
    # the orientation of the beam splitter is actually a rotated "X".
    @assert all(0 <= x < 90 for x in (angleL, angleR)) "All orientations are acute angles"

    # Note that the sine is non-negative in Quadrant 1, and the cosine
    # is non-negative in Quadrant 4. Thus we can use the following
    # estimates for cosine and sine. This is Equation (2.4) in the
    # reference. (Note, one can also use Quadrants 1 and 2 and reverse
    # the roles of cosine and sine. And so on like that.)
    N = length(data)
    NL1 = 0
    NL2 = 0
    NR1 = 0
    NR2 = 0
    for item in data
        NL1 += item.detectedL1
        NL2 += item.detectedL2
        NR1 += item.detectedR1
        NR2 += item.detectedR2
    end
    sinL = sqrt(NL1 / N)
    cosL = sqrt(NL2 / N)
    sinR = sqrt(NR1 / N)
    cosR = sqrt(NR2 / N)

    # Now we can apply the reference's Equation (2.3).
    cosLR = (cosR * cosL) + (sinR * sinL)
    sinLR = (sinR * cosL) - (cosR * sinL)

    # And then Equation (2.5).
    kappa = (cosLR * cosLR) - (sinLR * sinLR)

    return kappa
end

""" Read the raw data. Its order does not actually matter. """
function read_raw_data(stream)
    return [PulseData(row...) for row in filter(a -> length(a) == 7,
              [[parse(Int, s) for s in split(line)] for line in readlines(stream)])]
end

function run_analysis(inputstream)
    # Polarizing beam splitter orientations commonly used in actual
    # experiments. These must match the values used in the simulation
    # itself. They are by design all either zero degrees or in the
    # first quadrant.
    anglesL = [0.0, 45.0]
    anglesR = [22.5, 67.5]
    @assert all(0 <= x < 90 for x in [anglesL; anglesR])

    data = read_raw_data(inputstream) |> adjust_data_for_light_pulse_orientation
    dataL1R1, dataL1R2, dataL2R1, dataL2R2 = split_data_according_to_PBS_setting(data)

    kappaL1R1 = compute_correlation_coefficient(anglesL[1], anglesR[1], dataL1R1)
    kappaL1R2 = compute_correlation_coefficient(anglesL[1], anglesR[2], dataL1R2)
    kappaL2R1 = compute_correlation_coefficient(anglesL[2], anglesR[1], dataL2R1)
    kappaL2R2 = compute_correlation_coefficient(anglesL[2], anglesR[2], dataL2R2)

    chsh_contrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2

    # The nominal value of the CHSH contrast for the chosen polarizer
    # orientations is 2*sqrt(2).
    chsh_contrast_nominal = 2 * sqrt(2.0)

    @printf("\n   light pulse events   %9d\n\n", length(data))
    println("    correlation coefs")
    @printf("              %2d° %2d°   %+9.6f\n", anglesL[1], anglesR[1], kappaL1R1)
    @printf("              %2d° %2d°   %+9.6f\n", anglesL[1], anglesR[2], kappaL1R2)
    @printf("              %2d° %2d°   %+9.6f\n", anglesL[2], anglesR[1], kappaL2R1)
    @printf("              %2d° %2d°   %+9.6f\n\n", anglesL[2], anglesR[2], kappaL2R2)
    @printf("        CHSH contrast   %+9.6f\n", chsh_contrast)
    @printf("  2*sqrt(2) = nominal   %+9.6f\n", chsh_contrast_nominal)
    @printf("           difference   %+9.6f\n", chsh_contrast - chsh_contrast_nominal)

    # A "CHSH violation" occurs if the CHSH contrast is >2.
    # https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
    @printf("\n       CHSH violation   %+9.6f\n\n", chsh_contrast - 2)
end

run_analysis("datalog.log")
</syntaxhighlight>{{out}}
An example using data from the Julia implementation of ''[[Simulated optics experiment/Simulator]]''.
<pre>
        
   light pulse events     1000000

    correlation coefs
               0° 22°   -0.459472
               0° 68°   +0.955042
              45° 22°   +0.752309
              45° 68°   +0.749286

        CHSH contrast   +2.913088
  2*sqrt(2) = nominal   +2.828427
           difference   +0.084660

       CHSH violation   +0.913088

</pre>

=={{header|Nim}}==
{{trans|Python}}
<syntaxhighlight lang="Nim">import std/[hashes, math, os, sequtils, sets, strscans, strformat, strutils, sugar]

type
  PulseData = ref object
    logS: int
    logL: int
    logR: int
    detectedL1: int
    detectedL2: int
    detectedR1: int
    detectedR2: int
  PulseDataSet = HashSet[PulseData]


proc hash(pd: PulseData): Hash =
  ## Hash value of a PulseData object (needed for HashSet).
  hash(cast[int](pd[].addr))


func swapLRChannels(pd: PulseData) =
  ## Swap channels on both right and left. This is done if the
  ## light source was (1,90°) instead of (1,0°). For in that case
  ## the orientations of the polarizing beam splitters, relative to
  ## the light source, is different by 90°.
  swap pd.detectedL1, pd.detectedL2
  swap pd.detectedR1, pd.detectedR2


proc split(data: PulseDataSet; predicate: proc(pd: PulseData): bool): tuple[a, b: PulseDataSet] =
  ## Split the data into two subsets, according to whether a set
  ## item satisfies a predicate. The return value is a tuple, with
  ## those items satisfying the predicate in the first tuple entry, the
  ## other items in the second entry.

  let subset1 = collect:
                  for x in data:
                    if predicate(x): {x}
  let subset2 = data - subset1
  result = (subset1, subset2)


proc adjustForLightPulseOrientation(data: PulseDataSet): PulseDataSet =
  ## Some data items are for a (1,0°) light pulse. The others are
  ## for a (1,90°) light pulse. Thus the light pulses are oriented
  ## differently with respect to the polarizing beam splitters. We
  ## adjust for that distinction here.
  let (data0, data90) = data.split(proc(pd: PulseData): bool = pd.logS == 0)
  for item in data90:
    item.swapLRChannels()
  result = union(data0, data90)


proc splitAccordingToPbsSetting(data: PulseDataSet): tuple[a, b, c, d: PulseDataSet] =
  ## Split the data into four subsets: one subset for each
  ## arrangement of the two polarizing beam splitters.
  let (dataL1, dataL2) = data.split(proc(pd: PulseData): bool = pd.logL == 0)
  let (dataL1R1, dataL1R2) = dataL1.split(proc(pd: PulseData): bool = pd.logR == 0)
  let (dataL2R1, dataL2R2) = dataL2.split(proc(pd: PulseData): bool = pd.logR == 0)
  result = (dataL1R1, dataL1R2, dataL2R1, dataL2R2)


func computeCorrelationCoefficient(data: PulseDataSet; angleL, angleR: float): float =
    ## Compute the correlation coefficient for the subset of the data that
    ## corresponding to a particular setting of the polarizing beam splitters.

    # We make certain the orientations of beam splitters are
    # represented by perpendicular angles in the first and fourth
    # quadrant. This restriction causes no loss of generality, because
    # the orientation of the beam splitter is actually a rotated "X".
    assert angleL >= 0 and angleL < 90 and angleR >= 0 and angleR < 90
    #perpendicularL = angleL - 90 # In Quadrant 4.
    #perpendicularR = angleR - 90 # In Quadrant 4.

    # Note that the sine is non-negative in Quadrant 1, and the cosine
    # is non-negative in Quadrant 4. Thus we can use the following
    # estimates for cosine and sine. This is Equation (2.4) in the
    # reference. (Note, one can also use Quadrants 1 and 2 and reverse
    # the roles of cosine and sine. And so on like that.)
    let n = data.len
    var nL1, nL2, nR1, nR2 = 0
    for item in data:
      inc nL1, item.detectedL1
      inc nL2, item.detectedL2
      inc nR1, item.detectedR1
      inc nR2, item.detectedR2
    let sinL = sqrt(nL1 / n)
    let cosL = sqrt(nL2 / n)
    let sinR = sqrt(nR1 / n)
    let cosR = sqrt(nR2 / n)

    # Now we can apply the reference's Equation (2.3).
    let cosLR = (cosR * cosL) + (sinR * sinL)
    let sinLR = (sinR * cosL) - (cosR * sinL)

    # And then Equation (2.5).
    result = (cosLR * cosLR) - (sinLR * sinLR)


proc readRawData(filename: string): PulseDataSet =
  ## Read the raw data. Its order does not actually matter, so we
  ## return the data as a HashSet.

  func makeRecord(line: string): PulseData =
    new(result)
    discard line.scanf("$i $i $i $i $i $i $i", result.logS, result.logL, result.logR,
                       result.detectedL1, result.detectedL2, result.detectedR1, result.detectedR2)

  proc readData(f: File): PulseDataSet =
    let numPulses = f.readline().parseInt()
    for i in 1..numPulses:
      result.incl f.readLine().makeRecord()

  if filename != "-":
    let f = open(filename, fmRead)
    result = f.readData()
    f.close()
  else:
    result = stdin.readData()


when isMainModule:

  if paramCount() notin [0, 1]:
    quit &"Usage: {getAppFilename()} [RAW_DATA_FILE]", QuitFailure

  let rawDataFilename = if paramCount() == 1: paramStr(1) else: "-"

  # Polarizing beam splitter orientations commonly used in actual
  # experiments. These must match the values used in the simulation
  # itself. They are by design all either zero degrees or in the
  # first quadrant.
  const AnglesL = [0.0, 45.0]
  const AnglesR = [22.5, 67.5]
  assert allIt(AnglesL, it >= 0 and it < 90) and allIt(AnglesR, it >= 0 and it < 90)

  var data = readRawData(rawDataFilename)
  data = adjustForLightPulseOrientation(data)
  let (dataL1R1, dataL1R2, dataL2R1, dataL2R2) = data.splitAccordingToPbsSetting()

  let
    kappaL1R1 = dataL1R1.computeCorrelationCoefficient(AnglesL[0], AnglesR[0])
    kappaL1R2 = dataL1R2.computeCorrelationCoefficient(AnglesL[0], AnglesR[1])
    kappaL2R1 = dataL2R1.computeCorrelationCoefficient(AnglesL[1], AnglesR[0])
    kappaL2R2 = dataL2R2.computeCorrelationCoefficient(AnglesL[1], AnglesR[1])

  let chshContrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2

# The nominal value of the CHSH contrast for the chosen polarizer
# orientations is 2*sqrt(2).
let chshContrastNominal = 2 * sqrt(2.0)

echo()
echo &"   light pulse events   {data.len:9}"
echo()
echo "    correlation coefs"
echo &"          {AnglesL[0]:4.1f}° {AnglesR[0]:4.1f}°   {kappaL1R1:+9.6f}"
echo &"          {AnglesL[0]:4.1f}° {AnglesR[1]:4.1f}°   {kappaL1R2:+9.6f}"
echo &"          {AnglesL[1]:4.1f}° {AnglesR[0]:4.1f}°   {kappaL2R1:+9.6f}"
echo &"          {AnglesL[1]:4.1f}° {AnglesR[1]:4.1f}°   {kappaL2R2:+9.6f}"
echo()
echo &"        CHSH contrast   {chshContrast:+9.6f}"
echo &"  2*sqrt(2) = nominal   {chshContrastNominal:+9.6f}"
echo &"           difference   {chshContrast - chshContrastNominal:+9.6f}"

# A "CHSH violation" occurs if the CHSH contrast is >2.
# https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
echo()
echo &"       CHSH violation   {chshContrast - 2:+9.6f}"
echo()
</syntaxhighlight>

{{out}}
The result is identical to that of the Python version. For instance, using data produced by the Python version of the simulator:
<pre>
   light pulse events      100000

    correlation coefs
           0.0° 22.5°   -0.708395
           0.0° 67.5°   +0.705963
          45.0° 22.5°   +0.709600
          45.0° 67.5°   +0.705349

        CHSH contrast   +2.829306
  2*sqrt(2) = nominal   +2.828427
           difference   +0.000879

       CHSH violation   +0.829306</pre>

=={{header|ObjectIcon}}==
<syntaxhighlight lang="objecticon">
#!/bin/env -S oiscript
#
# If you run this program as a script, you might get something like
# the following as a "help" message:
#
#   Usage: /tmp/oixcache/f8f91465ab42189d5623040505fa7f12 [-|raw_data_file]
#
# But you can instead use "oit" to compile this program into a
# bytecode executable, which would print this sort of message instead:
#
#   Usage: ./my_program [-|raw_data_file]
#

#
# Reference:
#
#   A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
#   J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
#   https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)
#

# This is necessary because there are "degree" symbols (°) in the
# code:
$encoding UTF-8

# To get the following more specific set of imports, one can write
# "import io, util" and then run "fiximports" on the source code:
import
   io(BufferStream, FileStream, open, stop),
   ipl.printf(printf),
   util(Math)

$define angleL1 0.0
$define angleL2 45.0
$define angleR1 22.5
$define angleR2 67.5

procedure main (args)
  # "&why" does not exist in regular Icon, although "stop"
  # does.
  #
  # Fellow Gentoo users might notice an uncanny resemblance between
  # the use of "stop" here and the use of "die" in Gentoo ebuilds. In
  # both cases, if we left out the check, the error might be silently
  # ignored.
  case *args of
  {
    0 : sally_forth ("-") | stop (&why)
    1 : sally_forth (args[1]) | stop (&why)
    default :
    {
      io.write ("Usage: ", &progname, " [-|raw_data_file]")
      exit (1)
    }
  }
end

procedure sally_forth (input_file)
  local data, dataL1, dataL2, dataR1, dataR2
  local dataL1R1, dataL1R2, dataL2R1, dataL2R2
  local kappaL1R1, kappaL1R2, kappaL2R1, kappaL2R2
  local chsh_contrast, chsh_contrast_nominal

  data := read_raw_data (input_file)
  data := adjust_for_geometry (data)

  # The following means of dividing the set into four non-overlapping
  # groups has been designed for clarity rather than efficiency. (We
  # should not care about efficiency, anyway, in a program such as
  # this.)
  dataL1 := filter_set (logL_equals_0, data)
  dataL2 := filter_set (logL_equals_1, data)
  dataR1 := filter_set (logR_equals_0, data)
  dataR2 := filter_set (logR_equals_1, data)
  dataL1R1 := dataL1 ** dataR1  # Set intersection.
  dataL1R2 := dataL1 ** dataR2
  dataL2R1 := dataL2 ** dataR1
  dataL2R2 := dataL2 ** dataR2

  # The sizes of the groups should add up to the size of the entire
  # data.
  *dataL1R1 + *dataL1R2 + *dataL2R1 + *dataL2R2 = *data |
    {
      #
      # "&why" does not exist in regular Icon.
      #
      # In all the dialects of Icon, however, it is common for failure
      # of a procedure to indicate an error.
      #
      # This has a severe downside that, in many contexts, a failure
      # is silently ignored. This behavior is unlike an exception
      # (which is either caught or it terminates the
      # program). However, in some ways it is like the practice of
      # indicating an error through a return value. Differences from
      # that include: (a) failure will often cause the containing
      # context to fail as well (possibly resulting in assignments and
      # swaps/exchanges being undone!), and (b) there is no return
      # value in case of failure.
      #
      # BTW Object Icon does have exceptions:
      # https://chemoelectric.github.io/objecticon/ExceptionPackage.html
      #
      &why := "The four groups of data have been formed incorrectly."
      fail
    }

  kappaL1R1 := correlation_coefficient (angleL1, angleR1, dataL1R1)
  kappaL1R2 := correlation_coefficient (angleL1, angleR2, dataL1R2)
  kappaL2R1 := correlation_coefficient (angleL2, angleR1, dataL2R1)
  kappaL2R2 := correlation_coefficient (angleL2, angleR2, dataL2R2)

  chsh_contrast := (-kappaL1R1) + kappaL1R2 + kappaL2R1 + kappaL2R2

  printf ("\n   light pulse events   %9d\n", *data)
  printf ("\n    correlation coefs\n")

  # The u"" is necessary in Object Icon when the string is a Unicode
  # string.
  printf (u"          %4.1r° %4.1r°   %+9.6r\n",
          angleL1, angleR1, kappaL1R1)
  printf (u"          %4.1r° %4.1r°   %+9.6r\n",
          angleL1, angleR2, kappaL1R2)
  printf (u"          %4.1r° %4.1r°   %+9.6r\n",
          angleL2, angleR1, kappaL2R1)
  printf (u"          %4.1r° %4.1r°   %+9.6r\n",
          angleL2, angleR2, kappaL2R2)

  chsh_contrast_nominal := 2.0 * Math.sqrt (2.0)

  printf ("\n        CHSH contrast   %+9.6r\n", chsh_contrast)
  printf ("  2*sqrt(2) = nominal   %+9.6r\n", chsh_contrast_nominal)
  printf ("           difference   %+9.6r\n", (chsh_contrast -
                                             chsh_contrast_nominal))

  # A "CHSH violation" occurs if the CHSH contrast is >2.
  # https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
  #
  # Note CAREFULLY, however, that we are talking here about simulation
  # of a CLASSICAL optics model, using the CLASSICAL formula for
  # correlation coefficients (of the respective electromagnetic
  # fields). We are calculating those coefficients by inference from
  # detection events.
  #
  # Some might try to object that a CHSH contrast cannot be calculated
  # from correlation coefficients computed according to the classical
  # formula. But that would be assuming the conclusion: OF COURSE no
  # classical model can have a CHSH violation, if by fiat you do not
  # let a classical model and the classical formulas be used!
  #
  printf ("\n       CHSH violation   %+9.6r\n\n", (chsh_contrast - 2))

  # The following is necessary (it is equivalent to "return &null"),
  # because we are using failure to indicate an error.
  return
end

procedure correlation_coefficient (angleL, angleR, data_group)
  local N, NL1, NL2, NR1, NR2, item
  local cosL, sinL, cosR, sinR, cosLR, sinLR, kappa

  # We know angleL and angleR are in Quadrant 1. Its two
  # perpendiculars are thus in Quadrants 2 and 4. We can choose
  # whichever of the two we want. Choose Quadrant 4. The sine is
  # non-negative in Quadrant 1, and the cosine is non-negative in
  # Quadrant 4. Thus we can use the following estimates for cosine and
  # sine, and not have to deal with negative signs.
  #
  # This is Equation (2.4) in the reference.
  #
  N := real (*data_group)
  NL1 := 0.0
  NL2 := 0.0
  NR1 := 0.0
  NR2 := 0.0
  every item := !data_group do
  {
    NL1 +:= item.detectedL1()
    NL2 +:= item.detectedL2()
    NR1 +:= item.detectedR1()
    NR2 +:= item.detectedR2()
  }
  sinL := Math.sqrt (NL1 / N)
  cosL := Math.sqrt (NL2 / N)
  sinR := Math.sqrt (NR1 / N)
  cosR := Math.sqrt (NR2 / N)

  # Apply the reference's Equation (2.3).
  cosLR := (cosR * cosL) + (sinR * sinL)
  sinLR := (sinR * cosL) - (cosR * sinL)

  # Apply Equation (2.5).
  kappa := (cosLR * cosLR) - (sinLR * sinLR)

  return kappa
end

procedure read_raw_data (input_file)
  local data, f
  local num_lines
  local logS, logL, logR
  local detectedL1, detectedL2
  local detectedR1, detectedR2

  if type (input_file) ~== "string" then
  {
    f := input_file
    num_lines := integer (f.read())
    data := set ()

    # The "every" construct below goes through all the numbers
    # generated by "1 to num_lines", after which it quits because the
    # generator FAILS. (I leave out any loop variable, "every i :=",
    # because one is not needed here.) One might, especially at first,
    # accidentally write something like "while i := 1 to num_lines",
    # but this would actually be a forever-loop, with i always equal
    # to 1. The "while" construct, unlike "every", would restart the
    # generator from scratch each time, so that the generator ALWAYS
    # SUCCEEDS. But "while" is what you would want for something like
    # "while line := f.read()", because "read" is not a generator. It
    # is a procedure that either returns a line of input or fails.
    every 1 to num_lines do
    {
      # The following is a "string-scanning environment". It is a
      # series of operations on the string returned by f.read(). See
      # also
      # https://en.wikipedia.org/w/index.php?title=Icon_(programming_language)&oldid=1146573004
      #
      # Object Icon is not mentioned in that article, but right now I
      # do not feel like editing the Wikipedia to mention it.
      #
      # It IS mentioned that Python generators are supposedly inspired
      # by Icon generators (and presumably also by co-expressions). I
      # do not know if that is true, but there are similarities:
      # Python "yield" is similar to Icon "suspend". (Scheme’s
      # "call/cc" is far more general, however, and can be used to
      # implement yield/suspend in a few lines of code.)
      f.read() ?
      {
        tab(many(' \t'))
        logS := integer (tab(many(&digits)))
        tab(many(' \t'))
        logL := integer (tab(many(&digits)))
        tab(many(' \t'))
        logR := integer (tab(many(&digits)))
        tab(many(' \t'))
        detectedL1 := integer (tab(many(&digits)))
        tab(many(' \t'))
        detectedL2 := integer (tab(many(&digits)))
        tab(many(' \t'))
        detectedR1 := integer (tab(many(&digits)))
        tab(many(' \t'))
        detectedR2 := integer (tab(many(&digits)))
      }
      insert (data, PulseEventData (logS, logL, logR,
                                    detectedL1, detectedL2,
                                    detectedR1, detectedR2))
    }
  }
  else if input_file == "-" then
    data := read_raw_data (BufferStream (FileStream.stdin))
  else
  {
    f := open (input_file, "r") | stop (&why)
    data := read_raw_data (f) | stop (&why)
    f.close ()
  }
  return data
end

procedure adjust_for_geometry (data)
  local S, S1, S2, S2a, item
  S := data
  S1 := filter_set (logS_equals_0, S)
  S2 := S -- S1
  S2a := set ()
  every item := !S2 do
    #
    # Swap channels, to account for geometry.
    #
    # This step has been objected to as introducing "quantum magic",
    # but note that we are doing a CLASSICAL data analysis, having
    # nought to do with quantum mechanics. The light pulse source has
    # recorded for us which type of light pulses were emitted. One
    # might or might not be able to imagine a device for doing this
    # with "photons", but we do not need to: (a) this is a SIMULATION,
    # so we CAN use information we do not KNOW HOW to get, as long as
    # the information is moving strictly forwards in time; besides
    # which, (b) we can INSTEAD imagine a device that randomly selects
    # from one of two pairs of light beams, and which forms pulses by
    # opening and closing a shutter.
    #
    # According to "Bell's Theorem", it is impossible for such a
    # classical model--one that works solely by contact action moving
    # forwards in time--to produce the results we get. In my opinion,
    # this simulation should be enough to end the entire debacle this
    # "theorem" has helped create. But, in fact, it is easy to show
    # that John Bell's "proof" does the equivalent of using the
    # Angle-Side-Side "theorem" in geometry class. It is fake
    # mathematics. (See the "Talk" page of the task, where I
    # explain Bell's mistake in just the few words it takes.)
    #
    # One might or might not have noticed we are also including
    # non-detection events in denominators, even though they are not
    # part of the measured data set, in experiments that do not record
    # when pulses were emitted in the first place. This is a
    # SIMULATION, however, so we CAN include them.
    #
    #    -- Chemoelectric, 30 May 2023
    #
    # In any case, we are demonstrating how to use sets in Object
    # Icon. Sets are very similar in regular Icon and in Unicon,
    # although regular Icon does not have "classes". (It DOES have
    # "records".)
    #
    insert (S2a, PulseEventData (item.logS(),
                                 item.logL(), item.logR(),
                                 item.detectedL2(),
                                 item.detectedL1(),
                                 item.detectedR2(),
                                 item.detectedR1()))
  return S1 ++ S2a
end

procedure filter_set (semidet_predicate, S)
  # The predicate is "semi-deterministic". It does not return a
  # boolean. Instead, it either succeeds (and in Object Icon has to
  # return something, which we ignore) or it fails. You will see the
  # same approach to decision-making in Prolog and Mercury programs.
  local S1, item
  S1 := set ()
  every item := !S do
    if semidet_predicate (item) then
      insert (S1, item)
  return S1
end

procedure logS_equals_0 (event_data)
  # This is a semidet predicate. It succeeds if logS is zero,
  # otherwise it fails. The return value (which is returned only on
  # success and is always 0) will be ignored.
  return event_data.logS() = 0
end

#######################################
# Some more semidet predicates:

procedure logL_equals_0 (event_data)
  return event_data.logL() = 0
end

procedure logL_equals_1 (event_data)
  return event_data.logL() = 1
end

procedure logR_equals_0 (event_data)
  return event_data.logR() = 0
end

procedure logR_equals_1 (event_data)
  return event_data.logR() = 1
end

#######################################

final class PulseEventData ()
  # I have tried here to make PulseEventData objects both immutable
  # and unextendable, so no one can doubt the absence of "clever
  # tricks". They are merely containers of data: each of them unique,
  # so they can be stored in sets.

  private var_logS, var_logL, var_logR
  private var_detectedL1, var_detectedL2
  private var_detectedR1, var_detectedR2

  public new (logS, logL, logR,
              detectedL1, detectedL2,
              detectedR1, detectedR2)
    var_logS := logS
    var_logL := logL
    var_logR := logR
    var_detectedL1 := detectedL1
    var_detectedL2 := detectedL2
    var_detectedR1 := detectedR1
    var_detectedR2 := detectedR2
    return
  end

  public logS ()
    return var_logS
  end

  public logL ()
    return var_logL
  end

  public logR ()
    return var_logR
  end

  public detectedL1 ()
    return var_detectedL1
  end

  public detectedL2 ()
    return var_detectedL2
  end

  public detectedR1 ()
    return var_detectedR1
  end

  public detectedR2 ()
    return var_detectedR2
  end

  public to_string ()
    return ("PulseEventData(" || logS() ||
            ", " || logL() ||
            ", " || logR() ||
            ", " || detectedL1() ||
            ", " || detectedL2() ||
            ", " || detectedR1() ||
            ", " || detectedR2() || ")")
  end
end
</syntaxhighlight>

{{out}}
Using data generated by a 10000-point run of the Python version of the simulation.

Please be aware that one run of that Python program will not produce the same output as another. There are two reasons: (1) the random number generator is seeded randomly, and (2) the simulated devices are running as separate processes. (How Python handles random numbers in a program running as multiple processes I haven't the vaguest idea.)

A standard set of raw data might be published for people to use. Then we could always get the same numbers. But, to quote Cyrano Jones, ''"What would happen to Man's quest for knowledge?"''

<pre>

   light pulse events       10000

    correlation coefs
           0.0° 22.5°   -0.699476
           0.0° 67.5°   +0.714172
          45.0° 22.5°   +0.692368
          45.0° 67.5°   +0.711478

        CHSH contrast   +2.817493
  2*sqrt(2) = nominal   +2.828427
           difference   -0.010935

       CHSH violation   +0.817493

</pre>

=={{header|Phix}}==
Gone for as simple as possible, includes the generation phase, and for javascript compatability file save/load commented out.
<!--<syntaxhighlight lang="phix">(phixonline)-->
 <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
 <span style="color: #004080;">sequence</span> <span style="color: #000000;">datalog</span>
 <span style="color: #008080;">enum</span> <span style="color: #000000;">S</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">L</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">R</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">L1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">L2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">R1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">R2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">LL</span><span style="color: #0000FF;">=$</span>
 <span style="color: #008080;">constant</span> <span style="color: #000000;">angles</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- (this is indexed via L and R)</span>
                    <span style="color: #0000FF;">{</span><span style="color: #000000;">0.0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">45.0</span><span style="color: #0000FF;">},</span>
                    <span style="color: #0000FF;">{</span><span style="color: #000000;">22.5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">67.5</span><span style="color: #0000FF;">}}</span>
 
 <span style="color: #008080;">procedure</span> <span style="color: #000000;">detector_recieve</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">lr</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">amplitude</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">atom</span> <span style="color: #000000;">intensity</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">amplitude</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">amplitude</span><span style="color: #0000FF;">,</span>
          <span style="color: #000000;">randnum</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()</span>
     <span style="color: #000000;">datalog</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">:=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">randnum</span> <span style="color: #0000FF;"><=</span> <span style="color: #000000;">intensity</span><span style="color: #0000FF;">)</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 
 <span style="color: #008080;">procedure</span> <span style="color: #000000;">splitter_receive</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">lr</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">source_angle</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">integer</span> <span style="color: #000000;">rand_0_or_1</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">rand_range</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
     <span style="color: #000000;">datalog</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">rand_0_or_1</span>
     <span style="color: #004080;">atom</span> <span style="color: #000000;">my_angle</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">angles</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">][</span><span style="color: #000000;">rand_0_or_1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span>
          <span style="color: #000000;">relative_angle</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">my_angle</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">source_angle</span><span style="color: #0000FF;">,</span>
          <span style="color: #000000;">amplitude1</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">180</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">relative_angle</span><span style="color: #0000FF;">)),</span>
          <span style="color: #000000;">amplitude2</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">180</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">relative_angle</span><span style="color: #0000FF;">))</span>
     <span style="color: #000000;">detector_recieve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">amplitude1</span><span style="color: #0000FF;">)</span>
     <span style="color: #000000;">detector_recieve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">amplitude2</span><span style="color: #0000FF;">)</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 
 <span style="color: #008080;">procedure</span> <span style="color: #000000;">light_source</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">num_events</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">num_events</span> <span style="color: #008080;">do</span>
         <span style="color: #004080;">integer</span> <span style="color: #000000;">rand_0_or_1</span> <span style="color: #0000FF;">:=</span> <span style="color: #7060A8;">rand_range</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
         <span style="color: #000000;">datalog</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">S</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">rand_0_or_1</span>
         <span style="color: #008080;">if</span> <span style="color: #000000;">rand_0_or_1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
             <span style="color: #000000;">splitter_receive</span><span style="color: #0000FF;">(</span><span style="color: #000000;">L</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.0</span><span style="color: #0000FF;">)</span>
             <span style="color: #000000;">splitter_receive</span><span style="color: #0000FF;">(</span><span style="color: #000000;">R</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">90.0</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">else</span>
             <span style="color: #000000;">splitter_receive</span><span style="color: #0000FF;">(</span><span style="color: #000000;">L</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">90.0</span><span style="color: #0000FF;">)</span>
             <span style="color: #000000;">splitter_receive</span><span style="color: #0000FF;">(</span><span style="color: #000000;">R</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.0</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 
 <span style="color: #008080;">procedure</span> <span style="color: #000000;">generate</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">num_events</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">object</span> <span style="color: #000000;">output_file</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
     <span style="color: #000000;">datalog</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span><span style="color: #000000;">LL</span><span style="color: #0000FF;">),</span><span style="color: #000000;">num_events</span><span style="color: #0000FF;">)</span>
     <span style="color: #000000;">light_source</span><span style="color: #0000FF;">(</span><span style="color: #000000;">num_events</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">if</span> <span style="color: #004080;">string</span><span style="color: #0000FF;">(</span><span style="color: #000000;">output_file</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">or</span> <span style="color: #000000;">num_events</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">20</span> <span style="color: #008080;">then</span>
         <span style="color: #004080;">integer</span> <span style="color: #000000;">fn</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">iff</span><span style="color: #0000FF;">(</span><span style="color: #004080;">string</span><span style="color: #0000FF;">(</span><span style="color: #000000;">output_file</span><span style="color: #0000FF;">)?</span><span style="color: #7060A8;">open</span><span style="color: #0000FF;">(</span><span style="color: #000000;">output_file</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"w"</span><span style="color: #0000FF;">):</span><span style="color: #000000;">output_file</span><span style="color: #0000FF;">)</span>
         <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fn</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%d\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">num_events</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">l</span> <span style="color: #008080;">in</span> <span style="color: #000000;">datalog</span> <span style="color: #008080;">do</span>
             <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fn</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">:=</span><span style="color: #008000;">"%d"</span><span style="color: #0000FF;">))</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #008080;">if</span> <span style="color: #004080;">string</span><span style="color: #0000FF;">(</span><span style="color: #000000;">output_file</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">close</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fn</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 <span style="color: #000000;">generate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">100000</span><span style="color: #0000FF;">)</span>
 <span style="color: #000080;font-style:italic;">--generate(1e5,"datalog.log")</span>
 
 <span style="color: #008080;">function</span> <span style="color: #000000;">adjust_data_for_light_pulse_orientation</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">data</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span> <span style="color: #008080;">in</span> <span style="color: #000000;">data</span> <span style="color: #008080;">do</span>
         <span style="color: #008080;">if</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">S</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
             <span style="color: #000000;">data</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">reinstate</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">L1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">L2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R2</span><span style="color: #0000FF;">},</span><span style="color: #7060A8;">extract</span><span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">L2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">L1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R1</span><span style="color: #0000FF;">}))</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">data</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
 <span style="color: #008080;">function</span> <span style="color: #000000;">split_data_according_to_PBS_setting</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">data</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">sequence</span> <span style="color: #000000;">dataL1R1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{},</span> <span style="color: #000000;">dataL1R2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{},</span>
              <span style="color: #000000;">dataL2R1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{},</span> <span style="color: #000000;">dataL2R2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">d</span> <span style="color: #008080;">in</span> <span style="color: #000000;">data</span> <span style="color: #008080;">do</span>
         <span style="color: #004080;">sequence</span> <span style="color: #000000;">dLR</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">L</span><span style="color: #0000FF;">..</span><span style="color: #000000;">R</span><span style="color: #0000FF;">]</span>
         <span style="color: #008080;">if</span>    <span style="color: #000000;">dLR</span><span style="color: #0000FF;">={</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">then</span> <span style="color: #000000;">dataL1R1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL1R1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">elsif</span> <span style="color: #000000;">dLR</span><span style="color: #0000FF;">={</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">then</span> <span style="color: #000000;">dataL1R2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL1R2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">elsif</span> <span style="color: #000000;">dLR</span><span style="color: #0000FF;">={</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">then</span> <span style="color: #000000;">dataL2R1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL2R1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">elsif</span> <span style="color: #000000;">dLR</span><span style="color: #0000FF;">={</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">then</span> <span style="color: #000000;">dataL2R2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL2R2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
         <span style="color: #008080;">else</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">dataL1R1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dataL1R2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dataL2R1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dataL2R2</span><span style="color: #0000FF;">}</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
 <span style="color: #008080;">function</span> <span style="color: #000000;">compute_correlation_coefficient</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">data</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">integer</span> <span style="color: #000000;">N</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">data</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">NL1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">NL2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">NR1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">NR2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">item</span> <span style="color: #008080;">in</span> <span style="color: #000000;">data</span> <span style="color: #008080;">do</span>
         <span style="color: #000000;">NL1</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">item</span><span style="color: #0000FF;">[</span><span style="color: #000000;">L1</span><span style="color: #0000FF;">]</span>
         <span style="color: #000000;">NL2</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">item</span><span style="color: #0000FF;">[</span><span style="color: #000000;">L2</span><span style="color: #0000FF;">]</span>
         <span style="color: #000000;">NR1</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">item</span><span style="color: #0000FF;">[</span><span style="color: #000000;">R1</span><span style="color: #0000FF;">]</span>
         <span style="color: #000000;">NR2</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">item</span><span style="color: #0000FF;">[</span><span style="color: #000000;">R2</span><span style="color: #0000FF;">]</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #004080;">atom</span> <span style="color: #000000;">sinL</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">NL1</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">N</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">cosL</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">NL2</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">N</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">sinR</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">NR1</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">N</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">cosR</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">NR2</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">N</span><span style="color: #0000FF;">),</span>
         <span style="color: #000000;">cosLR</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">cosR</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">cosL</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">sinR</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">sinL</span><span style="color: #0000FF;">),</span>
         <span style="color: #000000;">sinLR</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">sinR</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">cosL</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">cosR</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">sinL</span><span style="color: #0000FF;">),</span>
         <span style="color: #000000;">kappa</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">cosLR</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">cosLR</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">sinLR</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">sinLR</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">kappa</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
 <span style="color: #000080;font-style:italic;">--function read_from_file(string input_file)
 --  sequence data = read_lines(input_file)
 --  assert(to_integer(data[1])+1=length(data))
 --  data = data[2..$]
 --  for i,d in data do
 --      data[i] = scanf(d,"%d %d %d %d %d %d %d")[1]
 --  end for
 --  return data
 --end function</span>
 
 <span style="color: #008080;">procedure</span> <span style="color: #000000;">run_analysis</span><span style="color: #0000FF;">()</span>
 <span style="color: #000080;font-style:italic;">--procedure run_analysis(string input_file)
 --  sequence data = read_from_file(input_file)
 --  assert(data=datalog)</span>
     <span style="color: #004080;">sequence</span> <span style="color: #000000;">data</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">datalog</span>
     <span style="color: #000000;">datalog</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> <span style="color: #000080;font-style:italic;">-- (kill refcount)</span>
     <span style="color: #000000;">data</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">adjust_data_for_light_pulse_orientation</span><span style="color: #0000FF;">(</span><span style="color: #000000;">data</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">sequence</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">dataL1R1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dataL1R2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dataL2R1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">dataL2R2</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">split_data_according_to_PBS_setting</span><span style="color: #0000FF;">(</span><span style="color: #000000;">data</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{{},{</span><span style="color: #000000;">angleL1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">angleL2</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">angleR1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">angleR2</span><span style="color: #0000FF;">}}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">angles</span><span style="color: #0000FF;">,</span>
          <span style="color: #000000;">kappaL1R1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">compute_correlation_coefficient</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL1R1</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">kappaL1R2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">compute_correlation_coefficient</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL1R2</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">kappaL2R1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">compute_correlation_coefficient</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL2R1</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">kappaL2R2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">compute_correlation_coefficient</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dataL2R2</span><span style="color: #0000FF;">),</span>
      <span style="color: #000000;">chsh_contrast</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">kappaL1R1</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">kappaL1R2</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">kappaL2R1</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">kappaL2R2</span><span style="color: #0000FF;">,</span>
     <span style="color: #000000;">chsh_contrast_nominal</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2.0</span><span style="color: #0000FF;">)</span>
 
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n   light pulse events   %9d\n\n"</span><span style="color: #0000FF;">,</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">data</span><span style="color: #0000FF;">))</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"    correlation coefs\n"</span><span style="color: #0000FF;">)</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"          %4.1f° %4.1f°   %+9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">angleL1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">angleR1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">kappaL1R1</span><span style="color: #0000FF;">})</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"          %4.1f° %4.1f°   %+9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">angleL1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">angleR2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">kappaL1R2</span><span style="color: #0000FF;">})</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"          %4.1f° %4.1f°   %+9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">angleL2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">angleR1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">kappaL2R1</span><span style="color: #0000FF;">})</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"          %4.1f° %4.1f°   %+9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">angleL2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">angleR2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">kappaL2R2</span><span style="color: #0000FF;">})</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"        CHSH contrast   %+9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">chsh_contrast</span><span style="color: #0000FF;">)</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"  2*sqrt(2) = nominal   %+9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">chsh_contrast_nominal</span><span style="color: #0000FF;">)</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"           difference   %+9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">chsh_contrast</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">chsh_contrast_nominal</span><span style="color: #0000FF;">)</span>
     <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n       CHSH violation   %+9.6f\n\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">chsh_contrast</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
 
 <span style="color: #000000;">run_analysis</span><span style="color: #0000FF;">()</span>
 <span style="color: #000080;font-style:italic;">--run_analysis("datalog.log")</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
   light pulse events      100000

    correlation coefs
           0.0° 22.5°   -0.708021
           0.0° 67.5°   +0.700617
          45.0° 22.5°   +0.709869
          45.0° 67.5°   +0.712739

        CHSH contrast   +2.831245
  2*sqrt(2) = nominal   +2.828427
           difference   +0.002818

       CHSH violation   +0.831245
</pre>

=={{header|Python}}==

The program takes one optional command line argument: the filename of the raw data saved by one of the ''[[Simulated optics experiment/Simulator]]'' implementations. If the argument is omitted or is "-", the raw data is read from standard input.

<syntaxhighlight lang="python">
#!/bin/env -S python3

#
# Reference:
#
#   A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
#   J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
#   https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)
#

import sys
from math import atan2, cos, floor, pi, radians, sin, sqrt

class PulseData():
    '''The data for one light pulse.'''

    def __init__(self, logS, logL, logR, detectedL1, detectedL2,
                 detectedR1, detectedR2):
        self.logS = logS
        self.logL = logL
        self.logR = logR
        self.detectedL1 = detectedL1
        self.detectedL2 = detectedL2
        self.detectedR1 = detectedR1
        self.detectedR2 = detectedR2

    def __repr__(self):
        return ("PulseData(" + str(self.logS) + ", "
                + str(self.logL) + ", "
                + str(self.logR) + ", "
                + str(self.detectedL1) + ", "
                + str(self.detectedL2) + ", "
                + str(self.detectedR1) + ", "
                + str(self.detectedR2) + ")")

    def swap_L_channels(self):
        '''Swap detectedL1 and detectedL2.'''
        (self.detectedL1, self.detectedL2) = \
            (self.detectedL2, self.detectedL1)

    def swap_R_channels(self):
        '''Swap detectedR1 and detectedR2.'''
        (self.detectedR1, self.detectedR2) = \
            (self.detectedR2, self.detectedR1)

    def swap_LR_channels(self):
        '''Swap channels on both right and left. This is done if the
        light source was (1,90°) instead of (1,0°). For in that case
        the orientations of the polarizing beam splitters, relative to
        the light source, is different by 90°.'''
        self.swap_L_channels()
        self.swap_R_channels()

def split_data(predicate, data):
    '''Split the data into two subsets, according to whether a set
    item satisfies a predicate. The return value is a tuple, with
    those items satisfying the predicate in the first tuple entry, the
    other items in the second entry.

    '''
    subset1 = {x for x in data if predicate(x)}
    subset2 = data - subset1
    return (subset1, subset2)

def adjust_data_for_light_pulse_orientation(data):
    '''Some data items are for a (1,0°) light pulse. The others are
    for a (1,90°) light pulse. Thus the light pulses are oriented
    differently with respect to the polarizing beam splitters. We
    adjust for that distinction here.'''
    data0, data90 = split_data(lambda item: item.logS == 0, data)
    for item in data90:
        item.swap_LR_channels()
    return (data0 | data90)     # Set union.

def split_data_according_to_PBS_setting(data):
    '''Split the data into four subsets: one subset for each
    arrangement of the two polarizing beam splitters.'''
    dataL1, dataL2 = split_data(lambda item: item.logL == 0, data)
    dataL1R1, dataL1R2 = \
        split_data(lambda item: item.logR == 0, dataL1)
    dataL2R1, dataL2R2 = \
        split_data(lambda item: item.logR == 0, dataL2)
    return (dataL1R1, dataL1R2, dataL2R1, dataL2R2)

def compute_correlation_coefficient(angleL, angleR, data):
    '''Compute the correlation coefficient for the subset of the data
    that corresponding to a particular setting of the polarizing beam
    splitters.'''

    # We make certain the orientations of beam splitters are
    # represented by perpendicular angles in the first and fourth
    # quadrant. This restriction causes no loss of generality, because
    # the orientation of the beam splitter is actually a rotated "X".
    assert (all(0 <= x and x < 90 for x in (angleL, angleR)))
    #perpendicularL = angleL - 90 # In Quadrant 4.
    #perpendicularR = angleR - 90 # In Quadrant 4.

    # Note that the sine is non-negative in Quadrant 1, and the cosine
    # is non-negative in Quadrant 4. Thus we can use the following
    # estimates for cosine and sine. This is Equation (2.4) in the
    # reference. (Note, one can also use Quadrants 1 and 2 and reverse
    # the roles of cosine and sine. And so on like that.)
    N = len(data)
    NL1 = 0
    NL2 = 0
    NR1 = 0
    NR2 = 0
    for item in data:
        NL1 += item.detectedL1
        NL2 += item.detectedL2
        NR1 += item.detectedR1
        NR2 += item.detectedR2
    sinL = sqrt(NL1 / N)
    cosL = sqrt(NL2 / N)
    sinR = sqrt(NR1 / N)
    cosR = sqrt(NR2 / N)

    # Now we can apply the reference's Equation (2.3).
    cosLR = (cosR * cosL) + (sinR * sinL)
    sinLR = (sinR * cosL) - (cosR * sinL)

    # And then Equation (2.5).
    kappa = (cosLR * cosLR) - (sinLR * sinLR)

    return kappa

def read_raw_data(filename):
    '''Read the raw data. Its order does not actually matter, so we
    return the data as a set.'''

    def make_record(line):
        x = line.split()
        record = PulseData(logS = int(x[0]),
                           logL = int(x[1]),
                           logR = int(x[2]),
                           detectedL1 = int(x[3]),
                           detectedL2 = int(x[4]),
                           detectedR1 = int(x[5]),
                           detectedR2 = int(x[6]))
        return record

    def read_data(f):
        num_pulses = int(f.readline())
        data = set()
        for i in range(num_pulses):
            data.add(make_record(f.readline()))
        return data

    if filename != "-":
        with open(filename, "r") as f:
            data = read_data(f)
    else:
        data = read_data(sys.stdin)
    return data

if __name__ == '__main__':

    if len(sys.argv) != 1 and len(sys.argv) != 2:
        print("Usage: " + sys.argv[0] + " [RAW_DATA_FILE]")
        sys.exit(1)
    raw_data_filename = (sys.argv[1] if len(sys.argv) == 2 else "-")

    # Polarizing beam splitter orientations commonly used in actual
    # experiments. These must match the values used in the simulation
    # itself. They are by design all either zero degrees or in the
    # first quadrant.
    anglesL = (0.0, 45.0)
    anglesR = (22.5, 67.5)
    assert (all(0 <= x and x < 90 for x in anglesL + anglesR))

    data = read_raw_data(raw_data_filename)
    data = adjust_data_for_light_pulse_orientation(data)
    dataL1R1, dataL1R2, dataL2R1, dataL2R2 = \
        split_data_according_to_PBS_setting(data)

    kappaL1R1 = \
        compute_correlation_coefficient (anglesL[0], anglesR[0],
                                         dataL1R1)
    kappaL1R2 = \
        compute_correlation_coefficient (anglesL[0], anglesR[1],
                                         dataL1R2)
    kappaL2R1 = \
        compute_correlation_coefficient (anglesL[1], anglesR[0],
                                         dataL2R1)
    kappaL2R2 = \
        compute_correlation_coefficient (anglesL[1], anglesR[1],
                                         dataL2R2)

    chsh_contrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2

    # The nominal value of the CHSH contrast for the chosen polarizer
    # orientations is 2*sqrt(2).
    chsh_contrast_nominal = 2 * sqrt(2.0)

    print()
    print("   light pulse events   %9d" % len(data))
    print()
    print("    correlation coefs")
    print("          %4.1f° %4.1f°   %+9.6f" % (anglesL[0], anglesR[0],
                                                kappaL1R1))
    print("          %4.1f° %4.1f°   %+9.6f" % (anglesL[0], anglesR[1],
                                                kappaL1R2))
    print("          %4.1f° %4.1f°   %+9.6f" % (anglesL[1], anglesR[0],
                                                kappaL2R1))
    print("          %4.1f° %4.1f°   %+9.6f" % (anglesL[1], anglesR[1],
                                                kappaL2R2))
    print()
    print("        CHSH contrast   %+9.6f" % chsh_contrast)
    print("  2*sqrt(2) = nominal   %+9.6f" % chsh_contrast_nominal)
    print("           difference   %+9.6f" % (chsh_contrast -
                                              chsh_contrast_nominal))

    # A "CHSH violation" occurs if the CHSH contrast is >2.
    # https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
    print()
    print("       CHSH violation   %+9.6f" % (chsh_contrast - 2))
    print()
</syntaxhighlight>

{{out}}
An example using data from the Python implementation of ''[[Simulated optics experiment/Simulator]]''.
<pre>

   light pulse events      100000

    correlation coefs
           0.0° 22.5°   -0.707806
           0.0° 67.5°   +0.705415
          45.0° 22.5°   +0.712377
          45.0° 67.5°   +0.718882

        CHSH contrast   +2.844480
  2*sqrt(2) = nominal   +2.828427
           difference   +0.016053

       CHSH violation   +0.844480

</pre>

=={{header|RATFOR}}==
There is a Ratfor preprocessor in C at https://sourceforge.net/p/chemoelectric/ratfor77 
It outputs FORTRAN 77, which you can then compile like any FORTRAN 77 program.

Ratfor itself is a mix of Fortran syntax and C syntax. It was once, in a sense, the "Unix dialect of FORTRAN". (It was expanded into "Extended Fortran Language" but I do not remember what that added. I encountered EFL in System V Release 3.)

This program does the analysis in what might be the most efficient way possible. It uses very small constant space.

<syntaxhighlight lang="ratfor">
# -*- mode: indented-text; tab-width: 2; -*-

# Reference:
#
#  A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
#  J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
#  https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)

# Standard FORTRAN 77 provides no way to allocate space at runtime.
# One would have to make a call to C, or do some such thing, or simply
# recompile if you need a bigger array. However, it is not
# necessary. This program computes the sums as it reads the data in.

# Once upon a time, because memory was precious and compilers made it
# hard to allocate memory, this would have been THE way to do such
# calculations. But the code ends up being less modular.

# In 2023 there is, for this particular program, no fathomable reason
# not to use double precision. So I shall use double precision.

define(angleL1, 0.0D0)          # 'D' for double precision.
define(angleL2, 45.0D0)
define(angleR1, 22.5D0)
define(angleR2, 67.5D0)

subroutine count (N00, N01, N10, N11,
                  N00L1, N00L2, N00R1, N00R2,
                  N01L1, N01L2, N01R1, N01R2,
                  N10L1, N10L2, N10R1, N10R2,
                  N11L1, N11L2, N11R1, N11R2)
  implicit none
  integer N00, N01, N10, N11
  integer N00L1, N00L2, N00R1, N00R2
  integer N01L1, N01L2, N01R1, N01R2
  integer N10L1, N10L2, N10R1, N10R2
  integer N11L1, N11L2, N11R1, N11R2
  integer fields(1:7)
  integer events, i, tmp
  continue
  N00 = 0;  N01 = 0;  N10 = 0;  N11 = 0
  N00L1 = 0;  N00L2 = 0;  N00R1 = 0;  N00R2 = 0
  N01L1 = 0;  N01L2 = 0;  N01R1 = 0;  N01R2 = 0
  N10L1 = 0;  N10L2 = 0;  N10R1 = 0;  N10R2 = 0
  N11L1 = 0;  N11L2 = 0;  N11R1 = 0;  N11R2 = 0
  read (*,*) events
  for (i = 0; i != events; i = i + 1)
    {
      read (*,*) fields
      if (fields(1) == 1)
        {
          # Adjust for geometry.
          tmp = fields(4); fields(4) = fields(5); fields(5) = tmp
          tmp = fields(6); fields(6) = fields(7); fields(7) = tmp
        }
      if (fields(2) == 0 && fields(3) == 0)
        {
          N00 = N00 + 1
          N00L1 = N00L1 + fields(4)
          N00L2 = N00L2 + fields(5)
          N00R1 = N00R1 + fields(6)
          N00R2 = N00R2 + fields(7)
        }
      else if (fields(2) == 0 && fields(3) == 1)
        {
          N01 = N01 + 1
          N01L1 = N01L1 + fields(4)
          N01L2 = N01L2 + fields(5)
          N01R1 = N01R1 + fields(6)
          N01R2 = N01R2 + fields(7)
        }
      else if (fields(2) == 1 && fields(3) == 0)
        {
          N10 = N10 + 1
          N10L1 = N10L1 + fields(4)
          N10L2 = N10L2 + fields(5)
          N10R1 = N10R1 + fields(6)
          N10R2 = N10R2 + fields(7)
        }
      else
        {
          N11 = N11 + 1
          N11L1 = N11L1 + fields(4)
          N11L2 = N11L2 + fields(5)
          N11R1 = N11R1 + fields(6)
          N11R2 = N11R2 + fields(7)
        }
    }
end

function kappa (N, NL1, NL2, NR1, NR2)
  implicit none
  integer N, NL1, NL2, NR1, NR2
  double precision kappa
  double precision denom
  double precision sinL, cosL, sinR, cosR
  double precision sinLR, cosLR
  continue
  denom = N
  # Equation (2.4) in the reference.
  sinL = dsqrt (NL1 / denom)
  cosL = dsqrt (NL2 / denom)
  sinR = dsqrt (NR1 / denom)
  cosR = dsqrt (NR2 / denom)
  # Equation (2.3).
  cosLR = (cosR * cosL) + (sinR * sinL)
  sinLR = (sinR * cosL) - (cosR * sinL)
  # Equation (2.5).
  kappa = (cosLR * cosLR) - (sinLR * sinLR)
end

program OpSimA
  implicit none
  integer N00, N01, N10, N11
  integer N00L1, N00L2, N00R1, N00R2
  integer N01L1, N01L2, N01R1, N01R2
  integer N10L1, N10L2, N10R1, N10R2
  integer N11L1, N11L2, N11R1, N11R2
  double precision kappa
  double precision kappa00, kappa01, kappa10, kappa11
  double precision chsh
  continue
  call count (N00, N01, N10, N11,
              N00L1, N00L2, N00R1, N00R2,
              N01L1, N01L2, N01R1, N01R2,
              N10L1, N10L2, N10R1, N10R2,
              N11L1, N11L2, N11R1, N11R2)
  kappa00 = kappa (N00, N00L1, N00L2, N00R1, N00R2)
  kappa01 = kappa (N01, N01L1, N01L2, N01R1, N01R2)
  kappa10 = kappa (N10, N10L1, N10L2, N10R1, N10R2)
  kappa11 = kappa (N11, N11L1, N11L2, N11R1, N11R2)
  chsh = -kappa00 + kappa01 + kappa10 + kappa11
  write (*,'()')
  write (*,'("   light pulse events   ", I9)') n00 + n01 + n10 + n11
  write (*,'()')
  write (*,'("    correlation coefs")')
  write (*,'("          ", F4.1, "° ", F4.1, "°   ", SPF9.6)') _
    angleL1, angleR1, kappa00
  write (*,'("          ", F4.1, "° ", F4.1, "°   ", SPF9.6)') _
    angleL1, angleR2, kappa01
  write (*,'("          ", F4.1, "° ", F4.1, "°   ", SPF9.6)') _
    angleL2, angleR1, kappa10
  write (*,'("          ", F4.1, "° ", F4.1, "°   ", SPF9.6)') _
    angleL2, angleR2, kappa11
  write (*,'()')
  write (*,'("        CHSH contrast   ", SPF9.6)') chsh
  write (*,'("  2*sqrt(2) = nominal   ", SPF9.6)') 2 * dsqrt (2.0D0)
  write (*,'("           difference   ", SPF9.6)') _
    chsh - (2 * dsqrt (2.0D0))
  write (*,'()')
  write (*,'("       CHSH violation   ", SPF9.6)') chsh - 2
  write (*,'()')
end
</syntaxhighlight>

{{out}}
Outputs may look different, due to implementation-dependent behavior of formatted output.

Here is output after compiling with gfortran:
<pre>

   light pulse events      100000

    correlation coefs
           0.0° 22.5°   -0.707806
           0.0° 67.5°   +0.705415
          45.0° 22.5°   +0.712377
          45.0° 67.5°   +0.718882

        CHSH contrast   +2.844480
  2*sqrt(2) = nominal   +2.828427
           difference   +0.016053

       CHSH violation   +0.844480

</pre>
But here is output after compiling with f2c:
<pre>

   light pulse events      100000

    correlation coefs
            .0° 22.5°    -.707806
            .0° 67.5°    +.705415
          45.0° 22.5°    +.712377
          45.0° 67.5°    +.718882

        CHSH contrast   +2.844480
  2*sqrt(2) = nominal   +2.828427
           difference    +.016053

       CHSH violation    +.844480

</pre>
Getting f2c formatted output to print what I would prefer to see is not worth my effort. Probably, if one is using f2c, it is better simply to remove the "SP" that add + signs ahead of numbers. Or, frankly, to hack libf2c more to one's liking. Maybe it was done this way simply to save printer ribbon ink.

=={{header|Wren}}==
{{trans|Python}}
{{libheader|Wren-seq}}
{{libheader|Wren-assert}}
{{libheader|Wren-fmt}}
Unlike the Python example, there is no option to read data from standard input.
<syntaxhighlight lang="wren">import "io" for File
import "os" for Process
import "./seq" for Lst
import "./assert" for Assert
import "./fmt" for Fmt

/* The data for one light pulse. */
class PulseData {
    construct new( logS, logL, logR, detectedL1, detectedL2, detectedR1, detectedR2) {
        _logS = logS
        _logL = logL
        _logR = logR
        _detectedL1 = detectedL1
        _detectedL2 = detectedL2
        _detectedR1 = detectedR1
        _detectedR2 = detectedR2
    }

    logS { _logS }
    logL { _logL }
    logR { _logR }

    detectedL1 { _detectedL1 }
    detectedL2 { _detectedL2 }
    detectedR1 { _detectedR1 }
    detectedR2 { _detectedR2 }

    toString {
        return "PulseData(" + "%(_logS), %(_logL), %(_logR)" +
               "%(_detectedL1), %(_detectedL2), %(_detectedR1), %(_detectedR2))"
    }

    // Swap detectedL1 and detectedL2.
    swapLchannels() {
        var t = _detectedL1
        _detectedL1 = _detectedL2
        _detectedL2 = t
    }

    // Swap detectedR1 and detectedR2.
    swapRchannels() {
        var t = _detectedR1
        _detectedR1 = _detectedR2
        _detectedR2 = t
    }

    // Swap channels on both right and left. This is done if the
    // light source was (1,90°) instead of (1,0°). For in that case
    // the orientations of the polarizing beam splitters, relative to
    // the light source, is different by 90°.
    swapLRchannels() {
        swapLchannels()
        swapRchannels()
    }
}

// Split the data into two subsets, according to whether a set
// item satisfies a predicate. The return value is a tuple, with
// those items satisfying the predicate in the first tuple entry, the
// other items in the second entry.
var splitData = Fn.new { |predicate, data| Lst.partitions(data, predicate) }

// Some data items are for a (1,0°) light pulse. The others are
// for a (1,90°) light pulse. Thus the light pulses are oriented
// differently with respect to the polarizing beam splitters. We
// adjust for that distinction here.
var adjustDataForLightPulseOrientation = Fn.new { |data|
    var pf = Fn.new { |item| item.logS == 0 }
    var split = splitData.call(pf, data)
    var data0  = split[0]
    var data90 = split[1]
    for (item in data90) item.swapLRchannels()
    return data0 + data90
}

// Split the data into four subsets: one subset for each
// arrangement of the two polarizing beam splitters. 
var splitDataAccordingToPBSSetting = Fn.new { |data|
    var pfL = Fn.new { |item| item.logL == 0 }
    var split = splitData.call(pfL, data)
    var dataL1 = split[0]
    var dataL2 = split[1]
    var pfR = Fn.new { |item| item.logR == 0 }
    var splitL1 = splitData.call(pfR, dataL1)
    var dataL1R1 = splitL1[0]
    var dataL1R2 = splitL1[1]
    var splitL2 = splitData.call(pfR, dataL2)
    var dataL2R1 = splitL2[0]
    var dataL2R2 = splitL2[1]
    return [dataL1R1, dataL1R2, dataL2R1, dataL2R2]
}

// Compute the correlation coefficient for the subset of the data
// that corresponding to a particular setting of the polarizing beam splitters.
var computeCorrelationCoefficient = Fn.new { |angleL, angleR, data|
    // We make certain the orientations of beam splitters are
    // represented by perpendicular angles in the first and fourth
    // quadrant. This restriction causes no loss of generality, because
    // the orientation of the beam splitter is actually a rotated "X".
    Assert.ok([angleL, angleR].all { |x| 0 <= x && x < 90 })
    var perpendicularL = angleL - 90 // In Quadrant 4.
    var perpendicularR = angleR - 90 // In Quadrant 4.

    // Note that the sine is non-negative in Quadrant 1, and the cosine
    // is non-negative in Quadrant 4. Thus we can use the following
    // estimates for cosine and sine. This is Equation (2.4) in the
    // reference. (Note, one can also use Quadrants 1 and 2 and reverse
    // the roles of cosine and sine. And so on like that.)
    var N = data.count
    var NL1 = 0
    var NL2 = 0
    var NR1 = 0
    var NR2 = 0
    for (item in data) {
        NL1 = NL1 + item.detectedL1
        NL2 = NL2 + item.detectedL2
        NR1 = NR1 + item.detectedR1
        NR2 = NR2 + item.detectedR2
    }
    var sinL = (NL1 / N).sqrt
    var cosL = (NL2 / N).sqrt
    var sinR = (NR1 / N).sqrt
    var cosR = (NR2 / N).sqrt

    // Now we can apply the reference's Equation (2.3).
    var cosLR = (cosR * cosL) + (sinR * sinL)
    var sinLR = (sinR * cosL) - (cosR * sinL)

    // And then Equation (2.5).
    return (cosLR * cosLR) - (sinLR * sinLR)
}

// Read the raw data into a list.
var readRawData = Fn.new { |filename|
    var makeRecord = Fn.new { |line|
        var x = line.split(" ")
        x = x.map { |s| Num.fromString(s).truncate }.toList
        return PulseData.new(x[0], x[1], x[2], x[3], x[4], x[5], x[6])
    }

    var data = []
    var lines = File.read(filename).split("\n")
    var numPulses = Num.fromString(lines[0])
    for (i in 1..numPulses) {
        data.add(makeRecord.call(lines[i]))
    }
    return data
}

var args = Process.arguments
if (args.count != 1) {
    System.print("Please provide 1 command line argument: RAW_DATA_FILENAME")
    return
}
var filename = args[0]

// Polarizing beam splitter orientations commonly used in actual
// experiments. These must match the values used in the simulation
// itself. They are by design all either zero degrees or in the
// first quadrant.
var anglesL = [0, 45]
var anglesR = [22.5, 67.5]
Assert.ok((anglesL + anglesR).all { |x| 0 <= x && x < 90 })

var data = readRawData.call(filename)
data = adjustDataForLightPulseOrientation.call(data)
var split = splitDataAccordingToPBSSetting.call(data)
var dataL1R1 = split[0]
var dataL1R2 = split[1]
var dataL2R1 = split[2]
var dataL2R2 = split[3]

var kappaL1R1 = computeCorrelationCoefficient.call(anglesL[0], anglesR[0], dataL1R1)
var kappaL1R2 = computeCorrelationCoefficient.call(anglesL[0], anglesR[1], dataL1R2)
var kappaL2R1 = computeCorrelationCoefficient.call(anglesL[1], anglesR[0], dataL2R1)
var kappaL2R2 = computeCorrelationCoefficient.call(anglesL[1], anglesR[1], dataL2R2)

var chshContrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2

// The nominal value of the CHSH contrast for the chosen polarizer
// orientations is 2*sqrt(2).
var chshContrastNominal = 2.sqrt * 2

Fmt.print()
Fmt.print("   light pulse events      $,10d", data.count)
Fmt.print()
Fmt.print("    correlation coefs")
Fmt.print("              $4.1f° $4.1f°   $+9.6f", anglesL[0], anglesR[0], kappaL1R1)
Fmt.print("              $4.1f° $4.1f°   $+9.6f", anglesL[0], anglesR[1], kappaL1R2)
Fmt.print("              $4.1f° $4.1f°   $+9.6f", anglesL[1], anglesR[0], kappaL2R1)
Fmt.print("              $4.1f° $4.1f°   $+9.6f", anglesL[1], anglesR[1], kappaL2R2)
Fmt.print()
Fmt.print("            CHSH contrast   $+9.6f", chshContrast)
Fmt.print("      2*sqrt(2) = nominal   $+9.6f", chshContrastNominal)
Fmt.print("               difference   $+9.6f", chshContrast - chshContrastNominal)

// A "CHSH violation" occurs if the CHSH contrast is > 2.
// https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
Fmt.print()
Fmt.print("           CHSH violation   $+9.6f", chshContrast - 2)
Fmt.print()</syntaxhighlight>

{{out}}
An example using data from the Wren implementation of ''[[Simulated optics experiment/Simulator]]''.
<pre>

   light pulse events         100,000

    correlation coefs
               0.0° 22.5°   -0.702281
               0.0° 67.5°   +0.702892
              45.0° 22.5°   +0.696311
              45.0° 67.5°   +0.710296

            CHSH contrast   +2.811781
      2*sqrt(2) = nominal   +2.828427
               difference   -0.016646

           CHSH violation   +0.811781
</pre>