/*
 * Tests Joy Christian's 3-sphere model for the EPR-Bohm correlation
 */

/**
 *
 * @author Chantal Roth-Korosten, 2013 (Java version)
 * Translated into javascript by Daniel Sabsay <danielsabsay@pacbell.net> 2013
 *
**/
 
  /* whether to Math.Pick a and b on a plane or on the sphere */
    PICK_AB_VECTORS_ON_SPHERE = false;    
  /* whether to Math.PIck e on a plane or on the sphere - default must be SPHERE! */
    PICK_E_VECTORS_ON_SPHERE = true;
  /* constant phase shifts */
    PHI_OP = 0.000;  // in radians
    PHI_OQ = 0.000;  // in radians
    PHI_OR = -1.517;  // in radians
    PHI_OS = 0.663;  // in radians
    PLUS_TYPE = 0;
    TOTAL_TYPE = 1;
    MINUS_TYPE = 2;
    NEITHER_TYPE = 3;
    
    max_angle = 2*Math.PI;
    deltaAngle = 1.0;  // in degrees
    buckets = parseInt(max_angle/deg2rad(deltaAngle));        
    databuckets = [];  // [PLUS_TYPE,TOTAL_TYPE]
    for (j=0;j<=buckets;j++) databuckets[databuckets.length] = [0,0]; 
  
  /*
  # of trials per experiment
  */
    number_of_trials = 5000000;
  
    document.getElementById("count").innerHTML = number_of_trials;

    totalcount = 0;
    for (i=0;i<number_of_trials;i++) oneTrial();                            
   
    drawResults();


function oneTrial() {
    if (PICK_AB_VECTORS_ON_SPHERE) {
        a = randVectSphere();            
        b = randVectSphere();
    } else {
        a = [1.0,0.0,0.0];
        b = randVectPlane(Math.random()*max_angle);            
    }
    if (PICK_E_VECTORS_ON_SPHERE) {
        e = randVectSphere();  
    } else {
        e = randVectPlane(2*Math.PI*Math.random());
    }
    
    eta_ab = angle(a,b);
    ae = cross(a,e);
    be = cross(b,e);  
    eta_ae = angle(a,e);
    eta_be = angle(b,e);            
    eta_cross = angle(ae,be);
    
    N_a = Math.sqrt(Math.cos(eta_ae+PHI_OP) * Math.cos(eta_ae+PHI_OP) + Math.sin(eta_ae+PHI_OQ) * Math.sin(eta_ae+PHI_OQ)); 
    N_b = Math.sqrt(Math.cos(eta_be+PHI_OR) * Math.cos(eta_be+PHI_OR) + Math.sin(eta_be+PHI_OS) * Math.sin(eta_be+PHI_OS));
    
/*    
    C_a1 = Math.cos(eta_ae+PHI_OP)/N_a; // ordinary channel; lambda = +1
    C_a2 = Math.cos(eta_ae+PHI_OP + Math.PI)/N_a; // ordinary channel; lambda = -1    
    C_b1 = Math.cos(eta_be+PHI_OR+Math.PI/2)/N_b; // extraordinary channel; lambda = +1
    C_b2 = Math.cos(eta_be+PHI_OR+3*Math.PI/2)/N_b; // extraordinary channel; lambda = -1
*/    
    C_ab = (-Math.cos(eta_ae+PHI_OP) * Math.cos(eta_be+PHI_OR) + Math.cos(eta_cross)
              * Math.sin(eta_ae + PHI_OQ) * Math.sin(eta_be + PHI_OS))/(N_a*N_b);
    
    addExperimentCount(eta_ab);
    addPlusResult(eta_ab,throwDie(C_ab));
    addPlusResult(eta_ab,throwDie(C_ab));
    
    addExperimentCount(2*Math.PI-eta_ab);
    addPlusResult(2*Math.PI-eta_ab,throwDie(C_ab));
    addPlusResult(2*Math.PI-eta_ab,throwDie(C_ab));

}  
    
function throwDie(C) {
    probability = Math.abs(C); // probability is between 0 and 1   
    dieValue = Math.random(); // generates a random value beween 0 and 1
    // so if the probability is 0, the die throw will never pass
    // if the probabiity is 1, it will always pass
    // if the probabilty is 0.5, it will pass 50% of the time
    throwPasses = dieValue < probability;
    if (throwPasses) {
        return (C<0)? 1:-1;
    } else {
        return 0;
    }
}
function angle(a,b) {
    dot = a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
    a_mag = a[0]*a[0]+a[1]*a[1]+a[2]*a[2];
    b_mag = b[0]*b[0]+b[1]*b[1]+b[2]*b[2];
    return Math.acos(dot/Math.sqrt(a_mag*b_mag));
}

function cross(a,b) {return [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]}
   
function randVectSphere() {
    x = gaussianRand();  // 3D unit vector
    y = gaussianRand();  // in random direction
    z = gaussianRand();
    s = 1.0 / Math.sqrt(x*x + y*y + z*z);
    return [x*s,y*s,z*s];
}
    
function randVectPlane(eta_ab) { return [Math.cos(eta_ab),Math.sin(eta_ab),0.0]; }
    
function gaussianRand() {  // Interval is -0.5 to +0.5
 r6 = Math.random()+Math.random()+Math.random()+Math.random()+Math.random()+Math.random();
 return r6/3 - 1;
}
function addPlusResult(angleRad,value) {
    add(angleRad,value,true);                 
} 
function addMinusResult(angleRad,value) {
    add(angleRad,value,false);                
} 
function addExperimentCount(angleRad) {
    angleDeg = rad2deg(angleRad);
    bucket = getBucket(angleDeg); 
    databuckets[bucket][TOTAL_TYPE]++;
    totalcount++;
}
function add(angleRad,value,plus) {
    angleDeg = rad2deg(angleRad);
    bucket = getBucket(angleDeg);            
    if (Math.abs(value) == 1) {
        if (value > 0) { // if plus
            databuckets[bucket][PLUS_TYPE]++;
        } else { // if !plus
            databuckets[bucket][PLUS_TYPE]--;
        }            
    }              
}
function getBucket(angleDeg) {
    return parseInt(Math.abs(angleDeg)/deltaAngle);        
}
function deg2rad(deg) { return deg*Math.PI/180; }
function rad2deg(rad) { return rad*180/Math.PI; }
   
function drawResults()  {
  var canvas = document.getElementById("canvas");
  if (null==canvas || !canvas.getContext) return;
  var ctx = canvas.getContext("2d");
  var x0 = 0.5 + 0.5*canvas.width;  // x0 pixels from left to x=0
  var y0 = 0.5 + 0.5*canvas.height; // y0 pixels from top to y=0
  var max_bucket = databuckets.length-1;
  var xx, yy, dx = canvas.width/max_bucket;
  ctx.strokeStyle = "red";
  ctx.lineWidth = .4;
  ctx.beginPath();

  yscale = max_bucket/4*canvas.height/number_of_trials;
  for (var i=0;i<max_bucket;i++) {
    xx = i*dx;
    yy = yscale*databuckets[i][PLUS_TYPE];
    ctx.lineTo(xx,y0-yy);
  }
  ctx.stroke();

  ctx.beginPath();
  ctx.strokeStyle = "green"; 
  ctx.moveTo(0,y0); ctx.lineTo(ctx.canvas.width,y0);  // X axis
  ctx.moveTo(x0,0); ctx.lineTo(x0,ctx.canvas.height);  // Y axis
  ctx.stroke();
  ctx.strokeStyle = "blue";
  ctx.font = "6pt Arial";
  ctx.strokeText("-180",0,y0-2);
  ctx.strokeText("+180",274,y0-2);
}