View source for Simulated optics experiment/Simulator

{{draft task}}

{{alertbox|#DDDAE8|There is now an "extra credit" problem. One might wish to have a look at it.}}

===Introduction===

Simulation by computer modeling is done in the natural sciences for various reasons. Consider, for instance, simulation of the weather: it might be done to try to predict the actual weather for the next few days.

In this task, you will write a simulation of an experiment in optics, the physics of light. The simulation will be based on the following [[wp:Open_access|open access]] paper, but written independently of the paper's author's own simulation:
    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 purpose of the simulation is demonstrate flaws in some physicists' interpretations of certain actual experiments, by presenting a simulated counterexample to their reasoning.

===Task description===

In this task you should write the program or set of programs that simulate the experimental apparatus, so generating raw data. The ''[[Simulated optics experiment/Data analysis]]'' task handles analysis of the raw data. There should be no need for the simulator and the data analysis to be in the same programming language. Therefore any implementation of the data analysis can be used to check your implementation of the simulator.

Write simulations of the following experimental devices. The descriptions may seem complicated, but the [[#ObjectIcon|Object Icon]], [[#Python|Python]], and [[#RATFOR|RATFOR]] examples can serve as references. The RATFOR is certainly the simplest, but the Object Icon may seem more familiar to "object oriented" programmers. The two are similar. The Python reformulates some calculations as Gibbs vectors (the "arrow" kind of vector), and runs as multiple processes.

* A ''light source'', which occasionally emits two pulses of [[wp:Plane_of_polarization|plane-polarized light]], one towards the left, the other towards the right. Both pulses have an amplitude of 1. About half the time the pulses are polarized, respectively, at an angle of 0° on the left and 90° on the right. The rest of the time the angles are reversed: 90° on the left, 0° on the right. A random number generator should be used to select between the two settings. If the 0° angle is on the left, then a "0" should be recorded in a log. Otherwise a "1" is recorded.

* A [[wp:Polarizer|''polarizing beam splitter'']]. Actually you will need two of these devices, but the two should share their implementation. Each beam splitter has an angle setting, which for our own convenience is required to be greater than or equal to 0° but less than 90°. (The simulation in the paper does not impose this requirement, but imposing it does not change the results.) A polarizing beam splitter does the following. For input it receives one of the light pulses emitted by the light source, and for output it emits two new light pulses, whose directions of travel we will not worry about. One of the new light pulses has amplitude equal to the cosine of the difference between the angle of the incoming light and the angular setting of the beam splitter. The other new light pulse has amplitude equal to the absolute value of the sine of the difference between angles. The new light pulses are plane-polarized but this information is not used by the next device in line, so may be ignored if one wishes. (The reference [[#Python|Python]] example ''does'' compute the directions of polarization.) The beam splitters will actually have two angle settings, with the first of the two settings chosen randomly about half the time. If the first angle setting was chosen, "0" is recorded in a log. Otherwise a "1" is recorded.

* A [[wp:Photodetector|''light detector'']], or actually four of them that share their implementation. Each light detector receives as input, respectively, one of the four output pulses from the two polarizing beam splitters. A light detector first squares the amplitude of the incoming light pulse. This square is the [[wp:Intensity_(physics)|''intensity'']] of the pulse. Then it compares the intensity it just computed with a uniform random number between 0.0 and 1.0. If the random number is less than or equal to the intensity, the light detector outputs a "1", meaning that it has detected a light pulse. Otherwise the detector outputs a "0", representing the quiescent state of the detector. (That is, the detector has failed to detect the pulse.) The output of each light detector is recorded in a log.

The angle settings of the polarizing beam splitters will be as follows:
* On the left, the first setting angle is 0° and the second is 45°.
* On the right, the first angle is 22.5° and the second is 67.5°.

The simulation is run by having the light source emit some number of pulses and letting the other devices do their work. How you arrange this to work is up to you. The [[#Python|Python]] example, for instance, runs each device as a separate process, connected to each other by queues. But you can instead have, for instance, an event loop, coroutines, etc.--or even just ordinary, procedural calculation of the numbers. The last method is simplest, and perfectly correct. It is what the [[#ObjectIcon|Object Icon]] and [[#RATFOR|RATFOR]] examples do. However, surely all the methods have their place in the world of simulations.

The program must output its "raw data" results in the format shown here by example:
<pre>3
0 1 1 0 0 1 1
1 1 0 1 1 0 1
0 0 1 0 0 1 1</pre>
The first line is how many pulses were emitted by the light source. (You should have this be at least 1000. The number 3 here is for the sake of depicting the format.) This line is followed by that many more lines, each of which contains seven "0"s and "1"s, separated from each other by one space character. The seven entries, left to right, represent the following data, one line for each light pulse pair emitted by the source:
# The log recording of the light source setting.
# The log recording of the setting of the polarizing beam splitter that is on the left.
# The log recording of the setting of the polarizing beam splitter that is on the right.
# The output of the light detector on the left that is receiving the "cosine" pulses.
# The output of the light detector on the left that is receiving the "sine" pulses.
# The output of the light detector on the right that is receiving the "cosine" pulses.
# The output of the light detector on the right that is receiving the "sine" pulses.

You should feed this "raw data" output as input to one of the ''[[Simulated optics experiment/Data analysis]]'' programs and display your results.

===Extra "credit"===
The "credit" one gets is the feeling of having contributed to the community.

The simulation above uses classical optics theory, including the [https://en.wikipedia.org/w/index.php?title=Polarizer&oldid=1152014034 Law of Malus], and assumes light detectors are the source of "randomness" in detections. Instead write a simulation with these differences:

* The ''light source'' works exactly as above, but now we call it ''a source of photons-containing-hidden-variables'' and pay no attention to the amplitude of the light pulses. They are now simply little "pellets" of light, but containing some inner state that quantum mechanics pointedly ignores.
* A polarizing beam splitter, rather than emit light of reduced amplitude, emits up to two new photons-containing-hidden-variables, again of arbitrary amplitude. A photon-containing-hidden-variables possibly is emitted towards one of the light detectors, with probability equal to the square of the cosine of the difference in angle between the impinging photon and the beam splitter. The other light detector gets a photon-containing-hidden-variables with probability equal to the square of the sine.
* The ''light detector'' we will now call a ''photodetector''. It detects impinging photons-containing-hidden-variables with probability ''one''. It is a perfect photon-containing-hidden-variables detector.
* Output must be in the format described above, so the data analyzers can analyze them.

Though the simulation now includes "photons" instead of classical optics, it will be perfectly valid for the data analyzer to take into account which direction of travel got which photon. This is because these are not "photons" as physicists usually assume them to be, but ''photons-containing-hidden-variables''. The analysis therefore ''must'' take into account which photon is which.

All but a few physicists either claim photons-containing-hidden-variables cannot possibly explain actual experimental results, or take the word of other physicists for this. There supposedly is a "mathematical" proof of it. Our simulation should end up contradicting that claim, indicating that the "mathematics" ''must be'' incorrect in some fashion. How to interpret that result that I leave to the individual.

(It is not necessary to have an experimentally proven model of photons-containing-hidden-variables, for the simulation to achieve its goal. Such a model probably would explain deterministically why some photons get "re-emitted" and some do not, whereas we are assuming that polarizing beam splitters respond in some unknown fashion. But, in the realm of ''mathematics'', even a ''facetious'' counterexample would be absolute disproof of the "theorem".)

There is [[#Extra_"credit"_simulator_in_ATS|an example]] of the extra "credit" simulation in ML-like [[#ATS|ATS]].

===See also===
* ''[[Simulated optics experiment/Data analysis]]''

=={{header|ATS}}==
===Extra "credit" simulator in ATS===

The simulation is a synchronous event loop, with a small probability of photon-with-hidden-variables emission each time through the loop. The code is written in an ML-like style (with "ref" for the mutable variables, for instance). There is a fair bit of C code embedded in the source, including a good random number generator copied from Wikipedia.

<syntaxhighlight lang="ats">
#include "share/atspre_staload.hats"

(* Simulation is by synchronous event loop run for a given amount of
   clock time, then allowed to settle. *)

#define ATS_EXTERN_PREFIX ""

#define NIL list_nil ()
#define ::  list_cons

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

typedef ident = intGte 0  (* event identification, for syncing data *)
typedef zero_or_one = intBtwe (0, 1)
typedef log_entry = @(ident, zero_or_one)
typedef logbook = List0 log_entry

datatype pwhv (tk : tkind) =   (* "photon with hidden variables" *)
| pwhv of (ident, g0float tk)  (* angle of polarization, in degrees *)

datatype pwhv_channel (tk : tkind) = (* a queue *)
| pwhv_channel of
    ref (List0 (pwhv tk))

datatype pwhv_source (tk : tkind) = (* source of pwhv pairs *)
| pwhv_source of
    (g0float tk,                (* probability of emission *)
     ref ident,                 (* the next identification number *)
     ref logbook,               (* log of hidden variables *)
     pwhv_channel tk,           (* leftwards emission *)
     pwhv_channel tk)           (* rightwards emission *)

datatype pbs (tk : tkind) =     (* "polarizing beam splitter" *)
| pbs of (g0float tk,           (* angle setting 0, in degrees *)
          g0float tk,           (* angle setting 1, in degrees *)
          ref logbook,          (* log of settings *)
          pwhv_channel tk,      (* impinging pwhv *)
          pwhv_channel tk,      (* cos² emission *)
          pwhv_channel tk)      (* sin² emission *)

datatype photodetector (tk : tkind) =
| photodetector of
    (ref logbook,               (* log of detections *)
     pwhv_channel tk)           (* impinging pwhv *)

(*------------------------------------------------------------------*)
(* Some math functionality, for all the standard floating point
   types. The ats2-xprelude package includes this, and more, but one
   may wish to avoid the dependency. And there is support for math
   functions in libats/libc, but not with typekinds. *)

%{^
#include <math.h>
// cospi(3) and sinpi(3) would be better than cos(3) and sin(3), but
// I do not yet have cospi(3) or sinpi(3).
#define pi__ \
  3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214L
#define cospi_float(x) (cosf (((atstype_float) pi__) * (x)))
#define cospi_double(x) (cos (((atstype_double) pi__) * (x)))
#define cospi_ldouble(x) (cosl (((atstype_ldouble) pi__) * (x)))
#define sinpi_float(x) (sinf (((atstype_float) pi__) * (x)))
#define sinpi_double(x) (sin (((atstype_double) pi__) * (x)))
#define sinpi_ldouble(x) (sinl (((atstype_ldouble) pi__) * (x)))
%}

extern fn cospi_float : float -<> float = "mac#%"
extern fn cospi_double : double -<> double = "mac#%"
extern fn cospi_ldouble : ldouble -<> ldouble = "mac#%"
extern fn {tk : tkind} g0float_cospi : g0float tk -<> g0float tk
implement g0float_cospi<fltknd> x = cospi_float x
implement g0float_cospi<dblknd> x = cospi_double x
implement g0float_cospi<ldblknd> x = cospi_ldouble x

extern fn sinpi_float : float -<> float = "mac#%"
extern fn sinpi_double : double -<> double = "mac#%"
extern fn sinpi_ldouble : ldouble -<> ldouble = "mac#%"
extern fn {tk : tkind} g0float_sinpi : g0float tk -<> g0float tk
implement g0float_sinpi<fltknd> x = sinpi_float x
implement g0float_sinpi<dblknd> x = sinpi_double x
implement g0float_sinpi<ldblknd> x = sinpi_ldouble x

overload cospi with g0float_cospi
overload sinpi with g0float_sinpi

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

%{^
// PCG-XSH-RR random number generator, from Wikipedia.
// https://en.wikipedia.org/w/index.php?title=Permuted_congruential_generator&oldid=1136326070#Example_code

#include <stdint.h>
static uint64_t       state      = 0x4d595df4d0f33173;		// Or something seed-dependent
static uint64_t const multiplier = 6364136223846793005u;
static uint64_t const increment  = 1442695040888963407u;	// Or an arbitrary odd constant

static uint32_t rotr32(uint32_t x, unsigned r)
{
	return x >> r | x << (-r & 31);
}

uint32_t pcg32(void)
{
	uint64_t x = state;
	unsigned count = (unsigned)(x >> 59);		// 59 = 64 - 5

	state = x * multiplier + increment;
	x ^= x >> 18;								// 18 = (64 - 27)/2
	return rotr32((uint32_t)(x >> 27), count);	// 27 = 32 - 5
}

void pcg32_init(uint64_t seed)
{
	state = seed + increment;
	(void)pcg32();
}
%}

extern fn pcg32 : () -> uint32 = "mac#%"
extern fn pcg32_init : uint64 -> void = "mac#%"

fn {tk : tkind}
randnum () : g0float tk =
  (* Returns a random floating point number in [0,1) *)
  let
    extern castfn u32_to_ldouble : uint32 -<> ldouble

    (* THE FOLLOWING WILL WORK ONLY IF g0float tk IS OF A TYPE THE C
       COMPILER CAN CAST TO. This is not true of all the types in
       ats2-xprelude. OTOH you wouldn’t need this cast in the first
       place, if you were using ats2-xprelude. *)
    extern castfn ldouble2tk : ldouble -<> g0float tk
  in
    ldouble2tk ((u32_to_ldouble (pcg32 ())) / 4294967296.0L)
  end

%{^
#include <time.h>
#define get_time_as_uint64() ((atstype_uint64) (time (NULL)))
%}

fn {}
randomize () : void =
  (* Seed the random number generator with something variable enough
     to keep the program’s user entertained. *)
  let
    (* Seconds since the Epoch. *)
    extern fn get_time_as_uint64 : () -> uint64 = "mac#%"
  in
    pcg32_init (get_time_as_uint64 ())
  end

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

fn {tk : tkind}
pwhv_channel_put (channel : pwhv_channel tk,
                  pwhv    : pwhv tk) =
  let
    val+ pwhv_channel lst = channel
  in
    (* Inefficient if queues are large. Our queues will be tiny. *)
    !lst := !lst + (pwhv :: NIL)
  end

fn {tk : tkind}
pwhv_channel_get (channel : pwhv_channel tk) : Option (pwhv tk) =
  let
    val+ pwhv_channel lst = channel
  in
    case+ !lst of
    | NIL => None ()
    | hd :: tl => (!lst := tl; Some hd)
  end

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

fn {tk : tkind}
pwhv_source_poll (source : pwhv_source tk) : void =
  let
    val+ pwhv_source (probability_of_emission, ident, hv_log,
                      channel_L, channel_R) = source
  in
    if randnum<tk> () < probability_of_emission then
      let
        val id = !ident
      in
        !ident := succ id;
        if randnum<tk> () < g0f2f 0.5 then
          begin
            !hv_log := @(id, 0) :: !hv_log;
            pwhv_channel_put<tk> (channel_L, pwhv (id, g0i2f 0));
            pwhv_channel_put<tk> (channel_R, pwhv (id, g0i2f 90))
          end
        else
          begin
            !hv_log := @(id, 1) :: !hv_log;
            pwhv_channel_put<tk> (channel_L, pwhv (id, g0i2f 90));
            pwhv_channel_put<tk> (channel_R, pwhv (id, g0i2f 0))
          end
      end
  end

fn {tk : tkind}
pbs_poll (splitter : pbs tk) : void =
  let
    val+ pbs (angle0, angle1, pbs_log,
              channel_in, channel0, channel1) = splitter
  in
    case+ pwhv_channel_get (channel_in) of
    | None () => ()
    | Some (pwhv (id, angle)) =>
      if randnum<tk> () < g0f2f 0.5 then
        let
          val c = cospi ((angle - angle0) / g0i2f 180)
          and s = sinpi ((angle - angle0) / g0i2f 180)
        in
          !pbs_log := @(id, 0) :: !pbs_log;
          if randnum<tk> () < c * c then
            pwhv_channel_put<tk> (channel0, pwhv (id, angle0));
          if randnum<tk> () < s * s then
            pwhv_channel_put<tk> (channel1,
                                  pwhv (id, angle0 + g0i2f 90))
        end
      else
        let
          val c = cospi ((angle - angle1) / g0i2f 180)
          and s = sinpi ((angle - angle1) / g0i2f 180)
        in
          !pbs_log := @(id, 1) :: !pbs_log;
          if randnum<tk> () < c * c then
            pwhv_channel_put<tk> (channel0, pwhv (id, angle1));
          if randnum<tk> () < s * s then
            pwhv_channel_put<tk> (channel1,
                                  pwhv (id, angle1 + g0i2f 90))
        end
  end

fn {tk : tkind}
photodetector_poll (detector : photodetector tk) : void =
  let
    val+ photodetector (pd_log, channel_in) = detector
  in
    case+ pwhv_channel_get (channel_in) of
    | None () => ()
    | Some (pwhv (id, angle)) => (!pd_log := @(id, 1) :: !pd_log)
  end

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

fn {tk : tkind}
event_loop (num_events  : intGte 0,
            source      : pwhv_source tk,
            splitter_L  : pbs tk,
            splitter_R  : pbs tk,
            detector_Lc : photodetector tk,
            detector_Ls : photodetector tk,
            detector_Rc : photodetector tk,
            detector_Rs : photodetector tk) : void =
  let
    val+ pwhv_source (_, ident, _, _, _) = source
    fun
    main_loop () : void =
      if !ident < num_events then
        begin
          pwhv_source_poll (source);
          pbs_poll (splitter_L);
          pbs_poll (splitter_R);
          photodetector_poll (detector_Lc);
          photodetector_poll (detector_Ls);
          photodetector_poll (detector_Rc);
          photodetector_poll (detector_Rs);
          main_loop ()
        end
    fun
    finishing_loop (n : intGte 0) : void =
      if n <> 0 then
        begin
          pbs_poll (splitter_L);
          pbs_poll (splitter_R);
          photodetector_poll (detector_Lc);
          photodetector_poll (detector_Ls);
          photodetector_poll (detector_Rc);
          photodetector_poll (detector_Rs);
          finishing_loop (pred n)
        end
  in
    main_loop ();

    (* Actually some very small number is all that is needed here. But
       this bigger number does no harm. *)
    finishing_loop (100)
  end

fn {tk : tkind}
run_simulation (outf       : FILEref,
                angle_L0   : g0float tk,
                angle_L1   : g0float tk,
                angle_R0   : g0float tk,
                angle_R1   : g0float tk,
                probability_of_emission : g0float tk,
                num_events : intGte 0) : void =
  let
    val pbs2detector_Lc = pwhv_channel (ref NIL)
    and pbs2detector_Ls = pwhv_channel (ref NIL)
    and pbs2detector_Rc = pwhv_channel (ref NIL)
    and pbs2detector_Rs = pwhv_channel (ref NIL)
    and source2pbs_L = pwhv_channel (ref NIL)
    and source2pbs_R = pwhv_channel (ref NIL)

    val detector_Lc = photodetector (ref NIL, pbs2detector_Lc)
    and detector_Ls = photodetector (ref NIL, pbs2detector_Ls)
    and detector_Rc = photodetector (ref NIL, pbs2detector_Rc)
    and detector_Rs = photodetector (ref NIL, pbs2detector_Rs)
    and splitter_L = pbs (angle_L0, angle_L1, ref NIL,
                          source2pbs_L, pbs2detector_Lc,
                          pbs2detector_Ls)
    and splitter_R = pbs (angle_R0, angle_R1, ref NIL,
                          source2pbs_R, pbs2detector_Rc,
                          pbs2detector_Rs)
    and source = pwhv_source (probability_of_emission, ref 0,
                              ref NIL, source2pbs_L, source2pbs_R)

    val () = event_loop<tk> (num_events, source,
                             splitter_L, splitter_R,
                             detector_Lc, detector_Ls,
                             detector_Rc, detector_Rs);

    val+ pwhv_source (_, _, hv_log, _, _) = source
    val+ pbs (_, _, pbs_log_L, _, _, _) = splitter_L
    val+ pbs (_, _, pbs_log_R, _, _, _) = splitter_R
    val+ photodetector (pd_log_Lc, _) = detector_Lc
    val+ photodetector (pd_log_Ls, _) = detector_Ls
    val+ photodetector (pd_log_Rc, _) = detector_Rc
    val+ photodetector (pd_log_Rs, _) = detector_Rs

    val hv_log = !hv_log
    and pbs_log_L = !pbs_log_L
    and pbs_log_R = !pbs_log_R
    and pd_log_Lc = !pd_log_Lc
    and pd_log_Ls = !pd_log_Ls
    and pd_log_Rc = !pd_log_Rc
    and pd_log_Rs = !pd_log_Rs

    val n = length hv_log
    and n_L = length pbs_log_L
    and n_R = length pbs_log_R
    and n_Lc = length pd_log_Lc
    and n_Ls = length pd_log_Ls
    and n_Rc = length pd_log_Rc
    and n_Rs = length pd_log_Rs

    prval [n : int] EQINT () = eqint_make_gint n

    val () = assertloc (n_L = n)
    val () = assertloc (n_R = n)
    val () = assertloc (n_Lc <= n)
    val () = assertloc (n_Ls <= n)
    val () = assertloc (n_Rc <= n)
    val () = assertloc (n_Rs <= n)

    val data = mtrxszref_make_elt<zero_or_one> (i2sz n, i2sz 7, 0)
    fun
    transfer_data_from_logbook (log : logbook,
                                j   : intBtwe (0, 6)) : void =
      case+ log of
      | NIL => ()
      | @(id, entry) :: tl =>
        begin
          data[id, j] := entry;
          transfer_data_from_logbook (tl, j)
        end

    val () = transfer_data_from_logbook (hv_log, 0)
    val () = transfer_data_from_logbook (pbs_log_L, 1)
    val () = transfer_data_from_logbook (pbs_log_R, 2)
    val () = transfer_data_from_logbook (pd_log_Lc, 3)
    val () = transfer_data_from_logbook (pd_log_Ls, 4)
    val () = transfer_data_from_logbook (pd_log_Rc, 5)
    val () = transfer_data_from_logbook (pd_log_Rs, 6)

    fun
    fprint_data (outf : FILEref,
                 i    : intBtwe (0, n)) : void =
      if i <> n then
        begin
          fprintln! (outf, data[i, 0], " ", data[i, 1], " ",
                     data[i, 2], " ", data[i, 3], " ",
                     data[i, 4], " ", data[i, 5], " ", data[i, 6]);
          fprint_data (outf, succ i)
        end
  in
    fprintln! (outf, n);
    fprint_data (outf, 0)
  end

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

implement
main0 (argc, argv) =
  let
    val args = list_vt2t (listize_argc_argv (argc, argv))

    val num_events =
      if argc <= 1 then
        1000
      else
        (max ($extfcall (Int, "atoi", args[1]), 0)) : intGte 0

    val outf = stdout_ref

    val () = randomize ()

    val angle_L0 = 0.0
    and angle_L1 = 45.0
    and angle_R0 = 22.5
    and angle_R1 = 67.5
    and probability_of_emission = 0.05
  in
    run_simulation (outf, angle_L0, angle_L1, angle_R0, angle_R1,
                    probability_of_emission, num_events)
  end
</syntaxhighlight>

{{out}}
Though this is a low-efficiency ATS program, it still is immensely more quick than Python or Object Icon. So I did a one-million event simulation. As you can see, even if you model the light as little "pellets", you get an absolutely humongous CHSH violation of about 0.83. (In the context, 0.83 is enormous. :) )

The only thing necessary to get this result is to remember that you are trying to model photons ''with'' hidden variables, so in the analysis you actually ''are not allowed'' to leave out relevant information you have about the creation of the photons. It’s presumed to be part of the state of the physical object.

<pre>

   light pulse events     1000000

    correlation coefs
           0.0° 22.5°   -0.707575
           0.0° 67.5°   +0.707097
          45.0° 22.5°   +0.707377
          45.0° 67.5°   +0.707951

        CHSH contrast   +2.830000
  2*sqrt(2) = nominal   +2.828427
           difference   +0.001573

       CHSH violation   +0.830000

</pre>

=={{header|Julia}}==
Does direct calculation of cos and abs(sin) rather than the vector functions the Python example uses. This causes small differences in
output from differences in rounding error when the angle difference for the polarizing splitter is close to 0, but does not change the overall results.

<small>Actually one expects differences, because the simulation uses an unspecified random number generator. There is actually no ''good'' reason to use a random number generator rather than, say, evenly spaced numbers, but when trying to make a point one might pretend "randomness" matters. --[[User:Chemoelectric|Chemoelectric]] ([[User talk:Chemoelectric|talk]]) 18:42, 30 May 2023 (UTC)</small>

<syntaxhighlight lang=Julia>""" rosettacode.org/wiki/Simulated_optics_experiment/Simulator """

using Pipe

"""
    Simulate an output which is randomly either [0, 90] or [90, 0].
    Log to the first integer in outputlogline.
"""
function pairedlightsource!(outputlogline)
    zeroleft = rand([false, true])
    outputlogline[begin] = zeroleft
    return zeroleft ? [0, 90] : [90, 0]
end

"""
    Given either [0, 90] or [90, 0] output a 2 X 2 matrix of outputs
    using degrees for sind and cosd functions.
    Logs output of the random choices to 2nd and 3rd integers in outputlogline.
"""
function beamsplitters!(outputlogline, lightsources, settings)
    output, choices = zeros(length(lightsources), 2), rand([0, 1], length(lightsources))
    for i in eachindex(lightsources)
        outputlogline[begin + i] = choices[i]
        ang = settings[i][choices[i] + 1]
        output[i, :] .= (cosd(ang - lightsources[i]), abs(sind(ang - lightsources[i])))
    end
    return output
end

"""
    Given a 2 X 2 matrix of reals between 0 and 1, output a vector of 4 Bool values
    depending whether the squares of each value is <= or > an individually generated
    uniform random value between 0 and 1.
    Logs output of the random choices to 4th through 7th integers in outputlogline.
"""
function lightdetectors!(outputlogline, splitsources)
    for (i, src) in enumerate(splitsources)
        outputlogline[begin + 2 + i] = src * src <= rand()
    end
    return outputlogline
end

"""
    Run the simulation. Output is in the format as quoted from task:

    (quote)
    The program must output its "raw data" results in the format shown here by example:

        3
        0 1 1 0 0 1 1
        1 1 0 1 1 0 1
        0 0 1 0 0 1 1
        The first line is how many pulses were emitted by the light source. (You should have
        this be at least 1000. The number 3 here is for the sake of depicting the format.)
        This line is followed by that many more lines, each of which contains seven "0"s and "1"s,
        separated from each other by one space character. The seven entries, left to right,
        represent the following data, one line for each light pulse pair emitted by the source:

        The log recording of the light source setting.
        The log recording of the setting of the polarizing beam splitter that is on the left.
        The log recording of the setting of the polarizing beam splitter that is on the right.
        The output of the light detector on the left that is receiving the "cosine" pulses.
        The output of the light detector on the left that is receiving the "sine" pulses.
        The output of the light detector on the right that is receiving the "cosine" pulses.
        The output of the light detector on the right that is receiving the "sine" pulses.
    (end quote)

    If `datalogfile` is provided as a named argment, write output lines to that file.
"""
function lightsimulator(npulses; settings = [[0, 45], [22.5, 67.5]], datalogfile = "")
    simulatorloglines = zeros(Bool, npulses, 7) # log line entries must be either 0 or 1

    @Threads.threads for a in eachrow(simulatorloglines)
        @pipe pairedlightsource!(a) |> beamsplitters!(a, _, settings) |> lightdetectors!(a, _)
    end

    println(npulses)
    if npulses < 50 # don't print them all if this is a large number of data points
        foreach(a -> println(join(Int.(a), " ")), eachrow(simulatorloglines))
    end

    # if we need to save output to a data file
    if datalogfile != ""
        fh = open(datalogfile, write = true)
        println(fh, npulses)
        for a in eachrow(simulatorloglines)
            println(fh, join(string.(Int.(a)), " "))
        end
        close(fh)
    end
end

lightsimulator(1_000_000, datalogfile = "datalog.log")
</syntaxhighlight>{{out}}
<pre>
Sample output from Python or Julia data analysis examples:
1000000 

        
   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}}
This is a reasonably faithful translation of Python code with however some important differences.

We use threads instead of processes and communication is done via channels instead of queues. Note that as Nim
is a statically typed language, channels must specify the type of data exchanged. Fortunately, in the Python code each
queue is used with a single type: either integer, vector or a dictionary. So we were able to specify the type of data
to exchange between threads.

As there is no way to kill a thread, we added a function to stop a mechanism (and the associated thread). We could also
have used a global termination boolean but using a function is cleaner.

Due to the way channels work, we had to use pointers to transmit them between threads.

We kept the Mechanism class (defined as an inheritable object type in Nim). We added a termination boolean and a state
for the PRNG. Indeed, the default random generator is not thread safe, so we had to manage one for each thread (i.e.
each instance of Mechanism associated to a thread).

<syntaxhighlight lang="Nim">import std/[math, os, random, sequtils, strutils, tables]


####################################################################################################
# Vector.

# A simple implementation of Gibbs vectors, suited to our purpose.
# A vector is stored in polar form, with the angle in
# degrees between 0 (inclusive) and 360 (exclusive).
type Vector = object
  magnitude: float
  angle: float

func initVector*(magnitude, angle: float): Vector =
  Vector(magnitude: magnitude, angle: floorMod(angle, 360))

func `$`*(v: Vector): string =
  "Vector($1, $2)".format(v.magnitude, v.angle)

func fromRect*(x, y: float): Vector =
  ## Return a vector for the given rectangular coordinates.
  Vector(magnitude: hypot(x, y), angle: floorMod(arctan2(y, x).radToDeg, 360))

func toRect*(v: Vector): (float, float) =
  ## Return the x and y coordinates of the vector.
  result[0] = v.magnitude * cos(v.angle.degToRad)
  result[1] = v.magnitude * sin(v.angle.degToRad)

func `*`*(v: Vector; scalar: float): Vector =
  ## Multiply a vector by a scalar, returning a new vector.
  Vector(magnitude: v.magnitude * scalar, angle: v.angle)

func dotProduct*(u, v: Vector): float =
  ## Return the dot product of two vectors.
  u.magnitude * v.magnitude * cos(degToRad(u.angle - v.angle))

func `-`*(u, v: Vector): Vector =
  ## Return the difference of two vectors.
  let (xu, yu) = u.toRect()
  let (xv, yv) = v.toRect()
  result = fromRect(xu - xv, yu - yv)

func projection*(u, v: Vector): Vector =
  ## Return the projection of vector "u" onto vector "v".
  result = v * (dotProduct(u, v) / dotProduct(v, v))


####################################################################################################
# Channels.

type

  DataOut = Table[string, int]    # Representation of output data.

  IntChannel = ptr Channel[int]
  VectorChannel = ptr Channel[Vector]
  DataChannel = ptr Channel[DataOut]

proc newIntChannel(): IntChannel =
  cast[IntChannel](allocShared0(sizeof(Channel[int])))

proc newVectorChannel(): VectorChannel =
  cast[VectorChannel](allocShared0(sizeof(Channel[Vector])))

proc newDataChannel(): DataChannel =
  cast[DataChannel](allocShared0(sizeof(Channel[DataOut])))


####################################################################################################
# Mechanism.

# A Mechanism represents a part of the experimental apparatus.
type Mechanism = ref object of RootObj
  randState: Rand   # Random generator state (one per instance, so one per thread).
  operating: bool   # True if the mechanism is operating.

proc init(m: Mechanism) =
  ## Initialize a mechanism.
  m.randState = initRand()

method run(m: Mechanism) {.base, gcsafe.}=
  ## Virtual method to override.
  raise newException(CatchableError, "method without implementation override.")

proc start(m: Mechanism) {.thread.} =
  ## Run the mechanism.
  m.operating = true
  while m.operating:
    m.run()
    # A small pause of one ms to try not to overtax the computer.
    sleep(1)

func stop(m: Mechanism) =
  ## Stop the mechanism.
  m.operating = false


####################################################################################################
# LightSource.

# A LightSource occasionally emits oppositely plane-polarized
# light pulses, of fixed amplitude, polarized 90° with respect to
# each other. The pulses are represented by the vectors (1,0°) and
# (1,90°), respectively. One is emitted to the left, the other to
# the right. The probability is 1/2 that the (1,0°) pulse is emitted
# to the left.

# The LightSource records which pulse it emitted in which direction.
type LightSource = ref object of Mechanism
  cL: VectorChannel
  cR: VectorChannel
  cLog: IntChannel

proc newLightSource(cL, cR: VectorChannel; cLog: IntChannel): LightSource =
  ## Create a LightSource object.
  result = LightSource(cL: cL, cR: cR, cLog: cLog)
  result.init()

method run(ls: LightSource) =
  ## Emit a light pulse.
  let n = ls.randState.rand(1)
  let val = 90 * n.toFloat
  ls.cL[].send initVector(1, val)
  ls.cR[].send initVector(1, 90 - val)
  ls.cLog[].send(n)


####################################################################################################
# PolarizingBeamSplitter

# A polarizing beam splitter takes a plane-polarized light pulse
# and splits it into two plane-polarized pulses. The directions of
# polarization of the two output pulses are determined solely by the
# angular setting of the beam splitter—NOT by the angle of the
# original pulse. However, the amplitudes of the output pulses
# depend on the angular difference between the impinging light pulse
# and the beam splitter.

# Each beam splitter is designed to select randomly between one of
# two angle settings. It records which setting it selected (by
# putting that information into one of its output queues).

type PolarizingBeamSplitter = ref object of Mechanism
  cS: VectorChannel
  cS1: VectorChannel
  cS2: VectorChannel
  cLog: IntChannel
  angles: array[2, float]

proc newPolarizingBeamSplitter(cS, cS1, cS2: VectorChannel; cLog: IntChannel;
                               angles: array[2, float]): PolarizingBeamSplitter =
  ## Create a PolarizingBeamSplitter object.
  result = PolarizingBeamSplitter(cS: cS, cS1: cS1, cS2: cS2, cLog: cLog, angles: angles)
  result.init()

method run(pbs: PolarizingBeamSplitter) =
  ## Split a light pulse into two pulses. One of output pulses
  ## may be visualized as the vector projection of the input pulse
  ## onto the direction vector of the beam splitter. The other
  ## output pulse is the difference between the input pulse and the
  ## first output pulse: in other words, the orthogonal component.
  let angleSetting = pbs.randState.rand(1)
  pbs.cLog[].send(angleSetting)

  let angle = pbs.angles[angleSetting]
  assert angle >= 0 and angle < 90

  let v = pbs.cS[].recv()
  let v1 = projection(v, initVector(1, angle))
  let v2 = v - v1

  pbs.cS1[].send(v1)
  pbs.cS2[].send(v2)


####################################################################################################
# LightDetector.

# Our light detector is assumed to work as follows: if a
# uniformly distributed random number between 0 and 1 is less than
# or equal to the intensity (square of the amplitude) of an
# impinging light pulse, then the detector outputs a 1, meaning the
# pulse was detected. Otherwise it outputs a 0, representing the
# quiescent state of the detector.

type LightDetector = ref object of Mechanism
  cIn: VectorChannel
  cOut: IntChannel

proc newLightDetector(cIn: VectorChannel; cOut: IntChannel): LightDetector =
  ## Create a LightDetector object.
  result = LightDetector(cIn: cIn, cOut: cOut)
  result.init()

method run(ld: LightDetector) =
  ## When a light pulse comes in, either detect it or do not.
  let pulse = ld.cIn[].recv()
  let intensity = pulse.magnitude^2
  ld.cOut[].send ord(ld.randState.rand(1.0) <= intensity)


####################################################################################################
# DataSynchronizer.

# The data synchronizer combines the raw data from the logs and
# the detector outputs, putting them into dictionaries of
# corresponding data.

type DataSynchronizer = ref object of Mechanism
  cLogS: IntChannel
  cLogL: IntChannel
  cLogR: IntChannel
  cDetectedL1: IntChannel
  cDetectedL2: IntChannel
  cDetectedR1: IntChannel
  cDetectedR2: IntChannel
  cDataOut: DataChannel

proc newDataSynchronizer(cLogS, cLogL, cLogR, cDetectedL1, cDetectedL2,
                         cDetectedR1, cDetectedR2: IntChannel;
                         cDataOut: DataChannel): DataSynchronizer =
  ## Create a DataSynchronizer object.
  result = DataSynchronizer(cLogS: cLogs, cLogL: cLogL, cLogR: cLogR, cDetectedL1: cDetectedL1,
                            cDetectedL2: cDetectedL2, cDetectedR1: cDetectedR1,
                            cDetectedR2: cDetectedR2, cDataOut: cDataOut)

method run(ds: DataSynchronizer) =
  ## This method does the synchronizing.
  ds.cDataOut[].send DataOut {"logS": ds.cLogS[].recv(),
                              "logL": ds.cLogL[].recv(),
                              "logR": ds.cLogR[].recv(),
                              "detectedL1": ds.cDetectedL1[].recv(),
                              "detectedL2": ds.cDetectedL2[].recv(),
                              "detectedR1": ds.cDetectedR1[].recv(),
                              "detectedR2": ds.cDetectedR2[].recv()}.toTable

proc saveRawData(filename: string; data: seq[DataOut]) =

  proc saveData(f: File) =
    f.write $data.len, '\n'
    for entry in data:
      for i, key in ["logS", "logL", "logR", "detectedL1",
                     "detectedL2", "detectedR1", "detectedR2"]:
        f.write $entry[key]
        f.write if i == 6: '\n' else: ' '

  if filename != "-":
    var f = open(filename, fmWrite)
    f.saveData()
    f.close()
  else:
    stdout.saveData()


when isMainModule:

  if paramCount() notin [1, 2]:
    quit "Usage: $1 NUM_PULSES [RAW_DATA_FILE]" % getAppFilename(), QuitFailure
  let numPulses = paramStr(1).parseInt()
  let rawDataFilename = if paramCount() == 2: paramStr(2) else: "-"

  # Angles commonly used in actual experiments. Whatever angles you
  # use have to be zero degrees or placed in Quadrant 1. This
  # constraint comes with no loss of generality, because a
  # polarizing beam splitter is actually a sort of rotated
  # "X". Therefore its orientation can be specified by any one of
  # the arms of the X. Using the Quadrant 1 arm simplifies data
  # analysis.

  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)

  # Channels used for communications between the threads.
  const MaxSize = 100_000
  var
    cLogS = newIntChannel()
    cLogL = newIntChannel()
    cLogR = newIntChannel()
    cL = newVectorChannel()
    cR = newVectorChannel()
    cL1 = newVectorChannel()
    cL2 = newVectorChannel()
    cR1 = newVectorChannel()
    cR2 = newVectorChannel()
    cDetectedL1 = newIntChannel()
    cDetectedL2 = newIntChannel()
    cDetectedR1 = newIntChannel()
    cDetectedR2 = newIntChannel()
    cDataOut = newDataChannel()

  # Objects that will run in the various threads.
  var
    lightSource = newLightSource(cL, cR, cLogS)
    splitterL = newPolarizingBeamSplitter(cL, cL1, cL2, cLogL, AnglesL)
    splitterR = newPolarizingBeamSplitter(cR, cR1, cR2, cLogR, AnglesR)
    detectorL1 = newLightDetector(cL1, cDetectedL1)
    detectorL2 = newLightDetector(cL2, cDetectedL2)
    detectorR1 = newLightDetector(cR1, cDetectedR1)
    detectorR2 = newLightDetector(cR2, cDetectedR2)
    sync = newDataSynchronizer(cLogS, cLogL, cLogR, cDetectedL1, cDetectedL2,
                               cDetectedR1, cDetectedR2, cDataOut)

  # Open channels.
  for c in [cLogs, cLogL, cLogR, cDetectedL1, cDetectedL2, cDetectedR1, cDetectedR2]:
    c[].open(MaxSize)
  for c in [cL, cR, cL1, cL2, cR1, cR2]:
    c[].open(MaxSize)
  cDataOut[].open(MaxSize)

  # Threads.
  var lightsourceThread,
      splitterLThread,
      splitterRThread,
      detectorL1Thread,
      detectorL2Thread,
      detectorR1Thread,
      detectorR2Thread,
      syncThread: Thread[Mechanism]

  # Start the threads.
  lightsourceThread.createThread(start, Mechanism(lightSource))
  splitterLThread.createThread(start, Mechanism(splitterL))
  splitterRThread.createThread(start, Mechanism(splitterR))
  detectorL1Thread.createThread(start, Mechanism(detectorL1))
  detectorL2Thread.createThread(start, Mechanism(detectorL2))
  detectorR1Thread.createThread(start, Mechanism(detectorR1))
  detectorR2Thread.createThread(start, Mechanism(detectorR2))
  syncThread.createThread(start, Mechanism(sync))

  var data: seq[DataOut]
  for i in 1..numPulses:
    data.add cDataOut[].recv()
  saveRawData(rawDataFilename, data)

  # Stop the apparatus.
  lightSource.stop()
  splitterL.stop()
  splitterR.stop()
  detectorL1.stop()
  detectorL2.stop()
  detectorR1.stop()
  detectorR2.stop()
  sync.stop()

  # Wait for thread terminations.
  joinThread(lightsourceThread)
  joinThread(splitterLThread)
  joinThread(splitterRThread)
  joinThread(detectorL1Thread)
  joinThread(detectorL2Thread)
  joinThread(detectorR1Thread)
  joinThread(detectorR2Thread)
  joinThread(syncThread)

  # Close the channels.
  for c in [cLogs, cLogL, cLogR, cDetectedL1, cDetectedL2, cDetectedR1, cDetectedR2]:
    c[].close()
  for c in [cL, cR, cL1, cL2, cR1, cR2]:
    c[].close()
  cDataOut[].close()
</syntaxhighlight>

{{out}}
Here is an example of the output of the analysis of data produced by the simulator:
<pre>
   light pulse events      100000

    correlation coefs
           0.0° 22.5°   -0.709952
           0.0° 67.5°   +0.707978
          45.0° 22.5°   +0.699308
          45.0° 67.5°   +0.707870

        CHSH contrast   +2.825108
  2*sqrt(2) = nominal   +2.828427
           difference   -0.003319

       CHSH violation   +0.825108</pre>

=={{header|ObjectIcon}}==

This implementation is a bit simpler than the Python.

<syntaxhighlight lang="objecticon">
#!/bin/env -S oiscript

#
# 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 io
import ipl.random
import util

$encoding UTF-8

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

global rng

# I am shamelessly using global variables. You could put all these
# inside an object or a record, however, and pass them around that
# way.
global logS, logL, logR
global detectionsL1, detectionsL2
global detectionsR1, detectionsR2

procedure main (args)
  # Because main is run only once, there is no TECHNICAL reason to use
  # "initial" in main.  But it documents that this is initialization.
  initial
  {
    # Choose the RNG you prefer:

    rng := PCG32 ()
    #rng := MersenneTwister64 ()
    #rng := MersenneTwister32 ()

    # Use the traditional Icon RNG to randomize our actual RNG.
    randomize ()
    every 1 to ?1000 do
      rng.real()
  }

  case *args of
  {
    1 : sally_forth (integer (args[1]), "-") | stop (&why)
    2 : sally_forth (integer (args[1]), args[2]) | stop (&why)
    default :
    {
      io.write ("Usage: ", &progname, " num_events [-|raw_data_file]")
      exit (1)
    }
  }
end

procedure sally_forth (num_events, output_file)
  local detectorL1, detectorL2
  local detectorR1, detectorR2
  local splitterL, splitterR

  initial
  {
    logS := list (num_events)
    logL := list (num_events)
    logR := list (num_events)
    detectionsL1 := list (num_events)
    detectionsL2 := list (num_events)
    detectionsR1 := list (num_events)
    detectionsR2 := list (num_events)
  }

  detectorL1 := Light_Detector (detectionsL1)
  detectorL2 := Light_Detector (detectionsL2)
  detectorR1 := Light_Detector (detectionsR1)
  detectorR2 := Light_Detector (detectionsR2)

  splitterL := Beam_Splitter (angleL1, angleL2, logL,
                              detectorL1, detectorL2)
  splitterR := Beam_Splitter (angleR1, angleR2, logR,
                              detectorR1, detectorR2)

  # The light source runs the whole show.
  light_source (num_events, logS, splitterL, splitterR)

  write_raw_data (output_file)

  return
end

procedure write_raw_data (output_file)
  local f, i
  if type (output_file) ~== "string" then
  {
    f := output_file
    f.write(*logS)
    every i := 1 to *logS do
      f.write(logS[i], " ", logL[i], " ", logR[i], " ",
              detectionsL1[i], " ", detectionsL2[i], " ",
              detectionsR1[i], " ", detectionsR2[i])
  }
  else if output_file == "-" then
    write_raw_data (FileStream.stdout)
  else
  {
    f := open (output_file, "w") | stop (&why)
    write_raw_data (f) | stop (&why)
    f.close ()
  }
  return
end

procedure light_source (num_events, setting_log, splitterL, splitterR)
  # The light source is a loop that runs a finite number of times and
  # transmits angle data to the two polarizing beam splitters, through
  # their "receive()" methods. It records its randomly chosen setting
  # in a log.
  #
  # I recommend you do NOT think of the light source as "photons". To
  # do so is to invite prejudgments that hamper visualization of the
  # simulation. Think of the light source as two pairs of orthogonally
  # plane-polarized light beams, along with a system of mirrors and
  # shutters.
  #
  # But, even if the light source were "photons", there is no cause
  # for objection, if the simulation records the "photons" values in a
  # log and uses those values later. Remember, physicists claim
  # "Bell’s Theorem" proves that "hidden variables" theories are
  # impossible. To presume a SIMULATION cannot know what the "hidden
  # variables" of a "photon" are is to presume the conclusion: that no
  # "hidden variables" theory is possible!
  #
  # On the other hand, for our simulation we MUST presume otherwise,
  # because what WE are trying to show is that, if you presume there
  # ARE "hidden variables", which a SIMULATION certainly could take
  # into account, then we get the CHSH contrast predicted by quantum
  # mechanics. And this is SUPPOSED to be impossible to do!
  #
  # Which is a blatant contradiction of John Bell’s conclusion. If
  # even ONE counterexample to a "theorem" exists, then THERE IS NO
  # THEOREM. Here is a counterexample. Therefore Bell MUST have made
  # an error.
  #
  # Short answer to "What error did Bell make?": Bell presumed his a
  # and b were THE "hidden variables". This is wrong. Bell's a and b
  # must be treated as FUNCTIONS of the "hidden variables". One simply
  # does not presume that such variables are independent and not
  # functions. To do so is a symptom of "innumeracy".
  #
  # Now, the "hidden variables" obviously share a common origin in the
  # light source--and so, just as obviously, a and b SHOULD BE tightly
  # correlated with each other. They share the same "hidden
  # variables".
  #
  # That these facts escape notice among the supposedly greatest
  # "geniuses" of our society is a profound SOCIAL conundrum. And it
  # means that "photons" almost surely do have "hidden variables" that
  # mainstream physicists simply are not looking for. They have given
  # up trying.
  #
  # This all affects computer programmers: why expend so much resource
  # trying to make a "quantum computer" out of "superposition states",
  # and having computer science and computer engineering graduate
  # students study "quantum logic", if there is nothing "magical"
  # about "superposition states"? Would one stake their future on an
  # inherently decimal digital computer made with 10-voltage logic
  # instead of saturated transistors? You simply do not stake your
  # future on such an idea.

  local i, rand_0_or_1

  every i := 1 to num_events do
  {
    rand_0_or_1 := rng.range(2) - 1
    setting_log[i] := rand_0_or_1
    if rand_0_or_1 = 0 then
    {
      splitterL.receive(i, 0.0)
      splitterR.receive(i, 90.0)
    }
    else
    {
      splitterL.receive(i, 90.0)
      splitterR.receive(i, 0.0)
    }
  }
  return
end

class Beam_Splitter ()
  # A polarizing beam splitter receives light-source angles via its
  # "receive()" method (and assumes the light pulses have amplitude
  # 1.0). A beam splitter angle is randomly selected and recorded in a
  # log. Amplitudes of output pulses are computed, and transmitted to
  # the "receive()" methods of light detector objects.

  public angle1, angle2, setting_log
  public detector1, detector2

  public new (phi1, phi2, log, detec1, detec2)
    angle1 := phi1
    angle2 := phi2
    setting_log := log
    detector1 := detec1
    detector2 := detec2
    return
  end

  public receive (i, source_angle)
    local my_angle, relative_angle, rand_0_or_1
    local amplitude1, amplitude2

    static pi_over_180
    initial
    {
      pi_over_180 := Math.atan (1.0) / 45.0
    }

    rand_0_or_1 := rng.range(2) - 1
    setting_log[i] := rand_0_or_1
    my_angle := (if rand_0_or_1 = 0 then angle1 else angle2)
    relative_angle := my_angle - source_angle
    amplitude1 := abs (Math.cos (pi_over_180 * relative_angle))
    amplitude2 := abs (Math.sin (pi_over_180 * relative_angle))
    detector1.receive(i, amplitude1)
    detector2.receive(i, amplitude2)
    return
  end
end

class Light_Detector ()
  # A light detector receives light pulse data through its "receive()"
  # method, and records its results in the "detections" list.

  private detections

  public new (detections_list)
    detections := detections_list
    return
  end

  public receive (i, amplitude)
    local intensity, randnum
    intensity := amplitude * amplitude
    randnum := rng.real()
    detections[i] := (if randnum <= intensity then 1 else 0)
    return
  end
end
</syntaxhighlight>

{{out}}
A randomized 10000-event run as analyzed by the Object Icon analyzer.
<pre>

   light pulse events       10000

    correlation coefs
               0° 23°   -0.696414
               0° 68°   +0.702837
              45° 23°   +0.686932
              45° 68°   +0.733537

        CHSH contrast   +2.819720
  2*sqrt(2) = nominal   +2.828427
           difference   -0.008707

       CHSH violation   +0.819720

</pre>

=={{header|Phix}}==
See [[Simulated_optics_experiment/Data_analysis#Phix]] which includes the generation phase, and for javascript compatability has file save/load commented out.

=={{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 simulation in what might be the most efficient way possible. It uses very small constant space. However, one probably would want to use a better random number generator. And you have to recompile to use a different initial seed, or to generate a different number of lines, if you want different output.

<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 writes results to output as they are
# calculated.

# Random numbers are obtained from the DLARAN function of LAPACK,
# translated below into Ratfor. This is a multiplicative congruential
# generator, and so not especially good for doing
# simulations. Substitute something else if you plan to use this
# program for serious work :)

define(numlines, 100000)

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

function dlaran (iseed)
  #
  # DLARAN returns a random real number from a uniform (0,1)
  # distribution. See
  # https://netlib.org/lapack/explore-html/d4/dff/dlaran_8f_source.html
  # for full details.
  #
  double precision dlaran
  integer iseed(1:4)
  integer m1, m2, m3, m4
  parameter (m1 = 494, m2 = 322, m3 = 2508, m4 = 2549)
  double precision one
  parameter (one = 1.0d+0)
  integer ipw2
  double precision r
  parameter (ipw2 = 4096, r = one / ipw2)
  integer it1, it2, it3, it4
  double precision rndout
  intrinsic dble, mod
10 continue
  #
  # multiply the seed by the multiplier modulo 2**48
  #
  it4 = iseed( 4 )*m4
  it3 = it4 / ipw2
  it4 = it4 - ipw2*it3
  it3 = it3 + iseed( 3 )*m4 + iseed( 4 )*m3
  it2 = it3 / ipw2
  it3 = it3 - ipw2*it2
  it2 = it2 + iseed( 2 )*m4 + iseed( 3 )*m3 + iseed( 4 )*m2
  it1 = it2 / ipw2
  it2 = it2 - ipw2*it1
  it1 = it1 + iseed( 1 )*m4 + iseed( 2 )*m3 + iseed( 3 )*m2 _
            + iseed( 4 )*m1
  it1 = mod( it1, ipw2 )
  #
  # return updated seed
  #
  iseed(1) = it1
  iseed(2) = it2
  iseed(3) = it3
  iseed(4) = it4
  #
  # convert 48-bit integer to a real number in the interval (0,1)
  #
  rndout = r*( dble( it1 )+r*( dble( it2 )+r*( dble( it3 ) _
                +r*( dble( it4 ) ) ) ) )
  if (rndout == 1.0d+0)
    {
      # If a real number has n bits of precision, and the first n bits
      # of the 48-bit integer above happen to be all 1 (which will
      # occur about once every 2**n calls), then DLARAN will be
      # rounded to exactly 1.0. Since DLARAN is not supposed to return
      # exactly 0.0 or 1.0 (and some callers of DLARAN, such as
      # CLARND, depend on that), the statistically correct thing to do
      # in this situation is simply to iterate again. N.B. the case
      # DLARAN = 0.0 should not be possible.
      goto 10
    }
  #
  dlaran = rndout
end

function rndnum ()   # Random in (0,1). (Neither 0 nor 1 is returned.)
  implicit none
  double precision rndnum
  double precision dlaran
  integer iseed(1:4)
  save iseed
  data iseed / 1, 1, 1, 1 /
  rndnum = dlaran (iseed)
end

subroutine ltsrc (logval, thetaL, thetaR) # Light source.
  implicit none
  integer logval
  double precision thetaL, thetaR
  double precision rndnum
  continue
  if (rndnum () < 0.5D0)
    {
      logval = 0
      thetaL = 0.0D0
      thetaR = 90.0D0
    }
  else
    {
      logval = 1
      thetaL = 90.0D0
      thetaR = 0.0D0
    }
end

subroutine pbs (logval, theta, phi1, phi2, _
                outC, outS)     # Polarizing beam splitter.
  implicit none
  integer logval
  double precision theta, phi1, phi2, outC, outS
  double precision rndnum
  double precision piv180
  continue
  piv180 = datan (1.0D0) / 45
  #
  # The sine might be negative, but we will be squaring it, anyway.
  #
  if (rndnum () < 0.5D0)
    {
      logval = 0
      outC = dcos (piv180 * (theta - phi1))
      outS = dsin (piv180 * (theta - phi1))
    }
  else
    {
      logval = 1
      outC = dcos (piv180 * (theta - phi2))
      outS = dsin (piv180 * (theta - phi2))
    }
end

subroutine detec (inp, outp)    # Light detector.
  implicit none
  double precision inp
  integer outp
  double precision rndnum
  continue
  if (rndnum () <= inp * inp)
    outp = 1
  else
    outp = 0
end

subroutine event ()
  implicit none
  integer logS, logL, logR
  integer pingL1, pingL2, pingR1, pingR2
  double precision thetaL, thetaR
  double precision valL1, valL2, valR1, valR2
  continue
  call ltsrc (logS, thetaL, thetaR)
  call pbs (logL, thetaL, angleL1, angleL2, valL1, valL2)
  call pbs (logR, thetaR, angleR1, angleR2, valR1, valR2)
  call detec (valL1, pingL1)
  call detec (valL2, pingL2)
  call detec (valR1, pingR1)
  call detec (valR2, pingR2)
  write (*,'(I1, 6(X, I1))') logS, logL, logR, _
    pingL1, pingL2, pingR1, pingR2
end

program OpSimS
  implicit none
  integer i
  continue
  #
  # '(I0)' will not work in f2c, but I do not feel like writing out
  # code for f2c that avoids printing leading spaces.
  #
  write (*,'(I0)') numlines
  for (i = 0; i != numlines; i = i + 1)
    call event ()
end
</syntaxhighlight>

{{out}}
The formatted output (of the raw data) will not work with f2c, so the program was compiled with gfortran, and analyzed by whatever analyzer my tab key came up with first.
<pre>

   light pulse events      100000

    correlation coefs
           0.0° 22.5°   -0.709578
           0.0° 67.5°   +0.704456
          45.0° 22.5°   +0.713273
          45.0° 67.5°   +0.717752

        CHSH contrast   +2.845060
  2*sqrt(2) = nominal   +2.828427
           difference   +0.016633

       CHSH violation   +0.845060

</pre>

=={{header|Python}}==

The program takes one or two arguments. The first argument is the number of light pulses the source will emit. The second argument, if present, is a filename for the raw data output. If the second argument is left out or is "-", output will be to standard output.

The simulation puts each simulated device (and also a "data synchronizer") in its own process. The processes communicate through queues.

<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
import multiprocessing as mp
import time as tm
from math import atan2, cos, floor, pi, radians, sin, sqrt
from random import randint, seed, uniform

def modulo(a, p):
    '''This is like the MODULO function of Fortran.'''
    return (a - (floor(a / p) * p));

class Vector:
    '''A simple implementation of Gibbs vectors, suited to our
    purpose.'''

    def __init__(self, magnitude, angle):
        '''A vector is stored in polar form, with the angle in
        degrees between 0 (inclusive) and 360 (exclusive).'''
        self.magnitude = magnitude
        self.angle = modulo(angle, 360)

    def __repr__(self):
        return ("Vector(" + str(self.magnitude)
                + "," + str(self.angle) + ")")

    @staticmethod
    def from_rect(x, y):
        '''Return a vector for the given rectangular coordinates.'''
        return Vector(sqrt(x**2 + y**2),
                      modulo (180 * pi * atan2 (y, x), 360))

    def to_rect(self):
        '''Return the x and y coordinates of the vector.'''
        x = self.magnitude * cos ((pi / 180) * self.angle)
        y = self.magnitude * sin ((pi / 180) * self.angle)
        return (x, y)

    @staticmethod
    def scalar_product(vector, scalar):
        '''Multiply a vector by a scalar, returning a new vector.'''
        return Vector(vector.magnitude * scalar, vector.angle)

    @staticmethod
    def dot_product(u, v):
        '''Return the dot product of two vectors.'''
        return (u.magnitude * v.magnitude
                * cos ((pi / 180) * (u.angle - v.angle)))

    @staticmethod
    def difference(u, v):
        '''Return the difference of two vectors.'''
        (xu, yu) = u.to_rect()
        (xv, yv) = v.to_rect()
        return Vector.from_rect(xu - xv, yu - yv)

    @staticmethod
    def projection(u, v):
        '''Return the projection of vector u onto vector v.'''
        scalar = Vector.dot_product(u, v) / Vector.dot_product(v, v)
        return Vector.scalar_product(v, scalar)

class Mechanism:
    '''A Mechanism represents a part of the experimental apparatus.'''

    def __call__(self):
        '''Run the mechanism.'''
        while True:
            self.run()
            # A small pause to try not to overtax the computer.
            tm.sleep(0.001)

    def run(self):
        '''Fill this in with what the mechanism does.'''
        raise NotImplementedError()

class LightSource(Mechanism):
    '''A LightSource occasionally emits oppositely plane-polarized
    light pulses, of fixed amplitude, polarized 90° with respect to
    each other. The pulses are represented by the vectors (1,0°) and
    (1,90°), respectively. One is emitted to the left, the other to
    the right. The probability is 1/2 that the (1,0°) pulse is emitted
    to the left.

    The LightSource records which pulse it emitted in which direction.

    '''

    def __init__(self, L, R, log):
        Mechanism.__init__(self)
        self.L = L         # Queue gets (1,0°) or (1,90°)
        self.R = R         # Queue gets (1,90°) or (1,0°)
        self.log = log     # Queue gets 0 if (1,0°) sent left, else 1.

    def run(self):
        '''Emit a light pulse.'''
        n = randint(0, 1)
        self.L.put(Vector(1, n * 90))
        self.R.put(Vector(1, 90 - (n * 90)))
        self.log.put(n)

class PolarizingBeamSplitter(Mechanism):
    '''A polarizing beam splitter takes a plane-polarized light pulse
    and splits it into two plane-polarized pulses. The directions of
    polarization of the two output pulses are determined solely by the
    angular setting of the beam splitter—NOT by the angle of the
    original pulse. However, the amplitudes of the output pulses
    depend on the angular difference between the impinging light pulse
    and the beam splitter.

    Each beam splitter is designed to select randomly between one of
    two angle settings. It records which setting it selected (by
    putting that information into one of its output queues).

    '''

    def __init__(self, S, S1, S2, log, angles):
        Mechanism.__init__(self)
        self.S = S              # Vector queue to read from.
        self.S1 = S1            # One vector queue out.
        self.S2 = S2            # The other vector queue out.
        self.log = log          # Which angle setting was used.
        self.angles = angles

    def run(self):
        '''Split a light pulse into two pulses. One of output pulses
        may be visualized as the vector projection of the input pulse
        onto the direction vector of the beam splitter. The other
        output pulse is the difference between the input pulse and the
        first output pulse: in other words, the orthogonal component.'''

        angle_setting = randint(0, 1)
        self.log.put(angle_setting)

        angle = self.angles[angle_setting]
        assert (0 <= angle and angle < 90)

        v = self.S.get()
        v1 = Vector.projection(v, Vector(1, angle))
        v2 = Vector.difference(v, v1)

        self.S1.put(v1)
        self.S2.put(v2)

class LightDetector(Mechanism):
    '''Our light detector is assumed to work as follows: if a
    uniformly distributed random number between 0 and 1 is less than
    or equal to the intensity (square of the amplitude) of an
    impinging light pulse, then the detector outputs a 1, meaning the
    pulse was detected. Otherwise it outputs a 0, representing the
    quiescent state of the detector.

    '''

    def __init__(self, In, Out):
        Mechanism.__init__(self)
        self.In = In
        self.Out = Out

    def run(self):
        '''When a light pulse comes in, either detect it or do not.'''
        pulse = self.In.get()
        intensity = pulse.magnitude**2
        self.Out.put(1 if uniform(0, 1) <= intensity else 0)

class DataSynchronizer(Mechanism):
    '''The data synchronizer combines the raw data from the logs and
    the detector outputs, putting them into dictionaries of
    corresponding data.

    '''

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

    def run(self):
        '''This method does the synchronizing.'''
        self.dataout.put({"logS" : self.logS.get(),
                          "logL" : self.logL.get(),
                          "logR" : self.logR.get(),
                          "detectedL1" : self.detectedL1.get(),
                          "detectedL2" : self.detectedL2.get(),
                          "detectedR1" : self.detectedR1.get(),
                          "detectedR2" : self.detectedR2.get()})

def save_raw_data(filename, data):
    def save_data(f):
        f.write(str(len(data)))
        f.write("\n")
        for entry in data:
            f.write(str(entry["logS"]))
            f.write(" ")
            f.write(str(entry["logL"]))
            f.write(" ")
            f.write(str(entry["logR"]))
            f.write(" ")
            f.write(str(entry["detectedL1"]))
            f.write(" ")
            f.write(str(entry["detectedL2"]))
            f.write(" ")
            f.write(str(entry["detectedR1"]))
            f.write(" ")
            f.write(str(entry["detectedR2"]))
            f.write("\n")
    if filename != "-":
        with open(filename, "w") as f:
            save_data(f)
    else:
        save_data(sys.stdout)

if __name__ == '__main__':

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

    seed()             # Initialize random numbers with a random seed.

    # Angles commonly used in actual experiments. Whatever angles you
    # use have to be zero degrees or placed in Quadrant 1. This
    # constraint comes with no loss of generality, because a
    # polarizing beam splitter is actually a sort of rotated
    # "X". Therefore its orientation can be specified by any one of
    # the arms of the X. Using the Quadrant 1 arm simplifies data
    # analysis.
    anglesL = (0.0, 45.0)
    anglesR = (22.5, 67.5)
    assert (all(0 <= x and x < 90 for x in anglesL + anglesR))

    # Queues used for communications between the processes. (Note that
    # the direction of communication is always forwards in time. This
    # forwards-in-time behavior is guaranteed by using queues for the
    # communications, instead of Python's two-way pipes.)
    max_size = 100000
    logS = mp.Queue(max_size)
    logL = mp.Queue(max_size)
    logR = mp.Queue(max_size)
    L = mp.Queue(max_size)
    R = mp.Queue(max_size)
    L1 = mp.Queue(max_size)
    L2 = mp.Queue(max_size)
    R1 = mp.Queue(max_size)
    R2 = mp.Queue(max_size)
    detectedL1 = mp.Queue(max_size)
    detectedL2 = mp.Queue(max_size)
    detectedR1 = mp.Queue(max_size)
    detectedR2 = mp.Queue(max_size)
    dataout = mp.Queue(max_size)

    # Objects that will run in the various processes.
    lightsource = LightSource(L, R, logS)
    splitterL = PolarizingBeamSplitter(L, L1, L2, logL, anglesL)
    splitterR = PolarizingBeamSplitter(R, R1, R2, logR, anglesR)
    detectorL1 = LightDetector(L1, detectedL1)
    detectorL2 = LightDetector(L2, detectedL2)
    detectorR1 = LightDetector(R1, detectedR1)
    detectorR2 = LightDetector(R2, detectedR2)
    sync = DataSynchronizer(logS, logL, logR, detectedL1, detectedL2,
                            detectedR1, detectedR2, dataout)

    # Processes.
    lightsource_process = mp.Process(target=lightsource)
    splitterL_process = mp.Process(target=splitterL)
    splitterR_process = mp.Process(target=splitterR)
    detectorL1_process = mp.Process(target=detectorL1)
    detectorL2_process = mp.Process(target=detectorL2)
    detectorR1_process = mp.Process(target=detectorR1)
    detectorR2_process = mp.Process(target=detectorR2)
    sync_process = mp.Process(target=sync)

    # Start the processes.
    sync_process.start()
    detectorL1_process.start()
    detectorL2_process.start()
    detectorR1_process.start()
    detectorR2_process.start()
    splitterL_process.start()
    splitterR_process.start()
    lightsource_process.start()

    data = []
    for i in range(num_pulses):
        data.append(dataout.get())
    save_raw_data(raw_data_filename, data)

    # Shut down the apparatus.
    logS.close()
    logL.close()
    logR.close()
    L.close()
    R.close()
    L1.close()
    L2.close()
    R1.close()
    R2.close()
    detectedL1.close()
    detectedL2.close()
    detectedR1.close()
    detectedR2.close()
    dataout.close()
    lightsource_process.terminate()
    splitterL_process.terminate()
    splitterR_process.terminate()
    detectorL1_process.terminate()
    detectorL2_process.terminate()
    detectorR1_process.terminate()
    detectorR2_process.terminate()
    sync_process.terminate()
</syntaxhighlight>

{{out}}
An example of output from the Python implementation of ''[[Simulated optics experiment/Data analysis]]'':
<pre>

   light pulse events      100000

    correlation coefs
               0° 22°   -0.707806
               0° 67°   +0.705415
              45° 22°   +0.712377
              45° 67°   +0.718882

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

       CHSH violation   +0.844480

</pre>

=={{header|Wren}}==
{{trans|Python}}
{{libheader|Wren-queue}}
{{libheader|Wren-vector}}
{{libheader|Wren-assert}}
Please note the following main differences from the Python version:

1. We use our own Vector2 class rather than translate Python's as this can already be instantiated using polar as well as Cartesian coordinates. However, it doesn't have a 'projection' method and so that has been coded separately.

2. Wren CLI is unable to spawn separate processes at present. Whilst concurrency is supported via co-operatively scheduled fibers, the VM is single threaded and therefore only one fiber can be run at a time. As access to shared state does not need to be synchronized and yielding to other fibers mid-execution is not required for this task, there is no point running each mechanism in a separate fiber and we therefore just run the entire script under the 'main' fiber.

3. Output is always saved to a file - there is no option to print it to the terminal.
<syntaxhighlight lang="wren">import "random" for Random
import "io" for File
import "os" for Process
import "./queue" for Queue
import "./vector" for Vector2
import "./assert" for Assert

var Rand = Random.new()
Vector2.useDegrees = true

/* Returns the projection of vector u onto vector v. */
var Projection = Fn.new { |u, v|
    var scalar = u.dot(v) / v.dot(v)
    return v * scalar
}

/* A mechanism represents a part of the experimental apparatus. */
class Mechanism {
    run() { Fiber.abort("Not implemented.") }
}

/*  A LightSource occasionally emits oppositely plane-polarized
    light pulses, of fixed amplitude, polarized 90° with respect to
    each other. The pulses are represented by the vectors (1,0°) and
    (1,90°), respectively. One is emitted to the left, the other to
    the right. The probability is 1/2 that the (1,0°) pulse is emitted
    to the left.

    The LightSource records which pulse it emitted in which direction.
*/
class LightSource is Mechanism {
    construct new(L, R, log) {
        _L = L      // Queue gets (1,0°) or (1,90°)
        _R = R      // Queue gets (1,90°) or (1,0°)
        _log = log  // Queue gets 0 if (1,0°) sent left, else 1.
    }

    // Emit a light pulse.
    run() {
        var n = Rand.int(2) // 0 or 1
        _L.push(Vector2.fromPolar(1, n * 90))
        _R.push(Vector2.fromPolar(1, 90 - (n * 90)))
        _log.push(n)
    }
}

/*  A polarizing beam splitter takes a plane-polarized light pulse
    and splits it into two plane-polarized pulses. The directions of
    polarization of the two output pulses are determined solely by the
    angular setting of the beam splitter—NOT by the angle of the
    original pulse. However, the amplitudes of the output pulses
    depend on the angular difference between the impinging light pulse
    and the beam splitter.

    Each beam splitter is designed to select randomly between one of
    two angle settings. It records which setting it selected (by
    putting that information into one of its output queues).
*/
class PolarizingBeamSplitter is Mechanism {
    construct new(S, S1, S2, log, angles) {
        _S = S              // Vector queue to read from.
        _S1 = S1            // One vector queue out.
        _S2 = S2            // The other vector queue out.
        _log = log          // The other vector queue out.
        _angles = angles
    }

    // Split a light pulse into two pulses. One of output pulses
    // may be visualized as the vector projection of the input pulse
    // onto the direction vector of the beam splitter. The other
    // output pulse is the difference between the input pulse and the
    // first output pulse: in other words, the orthogonal component.
    run() {
        var angleSetting = Rand.int(2) // 0 or 1
        _log.push(angleSetting)

        var angle = _angles[angleSetting]
        Assert.ok(0 <= angle && angle < 90)
        var v = _S.pop()
        var v1 = Projection.call(v, Vector2.fromPolar(1, angle))
        var v2 = v - v1

        _S1.push(v1)
        _S2.push(v2)
    }
}

/*  Our light detector is assumed to work as follows: if a
    uniformly distributed random number between 0 and 1 is less than
    or equal to the intensity (square of the amplitude) of an
    impinging light pulse, then the detector outputs a 1, meaning the
    pulse was detected. Otherwise it outputs a 0, representing the
    quiescent state of the detector.
*/
class LightDetector is Mechanism {
    construct new(In, Out) {
        _In = In
        _Out = Out
    }

    // When a light pulse comes in, either detect it or do not.
    run() {
        var pulse = _In.pop()
        var intensity = pulse.square
        _Out.push(Rand.float(1) <= intensity ? 1 : 0)
    }
}

/*  The data synchronizer combines the raw data from the logs and
    the detector outputs, putting them into dictionaries of
    corresponding data.
*/
class DataSynchronizer is Mechanism {
    construct new(logS, logL, logR,
                  detectedL1, detectedL2,
                  detectedR1, detectedR2,
                  dataout) {
        _logS = logS
        _logL = logL
        _logR = logR
        _detectedL1 = detectedL1
        _detectedL2 = detectedL2
        _detectedR1 = detectedR1
        _detectedR2 = detectedR2
        _dataout = dataout
    }

    // This method does the synchronizing.
    run() {
        _dataout.push(
            {
                "logS" : _logS.pop(),
                "logL" : _logL.pop(),
                "logR" : _logR.pop(),
                "detectedL1" : _detectedL1.pop(),
                "detectedL2" : _detectedL2.pop(),
                "detectedR1" : _detectedR1.pop(),
                "detectedR2" : _detectedR2.pop()
            }
        )
    }
}

var saveRawData = Fn.new { |filename, data|
    File.create(filename) { |f|
        f.writeBytes(data.count.toString)
        f.writeBytes("\n")
        for (entry in data) {
            f.writeBytes(entry["logS"].toString)
            f.writeBytes(" ")
            f.writeBytes(entry["logL"].toString)
            f.writeBytes(" ")
            f.writeBytes(entry["logR"].toString)
            f.writeBytes(" ")
            f.writeBytes(entry["detectedL1"].toString)
            f.writeBytes(" ")
            f.writeBytes(entry["detectedL2"].toString)
            f.writeBytes(" ")
            f.writeBytes(entry["detectedR1"].toString)
            f.writeBytes(" ")
            f.writeBytes(entry["detectedR2"].toString)
            f.writeBytes("\n")
        }
    }
}

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

// Angles commonly used in actual experiments. Whatever angles you
// use have to be zero degrees or placed in Quadrant 1. This
// constraint comes with no loss of generality, because a
// polarizing beam splitter is actually a sort of rotated
// "X". Therefore its orientation can be specified by any one of
// the arms of the X. Using the Quadrant 1 arm simplifies data
// analysis.
var anglesL = [0, 45]
var anglesR = [22.5, 67.5]
Assert.ok((anglesL + anglesR).all { |x| 0 <= x && x < 90 })

// Queues used for communications between the mechanisms. (Note that
// the direction of communication is always forwards in time.)
var logS = Queue.new()
var logL = Queue.new()
var logR = Queue.new()
var L = Queue.new()
var R = Queue.new()
var L1 = Queue.new()
var L2 = Queue.new()
var R1 = Queue.new()
var R2 = Queue.new()
var detectedL1 = Queue.new()
var detectedL2 = Queue.new()
var detectedR1 = Queue.new()
var detectedR2 = Queue.new()
var dataout = Queue.new()

// Mechanisms to be run.
var mechanisms = [
    LightSource.new(L, R, logS),
    PolarizingBeamSplitter.new(L, L1, L2, logL, anglesL),
    PolarizingBeamSplitter.new(R, R1, R2, logR, anglesR),
    LightDetector.new(L1, detectedL1),
    LightDetector.new(L2, detectedL2),
    LightDetector.new(R1, detectedR1),
    LightDetector.new(R2, detectedR2),
    DataSynchronizer.new(logS, logL, logR, detectedL1, detectedL2, detectedR1, detectedR2, dataout)
]

var data = []
for (i in 0...numPulses) {
    for (m in mechanisms) m.run()
    data.add(dataout.pop())
}
saveRawData.call(filename, data)</syntaxhighlight>

{{out}}
An example of output from the Wren implementation of ''[[Simulated optics experiment/Data analysis]]'':
<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>