beta_gen

C/C++ source

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