// (c) Copyright 1997 Peter Sandborn ALL RIGHTS RESERVED


//------------------------------------------------------------------------------
// distRandomNumbers.java:
// This class contains functions that return a random number with the specified
// distribution.
//------------------------------------------------------------------------------

import java.lang.Math;
import java.util.Random;

public class distRandomNumbers
{
    private MonteCarlo monte;
	private boolean haveNextNextGaussian;
	private double nextNextGaussian;

    // Constructors
    public distRandomNumbers(SetupButton p)
    {
        monte = p.monte;
		haveNextNextGaussian = false;
    }

    // Triangular distribution
    public double triangleDist(double mean,double low,double high)
    {
        if(low == high)  // avoids divide by zero
            return(mean);

        double x = ran1();
        double z = 2.0/(high-low);

        // Check left side of the distribution
        double result = -1.0;
        //double area_left_side = a1*(0.5*(mean*mean+low*low)-low*mean);
		double area_left_side = 0.5*(mean-low)*z;
        if(x != 0.0 && x <= area_left_side)  {
			double a1 = z/(mean-low);
            double A = a1/2.0;
            double B = -a1*low;
            double C = A*low*low - x;
            result = (-B + Math.sqrt(B*B - 4.0*A*C))/(2.0*A);
            if(result >= low && result <= high)
                return(result);
            else  {
                result = (-B - Math.sqrt(B*B - 4.0*A*C))/(2.0*A);
                return(result);
            }
        }  else  {
			double b1 = z/(high-mean);
            double A = -b1/2.0;
            double B = b1*high;
            double C = area_left_side - x + b1*mean*(0.5*mean - high);
            result = (-B + Math.sqrt(B*B - 4.0*A*C))/(2.0*A);
            if(result >= low && result <= high)
                return(result);
            else  {
                result = (-B - Math.sqrt(B*B - 4.0*A*C))/(2.0*A);
                return(result);
            }
        }
    }

    // Uniform distribution
    // Returns a uniformly distributed deviate between low and high
    public double uniformDist(double mean,double range)
    {
        double result = ran1();

        // Scale deviate
        result = result*(range);

        // Shift deviate to the proper mean
        result += mean-range/2.0;

        return(result);
    }

    // Lognormal distribution
    // Returns a lognormal distributed, positive, random
    // deviate with a most likely value of "value"
    // Numerical Recipes in C, page 215
    public double lognormalDist(double value,double stddev)
    {
        return(Math.exp(normalDist(value,stddev)));
    }

/*    // Normal distribution
    // Returns a normally distributed deviate with a mean of "value" and a
    // standard deviation of "stddev"
    // Numerical Recipes in C, page 217
    public double normalDist(double value,double stddev)
    {
        int iset = 0;
        double gset = 0.0;
        double fac;
        double r;
        double v1;
        double v2;
        double result = -1.0;


        if(iset == 0)  {
            do  {
                v1 = 2.0*ran1() - 1.0;
                v2 = 2.0*ran1() - 1.0;
                r = v1*v1*v2*v2;
            }  while(r >= 1.0);
            fac = Math.sqrt(-2.0*Math.log(r)/r);
            gset = v1*fac;
            iset = 1;
            result = v2*fac;
        }  else  {
            iset = 0;
            result = gset;
        }

        // Scale variance to specified standard deviation (variance = (standard deviation)^2
        result = result*stddev;

        // Shift deviate to the proper mean
        result += value;

        return(result);
    }  */
	public double normalDist(double value,double stddev)
    {
		double result;
		Random r = new Random();
		
        if (haveNextNextGaussian) {
            haveNextNextGaussian = false;
            result = nextNextGaussian;
		} else {
            double v1, v2, s;
            do { 
                    v1 = 2 * r.nextDouble() - 1;   // between -1.0 and 1.0
                    v2 = 2 * r.nextDouble() - 1;   // between -1.0 and 1.0
                    s = v1 * v1 + v2 * v2;
            } while (s >= 1 || s == 0);
            double multiplier = Math.sqrt(-2 * Math.log(s)/s);
            nextNextGaussian = v2 * multiplier;
            haveNextNextGaussian = true;
            result = v1 * multiplier;
		}
		// result is a normally distributed deviate with zero mean and unit variance

        // Scale variance to specified standard deviation (variance = (standard deviation)^2
        result = result*stddev;

        // Shift deviate to the proper mean
        result += value;

        return(result);
    }

    // Uniform distribution (between 0 and 1)
    // Returns a uniformly distributed deviate between 0.0 and 1.0
    // Numerical Recipes in C, page 210-211
    private double ran1()
    {
        double temp;
        int j;

        int M1 = 259200;
        int IA1 = 7141;
        int IC1 = 54773;
        double RM1 = 1.0/M1;
        int M2 = 134456;
        int IA2 = 8121;
        int IC2 = 28411;
        double RM2 = 1.0/M2;
        int M3 = 243000;
        int IA3 = 4561;
        int IC3 = 51349;

        if(monte.idum < 0 || monte.iff == 0)  {
            monte.iff = 1;
            monte.ix1 = (IC1 - monte.idum) % M1;
            monte.ix1 = (IA1 * monte.ix1 + IC1) % M1;
            monte.ix2 = monte.ix1 % M2;
            monte.ix1 = (IA1 * monte.ix1 + IC1) % M1;
            monte.ix3 = monte.ix1 % M3;
            for(j = 1; j <= 97; j++)  {
                monte.ix1 = (IA1 * monte.ix1 + IC1) % M1;
                monte.ix2 = (IA2 * monte.ix2 + IC2) % M2;
                monte.r[j] = (monte.ix1 + monte.ix2 * RM2) * RM1;
            }
            monte.idum = 1;
        }

        monte.ix1 = (IA1 * monte.ix1 + IC1) % M1;
        monte.ix2 = (IA2 * monte.ix2 + IC2) % M2;
        monte.ix3 = (IA3 * monte.ix3 + IC3) % M3;
        j = 1 + (int)((97 * monte.ix3)/M3);
        if(j > 97 || j < 1)
            ;  // This can't happen
        temp = monte.r[j];
        monte.r[j] = (monte.ix1 + monte.ix2 * RM2) * RM1;
        return(temp);
    }
}