// EPR events simulation using Bryan Sanctuary model.
// The method is predefined by the LHV when emitting and is applied in the _click function to produce the event.
// Depending on the LHV k and pk, the correlations are defined using polarization or coherence.
//
// The display of the correlation is done using the sum of the coherence correlations added to the sum of the polarization correlations.
// 
// The program produces a data file to display a graph using the gnuplot program.
// The command line for gnuplot is commented out in the code.
// Author Pierre Leroy pierrel5@free.fr

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

// ----------------------------------------
// Constants
#define PI 3.14159265358979323846

// used to compare values
#define EPSILON 0.0000001

// Count of simulated events in EPR test for each angle combination
#define N_EPR 10000

// ----------------------------------------
// RNG random number generator

static unsigned long n_seed = 123;               // RNG seed

#define RBITS 15
#define R1 (1 << RBITS)                          // normalize 1 value for random max
#define RMAX (R1 - 1)                            // random max value

// return random number in [0..1[ range
float rand1(void)
{
  n_seed = n_seed * 214013L + 2531011L;
  return ((n_seed >> 16) & RMAX)*(1.0f/R1);
}

// return random integer value -1 or 1
int rand_11(void)
{
  n_seed = n_seed * 214013L + 2531011L;
  return (int)(n_seed) >= 0 ? 1 : -1;
}

// ----------------------------------------
// Simulate emission at source 

// This structure contain hidden variables associated to a particle
struct lhv_t
{
  int k;
  int kp;
  double thRad;
};

// Simulate emission of the 2 particles.
// Define LHV for particles Alice and Bob.
void source_emit(struct lhv_t *pa, struct lhv_t *pb)
{
  double rand_theta = rand1()*2*PI;              // random angle 0..2PI
  int k = rand_11();                             // +/- 1

  // kp select randomly correlation type: polarization/coherence
  // kp = k select polarization
  // kp = -k select coherence.
  // ratio is 50%/50%
  int kp = rand1() < 0.5 ? k : -k;

  pa->thRad = rand_theta - PI/2;                 // theta particle Alice
  pa->k = k;
  pa->kp = kp;

  pb->thRad = rand_theta + PI/2;                 // theta particle Bob
  pb->k = -k;
  pb->kp = -kp;
}

// ----------------------------------------
// Simulate interaction with filter
// return event value +1 or -1
// inputs: 
//   - local particle LHV
//   - filter angle

int click_(struct lhv_t *lhv, double thARad)
{
  int evt = 0;
  double phase = 0;

  if (lhv->k == lhv->kp)
  {
    // polarization used to produce event
    phase = cos(thARad) * cos(lhv->thRad);
  }
  else
  {
    // coherence used to produce event
  
    // choose the larger spin axis
    double phasecc = cos(thARad * 2) * (cos(lhv->thRad) - sin(lhv->thRad));
    double phasess = sin(thARad * 2) * (sin(lhv->thRad) + cos(lhv->thRad));
    double phasem = (phasecc >= phasess) ? phasecc : phasess;
  
    if (fabs(phasem - phasecc) < EPSILON)
      phase = cos(thARad) * cos(lhv->thRad);
    else
    if (fabs(phasem - phasess) < EPSILON)
      phase = sin(thARad) * cos(lhv->thRad);
  }

  // return event
  return (phase >= 0) ? 1 : -1;
}

// ----------------------------------------
// EPR test

// Rotate Alice and Bob filter between 0..2*PI
// Memorize event correlation counters to draw correlation graph.

#define AN_STEP_DEG 5                            // step in filters angles difference to produce graph

// events counter for polarization
int ctr_EQP[360 + 1] = { 0 };                    // ++ -- events
int ctr_NEQP[360 + 1] = { 0 };                   // +- -+ events

// events counter for coherence
int ctr_EQC[360 + 1] = { 0 };                    // ++ -- events
int ctr_NEQC[360 + 1] = { 0 };                   // +- -+ events

// Do EPR test
void epr_test(void)
{
  int a;
  for (a=0; a<=360; a+=AN_STEP_DEG)              // filter a angle
  {
    double an_filter_a = (a*PI)/180.0;           // Filter angle at Alice in radians
    int ad;                                      // angle difference between filters
    for (ad=0; ad<=360; ad+=AN_STEP_DEG)
    {
      double an_filter_b = ((a + ad)*PI)/180.0;  // define b filter angle at Bob using a and ad

      // produce N_EPR events
      int n;
      for (n=0; n<N_EPR; n++)
      {
        struct lhv_t pa;                         // Declare Alice LHV
        struct lhv_t pb;                         // Declare Bob LHV
        int click_a;                             // Event at Alice
        int click_b;                             // Event at Bob

        // Simulate source. This initialize Alice and Bob LHV
        source_emit(&pa, &pb);

        // Simulate filter Alice, this define output event at Alice
        click_a = click_(&pa, an_filter_a);

        // simulate filter Bob
        click_b = click_(&pb, an_filter_b);

        // update counters
        if (pa.k == pa.kp)                       // k=kp: polarization event (k and kp sames for Bob, defines at emission)
        {
          if (click_a * click_b == 1)
            ctr_EQP[ad] += 1;                    // ++ or --
          else 
            ctr_NEQP[ad] += 1;                   // +- or -+
        }
        else                                     // k=-kp: coherence event
        {
          if (click_a * click_b == 1)
            ctr_EQC[ad] += 1;
          else
            ctr_NEQC[ad] += 1;
        }
      }
    }
  }
}

// ----------------------------------------
// Produce graph data file for gnuplot

// To draw graph with gnuplot use following command line into gnuplot console:

// Draw correlations command line:
// set xrange [0:360];plot "c:/sim_bm/gr_co.txt" using 1:2 with linespoint title "co. pol.","c:/sim_bm/gr_co.txt" using 1:3 with linespoint title "co. coh.","c:/sim_bm/gr_co.txt" using 1:4 with linespoint title "co. pol.+coh.","c:/sim_bm/gr_co.txt" using 1:5 with linespoint title "mustache",-cos(x*pi/180) title "-cosine" lt rgb "blue", 0
// Draw counters only command line:
// plot "c:/sim_bm/gr_co.txt" using 1:6 with linespoint title "EQP","c:/sim_bm/gr_co.txt" using 1:7 with linespoint title "NEQP","c:/sim_bm/gr_co.txt" using 1:8 with linespoint title "EQC","c:/sim_bm/gr_co.txt" using 1:9 with linespoint title "NEQC"

// get only polarization correlation from event counters
double get_EABP(int ad)
{
  double EQP  = ctr_EQP[ad];
  double NEQP = ctr_NEQP[ad];
  double EABP = (EQP + NEQP) ? (double)(EQP - NEQP) / (EQP + NEQP) : 0;
  return EABP;
}

// get get coherence correlation from event counters
double get_EABC(int ad)
{
  double EQC  = ctr_EQC[ad];
  double NEQC = ctr_NEQC[ad];
  double EABC = (EQC + NEQC) ? (double)(EQC - NEQC) / (EQC + NEQC) : 0;
  return EABC;
}

// code to draw exact "moustache" function
static double fmod1(double x) { return (x-floor(x)); }
static double tri(double x)   { return 1-4*fabs(fmod1(x/(2*PI))-0.5); }

// produce datas file to draw correlations graph
static void gen_gnuplot_file_co(const char *name)
{
  int ad, detect_sum = 0;
  double CHSH;

  FILE *f = fopen(name, "wb");
  if (!f)                                        // check file created
  {
    printf("Error: cannot create file %s\n", name);
    return;
  }

  // Save datas to graph file
  for (ad=0; ad<=360; ad+=AN_STEP_DEG)
  {
    // correlations
    double EABP = get_EABP(ad);
    double EABC = get_EABC(ad);
    double EAB = EABP + EABC;
    double ad_rad = ad*PI/180;
    double mustache = tri(ad_rad+PI) - cos(ad_rad);

    // save result for graph
    fprintf(f, "%d %.3f %.3f %.3f %.3f %d %d %d %d\n", ad, EABP, EABC, EAB, mustache,
      ctr_EQP[ad], ctr_NEQP[ad], ctr_EQC[ad], ctr_NEQC[ad]);
  }
  fclose(f);

  // print CHSH
  CHSH = fabs( (get_EABP(45)+get_EABC(45))*3 - (get_EABP(135)+get_EABC(135)) );
  printf("CHSH = %.3f\n", CHSH);
}

// main program entry
int main(void)
{
  printf("Start EPR simulation for N = %d events..\n", N_EPR);
  epr_test();
  printf("Simulation done.\n");
  gen_gnuplot_file_co("gr_co.txt");
}
