Below is source code for JoyChristianSimulation.java file
===============================

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

import java.util.Random;
import java.util.logging.Logger;
import javax.vecmath.Vector3d;

/**
 *
 * @author Chantal Roth, 2013
 */
public class JoyChristianSimulation {
    
    /* logging */           
    private static final Logger log =  Logger.getLogger("JCSimulation");    
    
    /* whether to pick a and b on a plane or on the sphere */
    private static boolean PICK_AB_VECTORS_ON_SPHERE = false;
    
    /* whether to pick e on a plane or on the sphere - default must be SPHERE! */
    private static boolean PICK_E_VECTORS_ON_SPHERE = true;

    /* maximum angle in degrees that we are considering */
    private static final double MAX_ANGLE = Math.PI;
    
    /* constant phase shifts */
    private static final double phi_op = +0.000;  // in radians
    private static final double phi_oq = +0.000;  // in radians
    private static final double phi_or = -1.517;  // in radians
    private static final double phi_os = +0.663;  // in radians
    
    /* nr experiments per angle */
    private static final int total_nr_experiments = 9000000;
    
    /* results */
    private Results results;
    
    /* utility to generate random values */
    private Random random;
    
    public JoyChristianSimulation() {
        random = new Random();
        results = new Results();
    }

    public void runExperiment() {               
        for (int i = 0; i  < total_nr_experiments; i++) {           
            oneRun( );                            
        }       
        PlotPanel.show(results);        
        results.printResults();
    }

    private void oneRun() {        
        Vector3D a, b, e;
        if (PICK_AB_VECTORS_ON_SPHERE ){
            a = randomVectorOnSphere();            
            b = randomVectorOnSphere();
        }
        else {
            a = new Vector3D(1.0, 0, 0);
            b = this.randomVectorOnPlane(Math.random()*MAX_ANGLE);            
        }
        if (PICK_E_VECTORS_ON_SPHERE) {
            e = randomVectorOnSphere();  
        }
        else {
            e = this.randomVectorOnPlane(Math.random()*2*Math.PI);
        }                        
        
        double theta = a.angle(b);
     
        Vector3d ae = cross(a, e);
        Vector3d be = cross(b, e);
       
        double eta_ae = angle(a,e);
        double eta_be = angle(b,e);            
        double eta_cross = angle(ae,be);
        
        double 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)); 
        double 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));
        
        double C_a1 = Math.cos(eta_ae + phi_op)/N_a; // ordinary channel; lambda = +1
        double C_a2 = Math.cos(eta_ae + phi_op + Math.PI)/N_a; // ordinary channel; lambda = -1
        
        double C_b1 = Math.cos(eta_be + phi_or + Math.PI/2)/N_b; // extraordinary channel; lambda = +1
        double C_b2 = Math.cos(eta_be + phi_or + 3*Math.PI/2)/N_b; // extraordinary channel; lambda = -1
        
        double 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));
        
        results.addExperimentCount(theta);
        
        results.addAEPlus(theta, throwDie(C_a1));
        results.addAEPlus(theta, throwDie(C_a2));
        
        results.addAEMinus(theta, throwDie(C_b1));
        results.addAEMinus(theta, throwDie(C_b2));
        
        results.addPlusResult(theta, throwDie(C_ab));
        results.addMinusResult(theta, throwDie(C_ab));
        
    }  
    
    public int throwDie(double C) {
        // probability is between 0 and 1
        double probability = Math.abs(C);
        // generates a random value beween 0 and 1
        double dieValue = random.nextDouble();
        // 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
        boolean throwPasses = dieValue < probability;
        int detection; // is either 0, 1 or -1
        if (throwPasses) {
            detection = (int)-Math.signum(C);
        }
        else {
            detection = 0;
        }
        return detection;
    }
    private double angle(Vector3d a, Vector3d b) {
         double eta = a.angle(b) ;         
         return eta;
    }
    
    private Vector3d cross(Vector3d a, Vector3d b) {       
        return new Vector3d(
                a.y*b.z - a.z*b.y, 
                a.z*b.x - a.x*b.z,
                a.x*b.y - a.y*b.x);
    }
   
    private Vector3D randomVectorOnSphere() {
        double x, y, z, length;
        do {
            x = random.nextGaussian();
            y = random.nextGaussian();
            z = random.nextGaussian();            
            length = x * x + y * y + z * z;            
        } while (length <= Double.MIN_NORMAL);
        double s = Math.sqrt(1.0 / length);
        return new Vector3D(x * s, y * s, z * s);
    }
    
    private Vector3D randomVectorOnPlane(double theta) {
        double x = Math.cos(theta);
        double y = Math.sin(theta);        
        return new Vector3D(x, y, 0);
    }
    
     public static void main(String[] args) {
        JoyChristianSimulation test = new JoyChristianSimulation();
        test.runExperiment();
    }     
     
     /* logging, mainly for debugging */
    private void log(String msg) {
       log.info(msg);
    }

    public Vector3D randomVectorRelativeTo(Vector3D a) {
        Vector3D b = new Vector3D(a.x, a.y, a.z);
        b.rotateX(Math.random()*2*Math.PI);
        b.rotateY(Math.random()*2*Math.PI);
        b.rotateZ(Math.random()*2*Math.PI);
        return b;
    }

    
}

Below is source code for Vector3D.java file
======================
/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package joy;

import javax.vecmath.Vector3d;

/**
 *
 * @author Chantal
 */
public class Vector3D extends Vector3d {

    public Vector3D(double x, double y, double z) {
        super(x, y, z);
    }

    public void rotateX(double angle) {
        rotateX(Math.cos(angle),  Math.sin(angle));
    }

    public void rotateY(double angle) {
        rotateY( Math.cos(angle),  Math.sin(angle));
    }

    public void rotateZ(double angle) {
        rotateZ( Math.cos(angle),  Math.sin(angle));
    }

    public void rotateX(double cosAngle, double sinAngle) {
        double newY = y * cosAngle - z * sinAngle;
        double newZ = y * sinAngle + z * cosAngle;
        y = newY;
        z = newZ;
    }

    public void rotateY(double cosAngle, double sinAngle) {
        double newX = z * sinAngle + x * cosAngle;
        double newZ = z * cosAngle - x * sinAngle;
        x = newX;
        z = newZ;
    }

    public void rotateZ(double cosAngle, double sinAngle) {
        double newX = x * cosAngle - y * sinAngle;
        double newY = x * sinAngle + y * cosAngle;
        x = newX;
        y = newY;
    }
}

Below is source code for Utils.java file
====================================
package joy;

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.logging.Logger;


/**
 *
 * @author Chantal Roth, 2013
 */
public class Utils {
    
     static final Logger log =  Logger.getLogger("Utils");
     
     public static boolean writeStringToFile(File f, String content, boolean append) {
        PrintWriter fout = null;
        try {
            fout = new PrintWriter(new BufferedWriter(new FileWriter(f, append)));
            fout.print(content);
            fout.flush();
            fout.close();
            return true;
        } catch (FileNotFoundException e) {
            warn("File " + f + " not found: "+e.getMessage());
        } catch (IOException e) {
            warn("IO Exception: "+e.getMessage());
        } finally {
            if (fout != null) {
                fout.flush();
                fout.close();
            }
        }
        return false;
    }
     private static void warn(String msg) {
       log.warning(msg);
    }
    
     
}

Below is source code for Results.java file
=============================================
package joy;

import java.io.File;
import java.util.logging.Logger;

/**
 *
 * @author Chantal Roth, 2013
 */
public class Results {
    static Logger log =  Logger.getLogger("Results");
    
    private static double MAX_ANGLE = 360.0;
    private static int PLUS_TYPE = 0;
    private static int MINUS_TYPE = 1;   
    private static int AEPLUS_TYPE = 2;
    private static int AENEG_TYPE = 3;
     private static int NEITHER_TYPE = 4;
    private static int TOTAL_TYPE = 5;
    
    private double[][] databuckets;
    
    private double deltaAngle = 1.0;
    private int buckets = (int) (MAX_ANGLE/deltaAngle);       
    private int totalcount;
         
    public static String DESCRIPTIONS[] =  {
        "+1 counts for C, both channels",
        "-1 counts for C, both channels",        
        "+1 counts for C, single channel",
        "-1 counts for C, single channel",
        "neither +1 nor -1",
        "Total number"
    };    
    
    public Results() {
        databuckets = new double [buckets+1][TOTAL_TYPE+1];
        totalcount = 0;
    }
    
    public void addAEPlus(double angleRad, double value) { 
        double angleDeg = Math.toDegrees(angleRad);
        int bucket = getBucket(angleDeg);            
        if (value > 0) databuckets[bucket][AEPLUS_TYPE]+=value;        
    }
     public void addAEMinus(double angleRad, double value) { 
        double angleDeg = Math.toDegrees(angleRad);
        int bucket = getBucket(angleDeg);            
        if (value < 0) databuckets[bucket][AENEG_TYPE]+=-value;        
    }
     
    public void addPlusResult(double angleRad, double value) {
        add(angleRad, value, true);        
    } 
    public void addMinusResult(double angleRad, double value) {
        add(angleRad, value, false);        
    } 
    public void addExperimentCount(double angleRad) {
        double angleDeg = Math.toDegrees(angleRad);
        int bucket = getBucket(angleDeg);   
        databuckets[bucket][TOTAL_TYPE]++;
        totalcount++;
    }
    private void add(double angleRad, double value, boolean plus) { 
        double angleDeg = Math.toDegrees(angleRad);
        int bucket = getBucket(angleDeg);            
        if (Math.abs(value) == 1) {
            if (value > 0) {
                if (plus) databuckets[bucket][PLUS_TYPE]++;
           }
            else {
                if (!plus) databuckets[bucket][MINUS_TYPE]++;
            }            
        }              
    }
      /** Convert angle in degrees to bucket nr */
    public int getBucket(double angleDeg) {                          
        int bucket = (int) (Math.round(Math.abs(angleDeg) / deltaAngle));        
        return bucket;
    }
   
    public double getProbability(double angleDeg, int which) {
        int bucket = getBucket(angleDeg);
        double totalnr = ((double)(databuckets[bucket][TOTAL_TYPE]));
        double prob = (double)databuckets[bucket][which]/totalnr;      
        return prob;
    }
    private void log(String msg) {
       log.info(msg);
    }
   
    public int getNrBuckets() {
        return buckets;
    }
    public int getTotalcount() {
        return totalcount;
    } 
   
    public void printResults() {
        log("Total count: "+totalcount);
        
        String s= "Total count, "+totalcount+"\n;"
                + "angle degrees, plus counts, neg counts, neither, total\n";
        for (int a = 0; a < MAX_ANGLE; a++ ) {
            int b = getBucket(a);
            s +=a+", "+databuckets[b][PLUS_TYPE]+", "+databuckets[b][MINUS_TYPE]
                    +", "+databuckets[b][NEITHER_TYPE]+", "+databuckets[b][TOTAL_TYPE]+"\n";
        }    
        Utils.writeStringToFile(new File("results.csv"), s, false);
    }
}

Below is source code for PlotPanel.java file
=================================
package joy;

import java.awt.Color;
import java.awt.Dimension;
import java.awt.Font;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.RenderingHints;
import java.awt.image.BufferedImage;
import java.io.File;
import java.text.DecimalFormat;
import java.util.logging.Logger;
import javax.imageio.ImageIO;
import javax.swing.JFrame;
import javax.swing.JPanel;

/**
 * simple class that draws the correlation result from the experiment
 * @author Chantal Roth, 2013
 */
public class PlotPanel extends JPanel {
    
    static final Logger log = Logger.getLogger("PlotPanel");    
    static double MAX_PLOT_ANGLE = 180.0;
    static final int BORDER = 40;  
    
    private Results results;

    public PlotPanel(Results results) {
        this.results = results;
    }

    public static void show(Results results) {
        PlotPanel pan = new PlotPanel(results);
        JFrame f = new JFrame("Classical Simulation of the EPR-Bohm Correlation");
        f.getContentPane().add(pan);
        f.setSize(new Dimension(800, 600));
        f.setVisible(true);
        pan.saveImage("simulationresult.png");
    }

    public BufferedImage getImage(int w, int h) {
        BufferedImage bimage = new BufferedImage(w, h, BufferedImage.TYPE_INT_RGB);
        this.paintComponent(bimage.createGraphics());
        return bimage;
    }

    @Override
    public void paintComponent(Graphics g) {
        Graphics2D g2 = (Graphics2D) g;
        g2.setRenderingHint(RenderingHints.KEY_FRACTIONALMETRICS, RenderingHints.VALUE_FRACTIONALMETRICS_ON);
        g2.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
        super.paintComponent(g2);

        
        // clear
        int w = this.getWidth();
        int h = this.getHeight();
        g.setColor(Color.white);
        g.fillRect(0, 0, w, h);

        // x from 0 to MAXANGLE
        double dx = (double) (w - 2 * BORDER) / MAX_PLOT_ANGLE;
        double dy = (double) (h - 2 * BORDER);
        int my = (int) ((h - BORDER));
        g.setColor(Color.black);
        g.drawLine(BORDER, my, w - BORDER, my);
        g.drawLine(BORDER, BORDER, BORDER, h - BORDER);

        g.setFont(new Font("sans serif", Font.PLAIN, 12));
        g.drawString("P = sin^2(eta/2)/2 and P = cos^2(eta/2)/2;   Total counts: " + results.getTotalcount(), BORDER + 20, 20);
        drawCoordinateSystem(g, dx, my, w, dy);
        drawData(g, dx, my, dy);
        drawReferenceCurves(g, dx, my, dy);
    }

    public void saveImage(String file) {
        BufferedImage img = getImage(800, 600);
        File f = new File(file);
        try {
            f.createNewFile();
            boolean ok = ImageIO.write(img, "png", f);
            if (!ok) {
                log("Faile to write to image " + f.getAbsolutePath());
            }
        } catch (Exception ex) {
            log("Unable to save image to " + f+": "+ex.getMessage());
        }
    }

    public void drawData(Graphics g, double dx, int my, double dy) {
        Color[] colors = {
            Color.blue, Color.red, Color.green.darker(), Color.gray.darker()
        };
        //g2.setStroke(new BasicStroke(2));
        double yscale = 1.0; // in case we want to scale the y axis 
        for (int which = 0; which < 4; which++) {
            g.setColor(colors[which]);
            g.drawString(Results.DESCRIPTIONS[which], BORDER + 20, which * 16 + 40);
            double xold = 0;
            for (double x = 0; x <= MAX_PLOT_ANGLE; x += 2.0) {
                double y = results.getProbability(x, which) * yscale;
                double yold = results.getProbability(xold, which) * yscale;
                int guix = (int) (x * dx) + BORDER;
                int guiy = (int) (my - y * dy);
                int guixold = (int) (xold * dx) + BORDER;
                int guiyold = (int) (my - yold * dy);
                g.drawLine(guix, guiy, guixold, guiyold);
                xold = x;
            }
        }
    }

    public void drawReferenceCurves(Graphics g, double dx, int my, double dy) {
        double xold = 0;
        for (double x = 0; x <= MAX_PLOT_ANGLE; x += 1.0) {
            g.setColor(Color.black);
            double rad = Math.toRadians(x);
            double y = Math.pow(Math.sin(rad / 2.0), 2) / 2.0;
            double radold = Math.toRadians(xold);
            double yold = Math.pow(Math.sin(radold / 2.0), 2) / 2.0;
            int guix = (int) (x * dx) + BORDER;
            int guiy = (int) (my - y * dy);
            int guixold = (int) (xold * dx) + BORDER;
            int guiyold = (int) (my - yold * dy);
            g.drawLine(guix, guiy, guixold, guiyold);
            xold = x;
        }

        xold = 0;
        for (double x = 0; x <= MAX_PLOT_ANGLE; x += 1.0) {
            g.setColor(Color.black);
            double rad = Math.toRadians(x);
            double y = Math.pow(Math.cos(rad / 2.0), 2) / 2.0;
            double radold = Math.toRadians(xold);
            double yold = Math.pow(Math.cos(radold / 2.0), 2) / 2.0;
            int guix = (int) (x * dx) + BORDER;
            int guiy = (int) (my - y * dy);
            int guixold = (int) (xold * dx) + BORDER;
            int guiyold = (int) (my - yold * dy);
            g.drawLine(guix, guiy, guixold, guiyold);
            xold = x;
        }
    }

    public void drawCoordinateSystem(Graphics g, double dx, int my, int w, double dy) {
        DecimalFormat f = new DecimalFormat("0.0");
        // draw x axis
        g.setColor(Color.black);
        for (double x = 0; x <= MAX_PLOT_ANGLE; x += 5.0) {
            int guix = (int) (x * dx) + BORDER;
            int guiy = (int) (my + 15);
            if (x % 5 == 0) {
                g.drawLine(guix, my - 3, guix, my + 3);
            }
            if (x % 45 == 0) {
                g.drawString("" + f.format(x), guix - 10, guiy);
            }
        }
        g.drawString("Angle eta_ab", w / 2 - BORDER, my + 30);
        // draw y axis      
        for (double y = 0; y <= 1.0; y += 0.1) {
            int guix = BORDER - 20;
            int guiy = (int) (my - y * dy);
            g.drawLine(guix + 10, guiy, BORDER, guiy);
            g.drawString("" + f.format(y), guix - 10, guiy + 5);
        }
    }

    private void log(String msg) {
        log.info(msg);
    }
}