/******************************************************************************
 *  Compilation:  javac Torus.java
 *  Execution:    java Torus N
 *
 *  Estimate the center-of-mass of the intersection of a torus and two
 *  planes using Monte Carlo integration. We assume the density = 1.0.
 *
 *  Reference: Numerical Recipes in FORTRAN, p. 223
 *
 ******************************************************************************/

public class Torus {

    public static void main(String[] args) {
        int N = Integer.parseInt(args[0]);    // number of samples
        double sum   = 0.0;
        double sumx  = 0.0;
        double sumy  = 0.0;
        double sumz  = 0.0;
        double sum2  = 0.0;
        double sum2x = 0.0;
        double sum2y = 0.0;
        double sum2z = 0.0;
        double V     = 3.0 * 7.0 * 2.0;      // volume of enclosing region

        for (int i = 0; i < N; i++) {
            double x =  1.0 + 3.0 * Math.random();
            double y = -3.0 + 7.0 * Math.random();
            double z = -1.0 + 2.0 * Math.random();

            // is (x, y, z) in the torus?
            if (z*z + Math.pow((Math.sqrt(x*x + y*y) - 3), 2) <= 1) {
                sum   += 1.0;
                sumx  += x;
                sumy  += y;
                sumz  += z;
                sum2  += 1.0;
                sum2x += x * x;
                sum2y += y * y;
                sum2z += z * z;
            }
        }

        // print results
        double mass   = V * sum  / N;
        double massx  = V * sumx / N;
        double massy  = V * sumy / N;
        double massz  = V * sumz / N;
        double error  = V * Math.sqrt((sum2 /N - (sum /N)*(sum /N)) / N);
        double errorx = V * Math.sqrt((sum2x/N - (sumx/N)*(sumx/N)) / N);
        double errory = V * Math.sqrt((sum2y/N - (sumy/N)*(sumy/N)) / N);
        double errorz = V * Math.sqrt((sum2z/N - (sumz/N)*(sumz/N)) / N);

        StdOut.println("Mass     = " + mass  + " +- " + error);
        StdOut.println("Moment x = " + massx + " +- " + errorx);
        StdOut.println("Moment y = " + massy + " +- " + errory);
        StdOut.println("Moment z = " + massz + " +- " + errorz);

        StdOut.println("Center of mass = (" + (massx/mass) + ", " +
                                                  (massy/mass) + ", " +
                                                  (massz/mass) + ")");
    }
}