//---------------------------------------------------------------------------
#include <vcl.h>
#include <stddef.h>
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#pragma hdrstop

#include "globals.h"

//---------------------------------------------------------------------------
#pragma package(smart_init)

// frand() in (0,1)
float frand()
{
    return (((float)rand())/((float)RAND_MAX));
};



int gauss_help=0;
float gauss_savegaussdev=0;
const float PI=3.141592654;

/*
// returns values with N(0,1) (normal distribution, mi=0, sigma=1)
float Gauss1()
{
    float v1,v2,fac,r;
    if(gauss_help==0)
    {
        do
        {
            v1=2.0*frand()-1.0;
            v2=2.0*frand()-1.0;
            r=sqr(v1)+sqr(v2);
        }
        while (r>=1.0);
        fac=sqrt(-2.0*log(r)/r);
        gauss_savegaussdev=v1*fac;
        gauss_help=1;
        return v2*fac;
    }
    else
    {
        gauss_help=0;
        return gauss_savegaussdev;
    };
};
*/

// returns values with N(0,1) (normal distribution, mi=0, sigma=1)
float Gauss2()
{
    float a,b;
    if(gauss_help==0)
    {
        float frnd=frand();
        frnd=frnd>0?frnd:1e-10;
        a=sqrt(-2.0*log(frnd));
        b=2*PI*frand();
        gauss_savegaussdev=a*cos(b);
        gauss_help=1;
        return a*sin(b);
    }
    else
    {
        gauss_help=0;
        return gauss_savegaussdev;
    };
};

float getGauss(float mean, float sd)
{
    return Gauss2()*sd+mean;
};

float getGaussPosit(float mean, float sd)
{
    float v = Gauss2()*sd+mean;
    return max(v,0);
};

float getGaussPositButR(float mean, float sd,float R)
{
    float v = Gauss2()*sd+mean;
    v=(v>R)?R:v;
    return max(v,0);
};