gaussian_gen

C/C++ source

0001 /*
0002  
0003     gaussian_gen.C    MEX file for a test program describing a centered random variable with a Gaussian distribution function
0004                     
0005     [s0] = gaussian_gen(xn,halfmode,fwmh,xreset0)
0006 
0007     INPUT:
0008         - xn: number of trials [1,1]
0009         - halfmode: 0 -> x = (-inf:+Inf), 1 -> x = (0:+Inf)
0010         - fwmh: full width at half maximum
0011         - xreset0: initialization number for the random generator [1,1]  
0012 
0013     OUTPUT:
0014     
0015         - x : random value [xn,1]
0016 
0017 
0018 
0019 by Y.PEYSSON CEA-DRFC <yves.peysson@cea.fr>
0020 */
0021 
0022 #include <math.h>
0023 #include "mex.h"
0024 
0025 /*input arguments*/
0026 
0027 #define    xn_IN        prhs[0]
0028 #define    halfmode_IN    prhs[1]
0029 #define    fwhm_IN        prhs[2]
0030 #define    xreset0_IN    prhs[3]           
0031 
0032 /* Output Arguments */
0033 
0034 #define    s0_OUT        plhs[0]
0035     
0036 #if !defined(MAX)
0037 #define    MAX(A, B)    ((A) > (B) ? (A) : (B))
0038 #endif
0039 
0040 #if !defined(MIN)
0041 #define    MIN(A, B)    ((A) < (B) ? (A) : (B))
0042 #endif
0043 
0044 
0045 
0046 
0047 /*
0048       Uniform random variate between 0 and 1
0049       Long period (> 2e18) random number generator 
0050       (see Numerical Recipes, 2nd Edition, p282)
0051       
0052 */
0053 
0054 static double ran2(idum)
0055 long int *idum;
0056 {
0057 
0058     long int j;
0059     long int k;
0060      static long int idum2 = 123456789;
0061      static long int iy = 0;
0062      static long int iv[32];
0063      double temp;
0064  
0065      double eps = 1.2e-7,rnmx,am,ndiv;
0066      long int im1 = 2147483563,im2 = 2147483399,imm1,ia1 = 40014,ia2 = 40692,iq1 = 53668,iq2 = 52774;
0067      long int ir1 = 12211,ir2 = 3791,ntab = 32;
0068  
0069      imm1 = im1 - 1L;
0070      am = (double)1.0/(double)im1;
0071      ndiv = (double)1.0+(double)imm1/(double)ntab;
0072      rnmx = (double)1.0 - eps;
0073  
0074      if (*idum <=0L) {
0075          if(-(*idum) < 1L) *idum = 1L;
0076          else *idum = -(*idum);
0077          idum2 = (*idum);
0078          for (j = ntab+7L;j>=0L;j--) {
0079              k = (long int)(*idum)/iq1;
0080              *idum = ia1*(*idum-k*iq1)-k*ir1;
0081              if(*idum<0L) *idum += im1;
0082              if (j < ntab) iv[j] = *idum;
0083          }
0084          iy = iv[0];
0085      }
0086      
0087      k = (long int)(*idum)/iq1;
0088      *idum = ia1*(*idum-k*iq1)-k*ir1;
0089      if (*idum < 0L) *idum += im1;
0090      k = (long int)(idum2/iq2);
0091      idum2 = ia2*(idum2-k*iq2)-k*ir2;
0092      if (idum2 < 0L) idum2 += im2;
0093      j = (long)(iy/(long int)ndiv);
0094      iy = iv[j]-idum2;
0095      iv[j] = *idum;
0096      if(iy < 1) iy += imm1;
0097      if ((temp = am*(double)iy) > rnmx) return(rnmx);
0098      else return(temp);
0099 }
0100 
0101 
0102 
0103 void tirgaussian(reset,halfmode,fwhm,res)
0104 long int *reset;
0105 double *halfmode,*fwhm,*res;
0106 {
0107     long int test = 0;
0108     double rand1,rand2,rand3;
0109         
0110     while (test == 0) {
0111         rand1 = -(fwhm[0]/(2.0*sqrt(2.0*log(2.0))))*(fwhm[0]/(2.0*sqrt(2.0*log(2.0))))*log(ran2(reset));
0112         rand2 = -(fwhm[0]/(2.0*sqrt(2.0*log(2.0))))*(fwhm[0]/(2.0*sqrt(2.0*log(2.0))))*log(ran2(reset)); 
0113         rand3 = 0.5*(rand1-1.0)*(rand1-1.0);
0114         if (rand3 <= rand2) {
0115             if (halfmode[0] == 1) {
0116                 res[0] = rand1;
0117             } else {
0118                 if (ran2(reset) < 0.5) {
0119                     res[0] = rand1;
0120                 } else {
0121                     res[0] = -rand1;
0122                 }
0123             }
0124             test = 1;
0125         }
0126     }
0127 }
0128 
0129 
0130 
0131 
0132 void gaussian(s0out,xn,halfmode,fwhm,xreset0)
0133 double s0out[];
0134 double *xn,*xreset0,*halfmode,*fwhm;
0135 {
0136     double xx;
0137     long int i,n,reset,reset0;
0138 /*
0139     Initialisation generateur v.a. uniform
0140 */
0141     reset0 = (long int)xreset0[0];
0142     reset = -reset0;
0143 
0144     n = (long int)xn[0];
0145     
0146     for (i=0L; i < n; i++) {
0147         tirgaussian(&reset,halfmode,fwhm,&xx);
0148         s0out[i] = xx;
0149     }
0150 }
0151 
0152 
0153 
0154 void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
0155 
0156     double *s0out;
0157     double *xn,*xreset0,*fwhm,*halfmode;
0158 
0159     /* Check for proper number of arguments */
0160 
0161     if (nrhs < 4) mexErrMsgTxt("The function gaussian_gen requires at least four input arguments.");
0162     if (nlhs > 1) mexErrMsgTxt("The function gaussian_gen requires up to one output argument.");
0163     
0164     xn = mxGetPr(xn_IN);
0165     halfmode = mxGetPr(halfmode_IN);
0166     fwhm = mxGetPr(fwhm_IN);
0167     xreset0 = mxGetPr(xreset0_IN);
0168             
0169     /* Create matrix for the return arguments */
0170 
0171     s0_OUT = mxCreateDoubleMatrix((mwSize)xn[0], 1, mxREAL);
0172 
0173     /* Assign pointers to the various parameters */
0174 
0175     s0out = mxGetPr(s0_OUT);
0176 
0177     /* Do the actual computations in a subroutine */
0178 
0179     gaussian(s0out,xn,halfmode,fwhm,xreset0);
0180     
0181     return;
0182 }