Simulated optics experiment/Simulator: Difference between revisions
Simulated optics experiment/Simulator (view source)
Revision as of 11:27, 6 February 2024
, 3 months ago→{{header|Wren}}: Changed to Wren S/H
m (→Extra "credit") |
m (→{{header|Wren}}: Changed to Wren S/H) |
||
| (12 intermediate revisions by 4 users not shown) | |||
Line 1:
{{draft task}}
{{alertbox|#DDDAE8|There is now an "extra credit" problem. One might wish to have a look at it.}}
===Introduction===
Line 15 ⟶ 17:
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]]
* 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.
Line 27 ⟶ 29:
* 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]]
The program must output its "raw data" results in the format shown here by example:
Line 50 ⟶ 52:
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
* 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.
Line 59 ⟶ 61:
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.
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}}==
Line 184 ⟶ 639:
</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}}==
Line 461 ⟶ 1,314:
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>
Line 834 ⟶ 1,908:
3. Output is always saved to a file - there is no option to print it to the terminal.
<syntaxhighlight lang="
import "io" for File
import "os" for Process
Line 1,059 ⟶ 2,133:
<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>
| |||