← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to gui/XFilter/src/random.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:
gui/XFilter/src/random.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 Massachusetss Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Corresponding 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
.VERSION
 
28
 090826         last modif
 
29
 
 
30
===========================================================================*/
 
31
 
 
32
 
 
33
#include <stdio.h>
 
34
#include <stdlib.h>
 
35
#include <math.h>
 
36
 
 
37
#ifndef PI   
 
38
#define PI 3.141592654  /* as in stroustrup */
 
39
#endif
 
40
 
 
41
/* -------------------------------------------------------------------
 
42
float   random_local()
 
43
 
 
44
purpose: draw random numbers between 0 and 1 under an uniform law
 
45
-------------------------------------------------------------------*/
 
46
 
 
47
#define M1 259200
 
48
#define IA1 7141
 
49
#define IC1 54773
 
50
#define RM1 (1.0/M1)
 
51
#define M2 134456
 
52
#define IA2 8121
 
53
#define IC2 28841
 
54
#define RM2 (1.0/M2)
 
55
#define M3 243000
 
56
#define IA3 4561
 
57
#define IC3 51349
 
58
 
 
59
float random_local(idum)
 
60
int *idum;
 
61
 
 
62
{
 
63
        static long ix1,ix2,ix3;
 
64
        static float r[98];
 
65
        float temp;
 
66
        static int iff=0;
 
67
        int j;
 
68
        void nrerror();
 
69
 
 
70
        if (*idum<0 || iff==0)
 
71
        {
 
72
                iff=1;
 
73
                ix1=(IC1-(*idum)) % M1;
 
74
                ix1=(IA1*ix1+IC1) % M1;
 
75
                ix2=ix1%M2;
 
76
                ix1=(IA1*ix1+IC1) % M1;
 
77
                ix3=ix1%M3;
 
78
                for(j=1;j<=97;j++)
 
79
                {
 
80
                        ix1=(IA1*ix1+IC1) % M1;
 
81
                        ix2=(IA2*ix2+IC2) % M2;
 
82
                        r[j]=(ix1+ix2*RM2)*RM1; 
 
83
                }
 
84
        *idum=1;
 
85
        }
 
86
        ix1=(IA1*ix1+IC1) % M1;
 
87
        ix2=(IA2*ix2+IC2) % M2;
 
88
        ix3=(IA3*ix3+IC3) % M3;
 
89
        j=1+((97*ix3)/M3);
 
90
        if(j>97 || j<1) {printf("je me suis plantee\n"); j=abs(*idum);};
 
91
        temp=r[j];
 
92
        r[j]=(ix1+ix2*RM2)*RM1;
 
93
        return temp;
 
94
}
 
95
 
 
96
/* -------------------------------------------------------------------
 
97
float   poisson()
 
98
 
 
99
purpose: draw numbers under an poisson law
 
100
-------------------------------------------------------------------*/
 
101
 
 
102
float poisson(xm,idum)
 
103
 
 
104
float xm;
 
105
int *idum;
 
106
 
 
107
{
 
108
static float sq,alxm,g,oldm=(-1.0);
 
109
float em,t ,y;
 
110
float random_local(),lngam();
 
111
 
 
112
if(xm<12.0)
 
113
        {
 
114
        if(xm!=oldm)
 
115
                {
 
116
                oldm=xm;
 
117
                g=exp(-1.*xm);
 
118
                };
 
119
        em=(-1);
 
120
        t=1.0;
 
121
        do
 
122
                {
 
123
                em+=1.0;
 
124
                t*=random_local(idum);
 
125
                }while(t>g);
 
126
        }
 
127
else
 
128
        {
 
129
        if(xm!=oldm)
 
130
                {
 
131
                oldm=xm;
 
132
                sq=sqrt(2.0*xm);
 
133
                alxm=log(xm);
 
134
                g=xm*alxm-lngam(xm+1.0);
 
135
                };
 
136
        do
 
137
                {
 
138
                do
 
139
                        {
 
140
                        y=tan(PI*random_local(idum));
 
141
                        em=sq*y+xm;
 
142
                        }while(em<0.0);
 
143
                em=floor(em);
 
144
                t=.9*(1.0+y*y)*exp(em*alxm-lngam(em+1.0)-g);
 
145
                }while(random_local(idum)>t);
 
146
        };
 
147
return (em);
 
148
}
 
149
 
 
150
/******************************************************/
 
151
 
 
152
float lngam(a)
 
153
 
 
154
float a;
 
155
 
 
156
{
 
157
double cof[6],stp,half=.5,one=1.,fpf=5.5,x,tmp,ser;
 
158
double pi;
 
159
double atemp;
 
160
float gama,gama1;
 
161
int i;
 
162
 
 
163
pi=3.141592654;
 
164
atemp=10.;
 
165
stp=2.50662827465;
 
166
cof[0]=76.18009173;
 
167
cof[1]=(-86.50532033);
 
168
cof[2]=24.01409822;
 
169
cof[3]=(-1.231739516);
 
170
cof[4]=0.120858003e-2;
 
171
cof[5]=(-0.536382e-5);
 
172
 
 
173
if(a==1)
 
174
        gama=1.;
 
175
else
 
176
        {
 
177
        if(a<1.)
 
178
                {
 
179
                atemp=1.;
 
180
                a=2.-a; 
 
181
                };
 
182
        x=a-one;
 
183
        tmp=x+fpf;
 
184
        tmp=(x+half)*log(tmp)-tmp;
 
185
        ser=one;
 
186
        for(i=0;i<6;i++)
 
187
                {
 
188
                x+=one;
 
189
                ser+=cof[i]/x;
 
190
                };
 
191
        gama1=tmp+log(stp*ser);
 
192
        gama=gama1;
 
193
        if(atemp<1.)
 
194
                {
 
195
                pi=pi*(a-1.);
 
196
                gama=log(pi/sin(pi))-gama;
 
197
                a=atemp;
 
198
                };
 
199
        };
 
200
return(gama);
 
201
}
 
202
 
 
203
/* -------------------------------------------------------------------
 
204
float   gauss()
 
205
 
 
206
purpose: draw numbers under an gaussian law
 
207
-------------------------------------------------------------------*/
 
208
 
 
209
float gauss(sig,idum)
 
210
 
 
211
float sig;
 
212
int *idum;
 
213
 
 
214
{
 
215
float x,z,r;
 
216
float random_local();
 
217
 
 
218
while( (z = pow(x=random_local(idum)-0.5,2.0) + pow(random_local(idum)-0.5,2.0) ) > 0.25 );
 
219
while ((r=random_local(idum))<=0.0);
 
220
 
 
221
return  sig*sqrt(-2.0*log(r)/z)*x;
 
222
 
 
223
}