cauchylorentz_gen

C/C++ source

0001 /*
0002  
0003     gaussian_gen.C    MEX file for a test program describing a centered random variable with a Cauchy/Lorentz distribution function
0004                     
0005     [s0] = cauchylorentz_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 tircauchylorentz(reset,halfmode,fwhm,res)
0104 long int *reset;
0105 double *halfmode,*fwhm,*res;
0106 {
0107     long int test = 0;
0108     double pi = 3.141592653526;
0109     double rand1;
0110         
0111     while (test == 0) {
0112         rand1 = ran2(reset);
0113         if (halfmode[0] == 1.0) {
0114             if (rand1 >= 0.5) {
0115                 res[0] = fwhm[0]*tan(pi*(rand1 - 0.5))/2;
0116                 test = 1;
0117             }
0118         } else {
0119             res[0] = fwhm[0]*tan(pi*(rand1 - 0.5))/2;
0120             test = 1;
0121         }
0122     }
0123 }
0124 
0125 
0126 void cauchylorentz(s0out,xn,halfmode,fwhm,xreset0)
0127 double s0out[];
0128 double *xn,*xreset0,*halfmode,*fwhm;
0129 {
0130     double xx;
0131     long int i,n,reset,reset0;
0132 /*
0133     Initialisation generateur v.a. uniform
0134 */
0135     reset0 = (long int)xreset0[0];
0136     reset = -reset0;
0137 
0138     n = (long int)xn[0];
0139     
0140     for (i=0L; i < n; i++) {
0141         tircauchylorentz(&reset,halfmode,fwhm,&xx);
0142         s0out[i] = xx;
0143     }
0144 }
0145 
0146 
0147 
0148 void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
0149 
0150     double *s0out;
0151     double *xn,*xreset0,*fwhm,*halfmode;
0152 
0153     /* Check for proper number of arguments */
0154 
0155     if (nrhs < 4) mexErrMsgTxt("The function cauchylorentz_gen requires at least four input arguments.");
0156     if (nlhs > 1) mexErrMsgTxt("The function cauchylorentz_gen requires up to one output argument.");
0157     
0158     xn = mxGetPr(xn_IN);
0159     halfmode = mxGetPr(halfmode_IN);
0160     fwhm = mxGetPr(fwhm_IN);
0161     xreset0 = mxGetPr(xreset0_IN);
0162             
0163     /* Create matrix for the return arguments */
0164 
0165     s0_OUT = mxCreateDoubleMatrix((mwSize)xn[0], 1, mxREAL);
0166 
0167     /* Assign pointers to the various parameters */
0168 
0169     s0out = mxGetPr(s0_OUT);
0170 
0171     /* Do the actual computations in a subroutine */
0172 
0173     cauchylorentz(s0out,xn,halfmode,fwhm,xreset0);
0174     
0175     return;
0176 }