double gasdev() {
        double R, v1, v2;
        if (calculateGASDEV) {
            // if we don't have an extra deviate handy
            do {
                // pick two uniform numbers in the square extending from -1 to +1
                // in each direction
                v1 = 2 * Math.random() - 1;
                v2 = 2 * Math.random() - 1;
                // check if v1 and v2 are in the unit circle
                R = v1 * v1 + v2 * v2;
            } while (R >= 1); // if v1 and v2 are not in the unit circle

            // make the Box-Muller transformation to get two normal deviates
            double fac = Math.sqrt((-1) * ((2 * Math.log(R)) / R));
            gasdevValue2 = v1 * fac; // one of the two normal deviates. gasdevValue2 is returned
            // next time this function is called
            gasdevValue1 = v2 * fac; // the other normal deviate, which will be returned at this time
            calculateGASDEV = false;
            return gasdevValue1;
        } else {
            calculateGASDEV = true;
            return gasdevValue2;
        }
    }