← Back to branch summary

~ubuntu-branches/ubuntu/wily/eso-midas/wily-proposed

« back to all changes in this revision

Viewing changes to libsrc/math/randev.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
libsrc/math/randev.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1995-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
        Internet e-mail: midas@eso.org
 
21
        Postal address: European Southern Observatory
 
22
                        Data Management Division 
 
23
                        Karl-Schwarzschild-Strasse 2
 
24
                        D 85748 Garching bei Muenchen 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
.IDENT       randev.c
 
30
.LANGUAGE    C
 
31
.AUTHOR      P.Grosbol,  IPG/ESO
 
32
.ENVIRON     UNIX
 
33
.KEYWORDS    Random numbers, exponential, gaussian
 
34
.COMMENT     Algorithm is adapted from 'Numerical Recipes in C' p.214
 
35
.VERSION     1.0  1990-Nov-07 : Creation,  PJG
 
36
.VERSION     1.1  1994-May-28 : Change name of constant PI, PJG
 
37
.VERSION     1.2  1995-Mar-09 : Correct error + double gammln, PJG
 
38
 
 
39
090416          last modif
 
40
-----------------------------------------------------------------------*/
 
41
#include <math.h>                    /* Standard mathematical library  */
 
42
 
 
43
#include <proto_gen.h>
 
44
 
 
45
#define PICONST       3.141592654
 
46
 
 
47
double expdev(pseed)
 
48
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
49
.PURPOSE   compute random number taken from an Exponential disribution
 
50
           with unit mean.
 
51
.RETURN    exponential distributed random number
 
52
-----------------------------------------------------------------------*/
 
53
int    *pseed;                /* pointer to seed if <0 initiate seq.   */
 
54
{
 
55
  double  x;
 
56
 
 
57
  x = randm(pseed);
 
58
  return -log(x);
 
59
}
 
60
 
 
61
static  int     gflag = 0;
 
62
static  double  gval;
 
63
 
 
64
double gaudev(pseed)
 
65
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
66
.PURPOSE   compute random number taken from a Gaussian disribution
 
67
           with zero mean and unit standard deviation.
 
68
.RETURN    gaussian distributed random number
 
69
-----------------------------------------------------------------------*/
 
70
int    *pseed;                /* pointer to seed if <0 initiate seq.   */
 
71
{
 
72
  double  r, fac, v1, v2;
 
73
 
 
74
  if (gflag) { gflag = 0; return gval; }
 
75
  else {
 
76
     do {
 
77
        v1 = 2.0*randm(pseed) - 1.0;
 
78
        v2 = 2.0*randm(pseed) - 1.0;
 
79
        r  = v1*v1 + v2*v2;
 
80
      } while (1.0<=r);
 
81
     fac  = sqrt( -2.0*log(r)/r );
 
82
     gval = v1*fac;
 
83
     gflag = 1;
 
84
     return v2*fac;
 
85
   }
 
86
}
 
87
 
 
88
double gamdev(ia,pseed)
 
89
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
90
.PURPOSE   compute random number taken from a Gamma  disribution of
 
91
           interger order of 'ia'
 
92
.RETURN    gamma distributed random number
 
93
-----------------------------------------------------------------------*/
 
94
int    ia;                    /* integer order of gamma function       */
 
95
int    *pseed;                /* pointer to seed if <0 initiate seq.   */
 
96
{
 
97
  int   j;
 
98
  double am, e, s, v1, v2, x, y;
 
99
 
 
100
  if (ia < 1) return 0.0;
 
101
  if (ia < 6) {
 
102
     x = 1.0;
 
103
     for (j=0; j<ia; j++) x *= randm(pseed);
 
104
     x = -log(x);
 
105
   } 
 
106
  else {
 
107
     do {
 
108
        do {
 
109
           do {
 
110
              v1 = 2.0*randm(pseed) - 1.0;
 
111
              v2 = 2.0*randm(pseed) - 1.0;
 
112
            } while (v1*v1+v2*v2 > 1.0);
 
113
           y  = v2/v1;
 
114
           am = ia - 1;
 
115
           s  = sqrt(2.0*am+1.0);
 
116
           x  = s*y + am;
 
117
         } while (x <= 0.0);
 
118
        e = (1.0+y*y)*exp(am*log(x/am)-s*y);
 
119
      } while (randm(pseed) > e);
 
120
   }
 
121
 
 
122
  return x;
 
123
}
 
124
 
 
125
/***
 
126
 
 
127
static double sq, alxm, g, oldm=(-1.0);
 
128
 
 
129
double poidev(xm,pseed)
 
130
/.+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
131
.PURPOSE   compute random number taken from a Poisson disribution
 
132
           with mean 'xm'.
 
133
.RETURN    poisson distributed random number
 
134
-----------------------------------------------------------------------./
 
135
double     xm;                /. mean value of poisson distribution    ./
 
136
int        *pseed;            /. pointer to seed if <0 initiate seq.   ./
 
137
{
 
138
  double   em, t, y, gammln();
 
139
  
 
140
  if (xm < 12.0) {
 
141
     if (xm != oldm) {
 
142
        oldm = xm;
 
143
        g = exp(-xm);
 
144
      }
 
145
     em = -1;
 
146
     t = 1.0;
 
147
     do {
 
148
        em += 1.0;
 
149
        t  *= randm(pseed);
 
150
      } while (t > g);
 
151
   } else {
 
152
      if (xm != oldm) {
 
153
         oldm = xm;
 
154
         sq = sqrt(2.0*xm);
 
155
         alxm = log(xm);
 
156
         g = xm*alxm - gammln(xm+1.0);
 
157
       }
 
158
      do {
 
159
         do {
 
160
            y = tan(PICONST*randm(pseed));
 
161
            em = sq*y + xm;
 
162
          } while (em < 0.0);
 
163
         em = floor(em);
 
164
         t = 0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
 
165
       } while (randm(pseed) > t);
 
166
    }
 
167
  return em;
 
168
}
 
169
***/
 
170