This week only. Pastebin PRO Accounts Christmas Special! Don't miss out!Want more features on Pastebin? Sign Up, it's FREE!
Guest

Ree distrubution

By: a guest on Feb 26th, 2012  |  syntax: C  |  size: 12.29 KB  |  views: 27  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1.     #include <stdio.h>
  2.     #include <string.h>
  3.     #include <math.h>
  4.     // Initiallization routine for random number
  5.     void   RandomInitialise(int,int);
  6.     // Call to generate a random number between 0  and 1
  7.     double RandomUniform(void);
  8.     // Call to generate a Gaussian random number
  9.     double RandomGaussian(double,double);
  10.     // Call to generate a random integer between the two arguments
  11.     int    RandomInt(int,int);
  12.     // call to generate a random double between the two arguments
  13.     double RandomDouble(double,double);
  14.     FILE *ofp;
  15.     int main ()
  16.     {
  17.       int i,j, r = 10000,Bin[10000], l, Reeint[10000], N =100;
  18.      
  19.       double theta,phi,pi=3.14159;
  20.       double Rx[10000], Ry[10000], Rz[10000], b0=1.0, Ree[10000], temp=0, k=0.0, c=1, avgRee =0;  
  21.      
  22.       // First need to initiallize the random number by giving it 2 seeds
  23.       /* Initiallize the Random number by giving it 2 integers
  24.          The integers should be between 0 and 31328 and 0 and 30081 */
  25.      
  26.       RandomInitialise(1802,9373);
  27.  
  28.          ofp=fopen("data9.txt","w");
  29.  
  30.       for(i=0;i<r;i++)
  31.       {
  32.               Rx[0]=0;
  33.               Ry[0]=0;
  34.               Rz[0]=0;
  35.               temp=0;
  36.              // printf("New walk is...\n");
  37.               for(j=1;j<N+1;j++)
  38.               {
  39.                       theta = RandomDouble(0, 2*pi);
  40.                       phi = RandomDouble(0,pi);
  41.                       Rx[j] = Rx[j-1] + b0*cos(phi);
  42.                       Ry[j] = Ry[j-1] + b0*sin(phi)*sin(theta);
  43.                       Rz[j] = Rz[j-1] + b0*sin(phi)*cos(theta);
  44.                      
  45.                            
  46.                      
  47.                //      printf("x = %f, y = %f, z = %f\n", Rx[j], Ry[j], Rz[j]);
  48.                      
  49.      
  50.      
  51.               }
  52.              
  53.               Ree[i] = sqrt(Rx[N]*Rx[N] + Ry[N]*Ry[N] + Rz[N]*Rz[N]);
  54.              // printf("Ree is = %f\n", Ree[i]);
  55.  
  56.  
  57.  
  58.                           /*temp = Ree;
  59.               Ree = Ree+temp;*/
  60.              
  61.      
  62.       }
  63.           for(i=0; i<N; i++)
  64.                           {
  65.  
  66.                           avgRee = avgRee + Ree[i];
  67.  
  68.                           }
  69.                                 avgRee = avgRee/N;
  70.                         for(i=0; i<N; i++)
  71.                         {
  72.                                 printf("Ree Value = %f\n", Ree[i]);
  73.                         //      fprintf(ofp,"%f\n",Ree[i]);
  74.  
  75.  
  76.                         }
  77.                           printf("Average Ree is %f\n", avgRee);
  78.  
  79.                           fprintf(ofp,"Average Ree  = %f\n",avgRee);
  80.  
  81.                          
  82.  
  83.         for(j=0; j<r;j++)
  84.         {
  85.                 int tracker=0;
  86.                 for(i=0; i<r; i++)
  87.                 {
  88.                         if(Ree[i]>k && Ree[i]<c)
  89.                         {
  90.                                 tracker++;
  91.                         }              
  92.                 }
  93.                 k = k++;
  94.                 c = c++;
  95.                 Bin[j]=tracker;
  96.         }
  97.         for(i=0; i<r; i++)
  98.         {
  99.                 printf("Bins are...%d\n", Bin[i]);
  100.                 fprintf(ofp,"%d\n",Bin[i]);
  101.  
  102.         }
  103.  
  104.  
  105.  
  106.    
  107.      
  108.      
  109.       // Generate 2 random numbers (theta and pi)
  110.       // theta is between 0 and 2PI and phi is between 0 and PI
  111.      
  112.      
  113.      // printf("%f %f\n", theta, phi);
  114.      
  115.       // Now do the same but use the function RandomDouble
  116.      
  117.       theta = RandomDouble(0, 2*pi);
  118.       phi = RandomDouble(0,pi);
  119.      
  120.      // printf("%f %f\n", theta, phi);
  121.      
  122.      
  123.            
  124.     }
  125.      
  126.     #define FALSE 0
  127.     #define TRUE 1
  128.      
  129.     /*
  130.      
  131.      This Random Number Generator is based on the algorithm in a FORTRAN
  132.      version published by George Marsaglia and Arif Zaman, Florida State
  133.      
  134.      University; ref.: see original comments below.
  135.      
  136.      At the fhw (Fachhochschule Wiesbaden, W.Germany), Dept. of Computer
  137.      Science, we have written sources in further languages (C, Modula-2
  138.      Turbo-Pascal(3.0, 5.0), Basic and Ada) to get exactly the same test
  139.      results compared with the original FORTRAN version.
  140.      April 1989
  141.      Karl-L. Noell <NOELL@DWIFH1.BITNET>
  142.      and  Helmut  Weber <WEBER@DWIFH1.BITNET>
  143.      
  144.      This random number generator originally appeared in "Toward a Universal
  145.      Random Number Generator" by George Marsaglia and Arif Zaman.
  146.      Florida State University Report: FSU-SCRI-87-50 (1987)
  147.      It was later modified by F. James and published in "A Review of Pseudo-
  148.      random Number Generators"
  149.      
  150.      THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
  151.      
  152.      (However, a newly discovered technique can yield
  153.      a period of 64^600. But that is still in the development stage.)
  154.      It passes ALL of the tests for random number generators and has a period
  155.      of 2^144, is completely portable (gives bit identical results on all
  156.      machines with at least 24-bit mantissas in the floating point
  157.      representation).
  158.      
  159.      The algorithm is a combination of a Fibonacci sequence (with lags of 97
  160.      and 33, and operation "subtraction plus one, modulo one") and an
  161.      "arithmetic sequence" (using subtraction).
  162.      
  163.      
  164.      Use IJ = 1802 & KL = 9373 to test the random number generator. The
  165.      subroutine RANMAR should be used to generate 20000 random numbers.
  166.      Then display the next six random numbers generated multiplied by 4096*4096
  167.      If the random number generator is working properly, the random numbers
  168.      
  169.      should be:
  170.      
  171.      6533892.0  14220222.0  7275067.0
  172.      
  173.      6172232.0  8354498.0   64633180.0
  174.      
  175.      */
  176.      
  177.      
  178.     /* Globals */
  179.      
  180.     double u[97],c,cd,cm;
  181.      
  182.     int i97,j97;
  183.      
  184.     int test = FALSE;
  185.      
  186.      
  187.     /*
  188.      
  189.      This is the initialization routine for the random number generator.
  190.      NOTE: The seed variables can have values between:    0 <= IJ <= 31328
  191.      0 <= KL <= 30081
  192.      
  193.      The random number sequences created by these two seeds are of sufficient
  194.      length to complete an entire calculation with. For example, if sveral
  195.      different groups are working on different parts of the same calculation,
  196.      each group could be assigned its own IJ seed. This would leave each group
  197.      with 30000 choices for the second seed. That is to say, this random
  198.      number generator can create 900 million different subsequences -- with
  199.      each subsequence having a length of approximately 64^30.
  200.      
  201.      */
  202.      
  203.     void RandomInitialise(int ij,int kl)
  204.      
  205.     {
  206.            
  207.             double s,t;
  208.            
  209.             int ii,i,j,k,l,jj,m;
  210.      
  211.            
  212.      
  213.            
  214.            
  215.             /*
  216.              
  217.              Handle the seed range errors
  218.              
  219.              First random number seed must be between 0 and 31328
  220.              
  221.              Second seed must have a value between 0 and 30081
  222.              
  223.              */
  224.            
  225.             if (ij < 0 || ij > 31328 || kl < 0 || kl > 30081) {
  226.                    
  227.                     ij = 1802;
  228.                    
  229.                     kl = 9373;
  230.                    
  231.             }
  232.            
  233.            
  234.             i = (ij / 177) % 177 + 2;
  235.            
  236.             j = (ij % 177)       + 2;
  237.            
  238.             k = (kl / 169) % 178 + 1;
  239.            
  240.             l = (kl % 169);
  241.            
  242.            
  243.             for (ii=0; ii<97; ii++) {
  244.                    
  245.                     s = 0.0;
  246.                    
  247.                     t = 0.5;
  248.                    
  249.                     for (jj=0; jj<24; jj++) {
  250.                            
  251.                             m = (((i * j) % 179) * k) % 179;
  252.                            
  253.                             i = j;
  254.                            
  255.                             j = k;
  256.                            
  257.                             k = m;
  258.                            
  259.                             l = (53 * l + 1) % 169;
  260.                            
  261.                             if (((l * m % 64)) >= 32)
  262.                                    
  263.                                     s += t;
  264.                            
  265.                             t *= 0.5;
  266.                            
  267.                     }
  268.                    
  269.                     u[ii] = s;
  270.                    
  271.             }
  272.            
  273.            
  274.             c    = 362436.0 / 16777216.0;
  275.            
  276.             cd   = 7654321.0 / 16777216.0;
  277.            
  278.             cm   = 16777213.0 / 16777216.0;
  279.            
  280.             i97  = 97;
  281.            
  282.             j97  = 33;
  283.            
  284.             test = TRUE;
  285.            
  286.     }
  287.      
  288.      
  289.     /*
  290.      
  291.      This is the random number generator proposed by George Marsaglia in
  292.      
  293.      Florida State University Report: FSU-SCRI-87-50
  294.      
  295.      */
  296.      
  297.     double RandomUniform(void)
  298.      
  299.     {
  300.            
  301.             double uni;
  302.            
  303.            
  304.             /* Make sure the initialisation routine has been called */
  305.            
  306.             if (!test)
  307.                    
  308.                     RandomInitialise(1802,9373);
  309.            
  310.            
  311.             uni = u[i97-1] - u[j97-1];
  312.            
  313.             if (uni <= 0.0)
  314.                    
  315.                     uni++;
  316.            
  317.             u[i97-1] = uni;
  318.            
  319.             i97--;
  320.            
  321.             if (i97 == 0)
  322.                    
  323.                     i97 = 97;
  324.            
  325.             j97--;
  326.            
  327.             if (j97 == 0)
  328.                    
  329.                     j97 = 97;
  330.            
  331.             c -= cd;
  332.            
  333.             if (c < 0.0)
  334.                    
  335.                     c += cm;
  336.            
  337.             uni -= c;
  338.            
  339.             if (uni < 0.0)
  340.                    
  341.                     uni++;
  342.            
  343.            
  344.             return(uni);
  345.            
  346.     }
  347.      
  348.      
  349.     /*
  350.      
  351.      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
  352.      
  353.      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
  354.      
  355.      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
  356.      
  357.      The function returns a normally distributed pseudo-random number
  358.      
  359.      with a given mean and standard devaiation.  Calls are made to a
  360.      
  361.      function subprogram which must return independent random
  362.      
  363.      numbers uniform in the interval (0,1).
  364.      
  365.      The algorithm uses the ratio of uniforms method of A.J. Kinderman
  366.      
  367.      and J.F. Monahan augmented with quadratic bounding curves.
  368.      
  369.      */
  370.      
  371.     double RandomGaussian(double mean,double stddev)
  372.      
  373.     {
  374.            
  375.             double  q,u,v,x,y;
  376.            
  377.            
  378.             /*  
  379.              
  380.              Generate P = (u,v) uniform in rect. enclosing acceptance region
  381.              
  382.              Make sure that any random numbers <= 0 are rejected, since
  383.              
  384.              gaussian() requires uniforms > 0, but RandomUniform() delivers >= 0.
  385.              
  386.              */
  387.            
  388.             do {
  389.                    
  390.                     u = RandomUniform();
  391.                    
  392.                     v = RandomUniform();
  393.                    
  394.                     if (u <= 0.0 || v <= 0.0) {
  395.                            
  396.                             u = 1.0;
  397.                            
  398.                             v = 1.0;
  399.                            
  400.                     }
  401.                    
  402.                     v = 1.7156 * (v - 0.5);
  403.                    
  404.                    
  405.                     /*  Evaluate the quadratic form */
  406.                    
  407.                     x = u - 0.449871;
  408.                    
  409.                     y = fabs(v) + 0.386595;
  410.                    
  411.                     q = x * x + y * (0.19600 * y - 0.25472 * x);
  412.                    
  413.                    
  414.                     /* Accept P if inside inner ellipse */
  415.                    
  416.                     if (q < 0.27597)
  417.                            
  418.                             break;
  419.                    
  420.                    
  421.                     /*  Reject P if outside outer ellipse, or outside acceptance region */
  422.                    
  423.         } while ((q > 0.27846) || (v * v > -4.0 * log(u) * u * u));
  424.            
  425.            
  426.         /*  Return ratio of P's coordinates as the normal deviate */
  427.            
  428.         return (mean + stddev * v / u);
  429.            
  430.     }
  431.      
  432.      
  433.     /*
  434.      
  435.      Return random integer within a range, lower -> upper INCLUSIVE
  436.      
  437.      */
  438.      
  439.     int RandomInt(int lower,int upper)
  440.      
  441.     {
  442.            
  443.             return((int)(RandomUniform() * (upper - lower + 1)) + lower);
  444.            
  445.     }
  446.      
  447.      
  448.     /*
  449.      
  450.      Return random float within a range, lower -> upper
  451.      
  452.      */
  453.      
  454.     double RandomDouble(double lower,double upper)
  455.      
  456.     {
  457.            
  458.             return((upper - lower) * RandomUniform() + lower);
  459.            
  460.     }
clone this paste RAW Paste Data