diff options
| author | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-03-04 21:21:30 -0500 |
|---|---|---|
| committer | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-03-04 21:21:30 -0500 |
| commit | d54fe7c1f704a63824c5bfa0ece65245572e9b27 (patch) | |
| tree | afc52015ffc2c74e0266653eecef1c8ef8ba5d91 /src/slalib/gresid.f_pcm | |
| download | calfuse-d54fe7c1f704a63824c5bfa0ece65245572e9b27.tar.gz | |
Initial commit
Diffstat (limited to 'src/slalib/gresid.f_pcm')
| -rw-r--r-- | src/slalib/gresid.f_pcm | 73 |
1 files changed, 73 insertions, 0 deletions
diff --git a/src/slalib/gresid.f_pcm b/src/slalib/gresid.f_pcm new file mode 100644 index 0000000..b2c6f76 --- /dev/null +++ b/src/slalib/gresid.f_pcm @@ -0,0 +1,73 @@ + REAL FUNCTION sla_GRESID (S) +*+ +* - - - - - - - +* G R E S I D +* - - - - - - - +* +* Generate pseudo-random normal deviate ( = 'Gaussian residual') +* (single precision) +* +* Given: +* S real standard deviation +* +* The results of many calls to this routine will be +* normally distributed with mean zero and standard deviation S. +* +* The Box-Muller algorithm is used. This is described in +* Numerical Recipes, section 7.2. +* +* !!! Microsoft Fortran dependent - calls the RAN routine !!! +* !!! To seed the random-number generator, either call the !!! +* !!! Microsoft SEED routine specifying some INTEGER*2 !!! +* !!! seed or call the function sla_RANDOM specifying some !!! +* !!! REAL seed. !!! +* +* P.T.Wallace Starlink 1 April 1993 +* +* Copyright (C) 1995 Rutherford Appleton Laboratory +*- + + IMPLICIT NONE + + REAL S + + REAL X,Y,RV,R,W,GNEXT,G + LOGICAL FIRST + + SAVE GNEXT,RV,FIRST + DATA FIRST / .TRUE. / + + + +* Second normal deviate of the pair available? + IF (FIRST) THEN + +* No - generate two random numbers in range +/- 1 + 1 CONTINUE + CALL RANDOM(RV) !!! PC + X = 2.0*RV-1.0 + CALL RANDOM(RV) !!! PC + Y = 2.0*RV-1.0 + +* Try again if not in unit circle + R = X*X+Y*Y + IF (R.GE.1.0) GO TO 1 + +* Box-Muller transformation, generating two deviates + W = SQRT(-2.0*LOG(R)/MAX(R,1E-20)) + GNEXT = X*W + G = Y*W + +* Set flag to indicate availability of next deviate + FIRST = .FALSE. + ELSE + +* Return second deviate of the pair & reset flag + G = GNEXT + FIRST = .TRUE. + END IF + +* Scale the deviate by the required standard deviation + sla_GRESID = G*S + + END |