Simulated optics experiment/Simulator

Simulated optics experiment/Simulator is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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 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 Python example, which does quite a bit but is not very large, can serve as a reference implementation.

• A light source, which occasionally emits two pulses of 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 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 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 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 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 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. However, simulating the devices as independent processes, objects in an event loop, etc., has a psychological advantage, which may be important if you are trying to convince someone of your point.)

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:

1. The log recording of the light source setting.
2. The log recording of the setting of the polarizing beam splitter that is on the left.
3. The log recording of the setting of the polarizing beam splitter that is on the right.
4. The output of the light detector on the left that is receiving the "cosine" pulses.
5. The output of the light detector on the left that is receiving the "sine" pulses.
6. The output of the light detector on the right that is receiving the "cosine" pulses.
7. 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.

See also

• Simulated optics experiment/Data analysis

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.

#!/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()

Output:

An example of output from the Python implementation of Simulated optics experiment/Data analysis:

   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.850985
  2*sqrt(2) = nominal   +2.828427
           difference   +0.022558

       CHSH violation   +0.850985