----------------------------------------------------------------------

-- A simulation of an EPR-B experiment, demonstrating that ‘Bell’s
-- theorem’ is fallacious.

-- References:

-- [1] J. S. Bell, ‘Bertlmann’s socks and the nature of reality’,
--     preprint, CERN-TH-2926 (1980).
--     http://cds.cern.ch/record/142461/ (Open access, CC BY 4.0)

-- [2] 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)

-- The following simulation is written in Ada—which, if carefully
-- used, seems to me a good language for conveying scientific
-- algorithms. Also, a free software Ada compiler is widely available:
-- GCC. Many readers can compile this program with just the following
-- simple command:

--   gnatmake eprb_simulation

-- Ada, simply put, seemed the best choice, among the programming
-- languages I felt comfortable using.

-- However, if you know another language, such as Fortran or Python,
-- it would be very useful for understanding this program if you were
-- to translate it from Ada into the other language.

----------------------------------------------------------------------

pragma ada_2022;
pragma wide_character_encoding (utf8);

with ada.text_io;
with ada.containers.ordered_sets;
with ada.numerics;
with ada.numerics.generic_elementary_functions;

procedure eprb_simulation is

  subtype range_1_to_2 is integer range 1 .. 2;

  type scalar is digits 15;
  type scalar_pair is array (range_1_to_2) of scalar;

  package scalar_elementary_functions is
    new ada.numerics.generic_elementary_functions (scalar);

  use ada.text_io;
  use ada.numerics;
  use scalar_elementary_functions;

----------------------------------------------------------------------

  -- For the sake of reproducibility, let us write our own random
  -- number generator. It will be a simple linear congruential
  -- generator. The author has used one like it, in quicksorts and
  -- quickselects to select the pivot. It is good enough for our
  -- purpose.

  type uint64 is mod 2 ** 64;

  -- The multiplier lcg_a comes from Steele, Guy; Vigna, Sebastiano
  -- (28 September 2021). ‘Computationally easy, spectrally good
  -- multipliers for congruential pseudorandom number generators’.
  -- arXiv:2001.05304v3 [cs.DS]

  lcg_a : constant uint64 := 16#F1357AEA2E62A9C5#;

  -- The value of lcg_c is not critical, but should be odd.

  lcg_c : constant uint64 := 1;

  seed  : uint64 := 0;

  function random_scalar
  return scalar
  with post => 0.0 <= random_scalar'result
                and random_scalar'result < 1.0 is
    randval      : scalar;
  begin
    -- Take the high 48 bits of the seed and divide it by 2**48.
    randval := scalar (seed / (2**16)) / scalar (2**48);

    -- Update the seed.
    seed := (lcg_a * seed) + lcg_c;

    return randval;
  end random_scalar;

  function random_1_or_2
  return range_1_to_2 is
  begin
    return (if random_scalar < 0.5 then 1 else 2);
  end random_1_or_2;

----------------------------------------------------------------------

  -- The simulation is of a two-arm Clauser-Aspect type
  -- experiment. See [2] for description of a closely related
  -- simulation. We will be referring to equation numbers in that
  -- paper.

  -- The objects of our simulation are as follows.

  -- There are light pulses that we shall call photons. These are
  -- randomly generated in pairs, with either a vertically polarized
  -- pulse to the left and a horizontally polarized pulse to the
  -- right, or vice versa. Call the polarization state φ. In our
  -- simulation, these photons are assumed also to carry ‘hidden
  -- variables’ called ξ. A photon generated with vertical
  -- polarization has ξ=ξa, and a photon with horizontal polarization
  -- has ξ=ξb. Thus the following is our representation:

  type ξ_value is (ξa, ξb, ξ_quiescent);
  type photon is
    record
      φ : scalar;       -- The current polarization angle, in radians.
      ξ : ξ_value;      -- The hidden variable. Either ξa or ξb.
    end record;

  -- The following are the only possible initial photons:

  photon_vertical   : constant photon := (φ => pi / 2.0, ξ => ξa);
  photon_horizontal : constant photon := (φ => 0.0, ξ => ξb);

  -- Later we will introduce special photodetectors that can measure
  -- the value of a photon’s ξ. The quiescent state of such a
  -- photodetector will be ξ_quiescent.

  -- Please note something carefully here: by including a ‘hidden
  -- variable’ in the simulation of a photon, we are contradicting an
  -- assumption in [1], which Bell attempts to justify on pages 15 and
  -- 16. To quote: ‘[I]t may be that it is not permissible to regard
  -- the experimental settings a and b in the analyzers as independent
  -- variables, as we did.’ That one can even imagine hidden variables
  -- and photodetectors able to measure the hidden variables is a
  -- direct counterexample to Bell’s assumption, and therefore proves
  -- beyond doubt that Bell’s assumption is fallacious! Therefore his
  -- argument as a whole fails.

  -- So we have done enough already. But let us proceed with the
  -- simulation, nonetheless. Doing so will help illustrate why simply
  -- assuming the existence of a photon with hidden variables
  -- contradicts Bell’s assumption—or, in other words, how Bell’s
  -- argument is circular. Bell does the following:

  --   1. He assumes that a particle must be an ‘independent
  --   variable’.

  --   2. This is equivalent to assuming that a particle has no
  --   ‘hidden variables’.

  --   3. He performs various mathematical operations from this
  --   assumption, concluding that a particle cannot have any hidden
  --   variables.

  -- It is a classic example of someone accidentally assuming the
  -- conclusion, which is an informal fallacy. There is no deductive
  -- error, and yet the conclusion remains unproven.

  -- (Side note: the contrast is striking, between our simple disproof
  -- by counterexample, and the sophistries at the top of page 16 of
  -- [1]. Notable is Bell’s reliance on his physical intuitions,
  -- employing no principles of logic and mathematics. Also he
  -- overcomplicates the matter, seemingly in a struggle to visualize
  -- the situation. All that one needs to visualize is a hidden
  -- variable and a detector capable of measuring it! Bell obviously
  -- did not realize he had to visualize a photon WITH a hidden
  -- variable, rather than the photon WITHOUT a hidden variable that
  -- he was used to imagining.)

  -- Next after photons, there are two polarizing beam splitters (PBS)
  -- (https://en.wikipedia.org/w/index.php?title=Polarizer&oldid=1152014034#Beam-splitting_polarizers).
  -- These redirect photons in proportion to the Law of Malus
  -- (https://en.wikipedia.org/w/index.php?title=Polarizer&oldid=1152014034#Malus's_law_and_other_properties),
  -- altering the photons’ polarization angles φ, but not their hidden
  -- variables ξ.

  type pbs is
    record
      θ : scalar;               -- The angle of beam splitting.
    end record;

  -- These are the different PBS settings, in radians, making up four
  -- pairs of (left,right) settings, chosen for maximum nominal CHSH
  -- contrast (=2√2):

  θ_left  : constant scalar_pair := (0.0, pi / 4.0);
  θ_right : constant scalar_pair := (pi / 8.0, (3.0 * pi) / 8.0);

  -- For computational convenience, and without loss of generality, we
  -- have normalized the PBS angles to the first quadrant. PBS
  -- geometry is actually a sort of ×, where one of the two slopes
  -- represents one of the channels, and the crossing slope represents
  -- the other.

  -- Finally, there are four photodetectors
  -- (https://en.wikipedia.org/w/index.php?title=Photodetector&oldid=1168137030).
  -- These are fancy new-model photodetectors: not only do they detect
  -- incident photons without fail, but also they tell you the value
  -- of the photon’s ξ. Their output is a ξ_value: the quiescent state
  -- is ξ_quiescent, and they go to ξa or ξb to register a photon
  -- detection. They can detect only one photon at a time, but the
  -- experimental arrangement is such that they will encounter at most
  -- one at a time.

----------------------------------------------------------------------

  -- Here is a procedure that generates a (left,right) pair of
  -- photons.

  procedure generate_photon_pair (photon_left  : out photon;
                                  photon_right : out photon) is
  begin
    if random_1_or_2 = 1 then
      photon_left := photon_vertical;
      photon_right := photon_horizontal;
    else
      photon_left := photon_horizontal;
      photon_right := photon_vertical;
    end if;
  end generate_photon_pair;

----------------------------------------------------------------------

  -- The following procedure decides towards which of two
  -- photodetectors a PBS will transmit a given photon. The output of
  -- the procedure is the outputs of the two photodetectors. The
  -- photon will have its angle of polarization altered, but the
  -- photodetectors cannot measure angle of polarization. The outputs
  -- are of what the photodetectors CAN measure: hidden variable
  -- values.

  procedure split_beam (photon0   : photon;
                        pbs0      : pbs;
                        detector1 : out ξ_value;
                        detector2 : out ξ_value) is
  begin
    detector1 := ξ_quiescent;
    detector2 := ξ_quiescent;
    if random_scalar < (cos (photon0.φ - pbs0.θ)) ** 2 then
      detector1 := photon0.ξ;
    else
      detector2 := photon0.ξ;
    end if;
  end split_beam;

----------------------------------------------------------------------

  -- The function one_event, below, simulates an event of the
  -- simulation: from generation of a pair of photons to their
  -- measurement by four photodetectors. The settings of the two PBS
  -- are chosen randomly and included in the output. The output is
  -- gathered into a record and given a unique identifier (so it can
  -- be stored in set objects).

  current_identifier : long_integer := 1;

  type event_record is
    record
      identifier        : long_integer;
      pbs_left_setting  : range_1_to_2;
      pbs_right_setting : range_1_to_2;
      detector_left_1   : ξ_value;
      detector_left_2   : ξ_value;
      detector_right_1  : ξ_value;
      detector_right_2  : ξ_value;
    end record;

  function one_event
  return event_record is
    ev           : event_record;
    pbs_left     : pbs;
    pbs_right    : pbs;
    photon_left  : photon;
    photon_right : photon;
  begin
    ev.identifier := current_identifier;
    current_identifier := current_identifier + 1;

    ev.pbs_left_setting := random_1_or_2;
    pbs_left.θ := θ_left(ev.pbs_left_setting);

    ev.pbs_right_setting := random_1_or_2;
    pbs_right.θ := θ_right(ev.pbs_right_setting);

    generate_photon_pair (photon_left, photon_right);

    split_beam (photon_left, pbs_left,
                ev.detector_left_1, ev.detector_left_2);
    split_beam (photon_right, pbs_right,
                ev.detector_right_1, ev.detector_right_2);

    return ev;
  end one_event;

----------------------------------------------------------------------

  -- We wish to simulate some number of events, and to gather the
  -- output records in a set. Here is where the set type is created.

  function "<" (a, b : event_record)
  return boolean is
  begin
    return (a.identifier < b.identifier);
  end "<";

  package event_record_sets is
    new ada.containers.ordered_sets (element_type => event_record);

----------------------------------------------------------------------

  -- The following function simulates a given number of events and
  -- gathers the output records into a set.

  function simulate_events (number_of_events : natural)
  return event_record_sets.set is
    events : event_record_sets.set;
  begin
    for i in 1 .. number_of_events loop
      event_record_sets.include (events, one_event);
    end loop;
    return events;
  end simulate_events;

----------------------------------------------------------------------

  -- Now comes analysis.

  procedure analyze (events : event_record_sets.set) is

    use event_record_sets;

    events_L1R1    : set;
    events_L1R2    : set;
    events_L2R1    : set;
    events_L2R2    : set;
    κ_L1R1, κ_L1R2 : scalar;
    κ_L2R1, κ_L2R2 : scalar;
    chsh_contrast  : scalar;
    chsh_nominal   : constant scalar := 2.0 * sqrt (2.0);

    procedure separate_into_four_subsets_by_pbs_settings is
      curs              : cursor;
      ev                : event_record;
    begin
      curs := first (events);
      while has_element (curs) loop
        ev := element (curs);
        case ev.pbs_left_setting is
          when 1 =>
            case ev.pbs_right_setting is
              when 1 => include (events_L1R1, ev);
              when 2 => include (events_L1R2, ev);
            end case;
          when 2 =>
            case ev.pbs_right_setting is
              when 1 => include (events_L2R1, ev);
              when 2 => include (events_L2R2, ev);
            end case;
        end case;
        curs := next (curs);
      end loop;
    end separate_into_four_subsets_by_pbs_settings;

    function estimate_correlation_coefficient (subset : set)
    return scalar is
      curs                   : cursor;
      ev                     : event_record;
      n, nLc, nLs, nRc, nRs  : scalar;
      cosL, sinL, cosR, sinR : scalar;
      cosLR, sinLR           : scalar;
    begin

      -- We use equation (2.4) of [2] to estimate cosines and sines
      -- for calculation of the correlation coefficient by equations
      -- (2.3) and (2.5).
      --
      -- We take account of the hidden variable in this way: on
      -- account of the identity cos(π/2-x)≡sin(x), if ξ=ξa it
      -- indicates that the cosine in split_beam is ‘actually’ a
      -- sine. For the hidden variable tells us that the photon had a
      -- vertical angle of polarization when it entered the PBS. We
      -- must, therefore, count this detection towards the correct
      -- cosine/sine estimate.

      n := 0.0;
      nLc := 0.0;
      nLs := 0.0;
      nRc := 0.0;
      nRs := 0.0;

      curs := first (subset);

      while has_element (curs) loop

        ev := element (curs);
        if ev.detector_left_1 /= ξ_quiescent then
          if ev.detector_left_1 = ξa then
            nLs := nLs + 1.0;
          else
            nLc := nLc + 1.0;
          end if;
        else
          if ev.detector_left_2 = ξb then
            nLs := nLs + 1.0;
          else
            nLc := nLc + 1.0;
          end if;
        end if;

        if ev.detector_right_1 /= ξ_quiescent then
          if ev.detector_right_1 = ξa then
            nRs := nRs + 1.0;
          else
            nRc := nRc + 1.0;
          end if;
        else
          if ev.detector_right_2 = ξb then
            nRs := nRs + 1.0;
          else
            nRc := nRc + 1.0;
          end if;
        end if;

        n := n + 1.0;

        curs := next (curs);

      end loop;

      cosL := sqrt (nLc / n);
      sinL := sqrt (nLs / n);
      cosR := sqrt (nRc / n);
      sinR := sqrt (nRs / n);

      -- Reference [2], Equation (2.3).
      cosLR := (cosR * cosL) + (sinR * sinL);
      sinLR := (sinR * cosL) - (cosR * sinL);

      -- Reference [2], Equation (2.5).
      return (cosLR ** 2) - (sinLR ** 2);

    end estimate_correlation_coefficient;

  begin
    separate_into_four_subsets_by_pbs_settings;
    κ_L1R1 := estimate_correlation_coefficient (events_L1R1);
    κ_L1R2 := estimate_correlation_coefficient (events_L1R2);
    κ_L2R1 := estimate_correlation_coefficient (events_L2R1);
    κ_L2R2 := estimate_correlation_coefficient (events_L2R2);

put_line(κ_L1R1'image);
put_line(κ_L1R2'image);
put_line(κ_L2R1'image);
put_line(κ_L2R2'image);

    chsh_contrast := κ_L1R1 - κ_L1R2 + κ_L2R1 + κ_L2R2;

    put_line ("CHSH contrast: " & chsh_contrast'image);
    put_line ("      nominal: " & chsh_nominal'image);
  end analyze;

----------------------------------------------------------------------

begin
  -- Simulate a reasonable number of events.
  analyze (simulate_events (10000));
end eprb_simulation;

----------------------------------------------------------------------
-- Some instructions for the Emacs text editor.
-- local variables:
-- mode: indented-text
-- tab-width: 2
-- end: