samedi 5 septembre 2020

Plotting random numbers in a historgram in C language

I use a code from "Numerical Recipes in C book" to generate a uniform random number which is named as function float ran1(long *idum). It produces a uniform random number between 0 and 1. I am trying to bin these uniform random numbers and make a histogram. I wrote the below code, however when I plot the bins versus the frequency, I do not get the uniform distribution.

Can someone help me know what the problem is?

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)


float ran1(long *idum);

int main(){
    
    double delta, x;
    int i , ibin ;
    long idum ;
    int histog[500] ; 
    
    
    delta=0.02 ;
    
    idum=-1 ;// necessary to generate the uniform random number, must be negative integer
    
    FILE *file1;
    
    file1=fopen("guassians.dat","w") ;
    
    
    for (i=0 ; i<=10000 ; ++i) {
        
        x=ran1(&idum) ; // here ran1(&idum) returns a uniform random number
      
        
        ibin=round(x/delta);
        
        if ( abs(ibin) < 250 ) {
            
            histog[ibin]=histog[ibin] + 1 ;
        }
    }
    
    for (ibin=-250 ; ibin <= 250 ; ++ibin){
        
        fprintf(file1, "%0.10lf \t %d \n", ibin*delta , histog[ibin]);
        
    }
    

        
    fclose(file1);
        
}



float ran1(long *idum){
    
    int j ;
    long k;
    static long iy=0;
    static long iv[NTAB];
    float temp ;
    
    if (*idum <= 0 || !iy) {
        
        if (-(*idum) < 1) *idum=1 ;
        else *idum = -(*idum) ;
        for (j=NTAB+7 ; j>=0 ; j--) {
            k=(*idum)/IQ ;
            *idum=IA*(*idum-k*IQ)-IR*k ;
            if (*idum < 0) *idum += IM;
            if (j < NTAB) iv[j] = *idum;
        }
        iy=iv[0] ;
    }
    k=(*idum)/IQ ;
    *idum=IA*(*idum-k*IQ)-IR*k ;
    if (*idum < 0) *idum += IM ;
    j=iy/NDIV ;
    iy=iv[j] ;
    iv[j] = *idum ;
    if ((temp=AM*iy) > RNMX) return RNMX ;
    else return temp;
        
}


enter image description here




Aucun commentaire:

Enregistrer un commentaire