C @(#)random.for	17.1.1.1 (ES0-DMD) 01/25/02 17:54:07
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT      (c) 1988 European Southern Observatory
C.IDENT          RANDOM
C.LANGUAGE       ESO-FOR
C.AUTHOR         Preben Grosbol (ESO-IPG)
C.KEYWORDS       random numbers, uniform, gaussian
C.COMMENT        Generate pseudo-random numbers, see Handbook of
C                Mathematical Functions, Abramowitz and Stegun, p.949.
C.VERSION        1.0    1988-Oct-08 : Creation,      PJG
C---------------------------------------------------------------------
C
      SUBROUTINE USEED(IIA,IIB,IIX)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.PURPOSE     define seeds and starting number for new sequence of
C             pseudo-random numbers
C.RETURN      None
C---------------------------------------------------------------------
      INTEGER   IIA                 ! IN:  Multiplier seed
      INTEGER   IIB                 ! IN:  Adder seed
      INTEGER   IIX                 ! IN:  Initial value
C
      INTEGER   IT,IA,IB,IX
      REAL      RT
      COMMON /RANCOM/IA,IB,IX,IT,RT
C
C     INITIATE SEED VARIABLES IN COMMON
C
      IT = 2**15
      RT = 1.0 / FLOAT(IT)
      IA = IIA
      IB = IIB
      IX = IIX
      RETURN
      END
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
      REAL FUNCTION UNIRAN()
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.PURPOSE     compute next number of a uniformed randon distribution
C             in the interval [0;1]
C.ALGORITHM   Modulo formula: X = a*X + b (mod T)
C.RETURN      Uniformed pseudo-random number in interval [0;1]
C---------------------------------------------------------------------
      INTEGER   IT,IA,IB,IX
      REAL      RT
      COMMON    /RANCOM/IA,IB,IX,IT,RT
C
C     COMPUTE NEXT RANDOM NUMBER
C
      IX = MOD( IA*IX+IB , IT )
      UNIRAN = FLOAT(IX) * RT
      RETURN
      END
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
      REAL FUNCTION GA1RAN()
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.PURPOSE     compute next number of a gaussian randon distribution
C             with mean 0.0 and sigma 1.0.
C.ALGORITHM   Sum of 12 uniformed distributed pseudo-random numbers
C.RETURN      Gaussian pseudo-random number with [0,1]
C---------------------------------------------------------------------
      INTEGER   I,IT,IA,IB,IX
      REAL      RT,GX
      REAL      UNIRAN
      COMMON    /RANCOM/IA,IB,IX,IT,RT
C
C     COMPUTE NEXT RANDOM NUMBER
C
      GX = -6.0
      DO 100, I = 1,12
         GX = GX + UNIRAN()
  100 CONTINUE
      GA1RAN = GX
      RETURN
      END
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
      REAL FUNCTION GA2RAN()
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.PURPOSE     compute next number of a gaussian randon distribution
C             with mean 0.0 and sigma 1.0.
C.ALGORITHM   Rational Approx. of Xp where Q(Xp) = p using
C             formula 26.2.23 in Mathmatical Functions, p.933.
C.RETURN      Gaussian pseudo-random number with [0,1]
C---------------------------------------------------------------------
      INTEGER   IT,IA,IB,IX
      REAL      RT,GX,UX
      REAL      UNIRAN
      COMMON    /RANCOM/IA,IB,IX,IT,RT
C
C     COMPUTE NEXT RANDOM NUMBER
C
      UX = UNIRAN()
      GX = UX
      IF (UX.GT.0.5) GX = 1.0 - UX
      IF (GX.LE.0.0) GX = 1.0E-20
      GX = SQRT( -2.0*ALOG(GX) )
      GX = GX - (2.515517+GX*(0.802853+GX*0.01328))/
     /          (1.0+GX*(1.432788+GX*(0.189269+GX*0.001368)))
      IF (UX.GT.0.5) THEN
         GA2RAN = GX
      ELSE
         GA2RAN = -GX
      ENDIF
      RETURN
      END