Help with ran1 and gasdev - generating sets of standard deviates

This is a discussion on Help with ran1 and gasdev - generating sets of standard deviates within the C Programming forums, part of the General Programming Boards category; Hello! I am having trouble with a project I'm working on for my professor. I have to generate a few ...

  1. #1
    Registered User
    Join Date
    Jun 2011
    Posts
    2

    Help with ran1 and gasdev - generating sets of standard deviates

    Hello! I am having trouble with a project I'm working on for my professor. I have to generate a few sets of standard deviations for a physics project. My professor gave me the code for ran1 and gasdev from numerical recipes in C. They generate standard deviations fine, but every time I run the program, I get the exact same sets of standard deviations. My professor said to set the seed equal to something that changes, like the number of seconds since midnight. I'm not sure how to do that, and I'm new to programming. Could you help, please? Here's the whole program:

    Code:
    /*set20.c*/
    #include<stdio.h>
    #include<math.h>
    #include<stdlib.h>
    #include<time.h>
    #include<string.h>
    #define MAXLINE 400
    
    /*ran1.c from Numerical Recipes - This outputs random numbers to gasdev.c*/
    #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 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;
    	}
    
    	#undef IA
    	#undef IM
    	#undef AM
    	#undef IQ
    	#undef IR
    	#undef NTAB
    	#undef NDIV
    	#undef EPS
    	#undef RNMX						/* (C) Copr. 1986-92 Numerical Recipes Software i9k''3. */
    
    	
    
    /*gasdev.c - This returns pseudorandom numebers along a normal distribution.*/
    float gasdev(long *idum)
    	{
    	float ran1(long *idum);
    	static int iset=0;
    	static float gset;
    	float fac,rsq,v1,v2;
    	
    	if  (iset == 0) {
    		do 
    			{
    			v1=2.0*ran1(idum)-1.0;
    			v2=2.0*ran1(idum)-1.0;
    			rsq=v1*v1+v2*v2;
    			}
    		while (rsq >= 1.0 || rsq == 0.0);
    		fac=sqrt(-2.0*log(rsq)/rsq);
    		gset=v1*fac;
    		iset=1;
    		return v2*fac;
    		}
    	else 	{
    		iset=0;
    		return gset;
    		}
    	}							/* (C) Copr. 1986-92 Numerical Recipes Software i9k''3. */
    /*This part outputs the number of lines in the file.*/
    int linecount(FILE *fp) {
    	  char buff[MAXLINE];
    	  int count = 0;
    	  
    	  while(fgets(buff,MAXLINE,fp) != NULL)
    	    {
    	      count++;
    	    }
    	  return count;
    	}
    
    /*Main function*/
    int main ( int argc, char *argv[] )
    	{
    	int N;							/*Number of Pulsars*/
    	float P0m;						/*Initial Period Mean*/
    	float P0s;						/*Initial Period std. dev.*/
    	float B0m;						/*Initial Mag. Field Mean*/
    	float B0s;						/*Initial Mag. Field std. dev.*/
    	
    	
    /*This section Parses the command line arguments.*/		
        	int m, n;                              			/* Loop counters. */
            int l;                               			/* String length. */
            int x;                             		    	/* Exit code. */
            int ch;                               			/* Character buffer. */
        	char q[256];                           			/* String buffer. */
            float agesinyears[1000];
    	int FILETHERE;
    
    	if (argv[1]==0)
    	{
    	N=0;
    	printf("Set20.c\nGenerates Simulated Set of Pulsars.\nUSAGE:\n-n<# of Pulsars>\n-P<Initial Period mean (log10)>\n-p<Initial Period stdev.(log10)\n-B<Magnetic Field mean(log10)>\n-b<Magnetic field stdev.(log10)\n-f<age list file(optional)>\n-h Help\n");
    	}
    
    	for( n = 1; n < argc; n++ )            			/* Scan through args. */
    	     		{
    	      		switch( (int)argv[n][0] )            		/* Check for option character. */
    	      			{
    	       			 case '-':
    	       			case '/': x = 0;                   	/* Bail out if 1. */
    		         	l = strlen( argv[n] );
    		         	for( m = 1; m < l; ++m ) 		/* Scan through options. */
    		         		{
    		           		ch = (int)argv[n][m];
    		           		switch( ch )
    		           			{
    			   			case 'n': 	if( m + 1 >= l )
    						     			{
    						    			exit( 1 );
    									N=100;
    						    			}
    							     	else
    							     		{
    							       		strcpy( q, &argv[n][m+1] );
    									N=atoi(q);
    							    		}
    							     	x = 1;
    		                     		break;
    			  			case 'P': 	if( m + 1 >= l )
    								     	{
    									exit( 1 );
    								  	}
    		                     				else
    							     		{
    							       		strcpy( q, &argv[n][m+1] );
    									P0m=atof(q);
    							     		}
    							     	x = 1;
    		                     		break;
    					   	case 'p': 	if( m + 1 >= l )
    								     	{
    								       	exit( 1 );
    								     	}
    		                     				else
    							     		{
    							       		strcpy( q, &argv[n][m+1] );
    									P0s=atof(q);
    							     		}
    		                     				x = 1;
    		                     		break;
    			   			case 'B': 	if( m + 1 >= l )
    								     	{
    								       	exit( 1 );
    								     	}
    		                     				else
    							     		{
    							       		strcpy( q, &argv[n][m+1] );
    									B0m=atof(q);
    							     		}
    		                     				x = 1;
    		                     				break;
    			   			case 'b': 	if( m + 1 >= l )
    								     	{
    								       	exit( 1 );
    								     	}
    		                     				else
    								     	{
    								       	strcpy( q, &argv[n][m+1] );
    									B0s=atof(q);
    								     	}
    				        			x = 1;
    		                     				break;
    			   			case 'f': 	if( m + 1 >= l )
    								     	{
    								       	exit( 1 );
    									FILETHERE=0;
    								     	}
    		                     				else
    								     	{
    								       	strcpy( q, &argv[n][m+1] );
    									FILETHERE=1;		
    				
    									/*Stores file entries into an array*/
    				
    									  FILE *file = fopen( q, "r" );
    									  if ( file == 0 )
    									    {
    									      printf( "Can't open file!\n" );
    									      exit(0);
    									    }
    
    									  FILE *fp;
    									  int Na = 0;
    									  if ((fp=fopen(q,"r")) != NULL) {
    									    Na = linecount(fp);
    									  }
    									  fclose(fp);
    									  
    									  char line[MAXLINE];
    									  double agey;
    									  int ia=0;
    									  while  ( ( fgets( line, MAXLINE, file ) ) != NULL )
    										    {
    										    ia++;
    										    agey = atof(line);
    										    agesinyears[ia]= agey;
    										    }
    									fclose( file );
    									}  	
    				       	 			x = 1;
    		                     				break;
    			   			case 'h': 	N=0;
    								printf("Set20.c\nGenerates Simulated Set of Pulsars.\nUSAGE:\n-n<# of Pulsars>\n-P<Initial Period mean (log10)>\n-p<Initial Period stdev.(log10)\n-B<Magnetic Field mean(log10)>\n-b<Magnetic field stdev.(log10)\n-f<age list file(optional)>\n-h Help\n");
    								
    				        			x = 1;
    		                     				break;
    						default:  	
    							     x = 1;      /* Not legal option. */
    							     exit( 1 );
    						break;
    						}
    		           	if( x == 1 )
    		           		{
    		             		break;
    		           		}
    		         	}
    		         	break;
    	       			default:   /* Not option -- text. */
    		        break;
    	       		}
    	     	}
    
    
    /*Computes the X-Ray Luminosity*/
    	int i;							/*Individual pulsar number.  This counts up from one.*/
    	i=1;	
    	long double P0;						/*Initial period (s)*/
    	long double B0;						/*Initial mag. field (G)*/
    	long double t;						/*Age (s)*/
    	long double I;						/*Moment of Inertia (g*cm^2)*/
    	I=1e45;	
    	long double R;						/*Radius (cm)*/
    	R=1e6;
    	long double y;						/*Age of the Pulsar in years(for testing purposes)*/
    	float b;						/*Log-base-10 of initial magnetic field strength*/
    	float p;						/*Log-base-10 of initial magnetic field strength*/
    	float r;						/*Normal distribution output used for B0 calculation*/
    	float s;						/*Normal distribution output used for P0 calculation*/
    	long double pi;						/*Declares the values of pi and c.  c is in cgs units.*/
    	pi=3.1415926536;
    	long double c;					
    	c=3e10;
    	long utime;
    	
    	utime=(long)time(NULL);					/*Unix time*/
    	long seed=utime;						/*Random seed for P0 and B0 distributions*/					
    	
    	float bsigma;						/*Mean and standard deviation for log-10 of magnetic field*/
    	float bmu;
    	bsigma=B0s;
    	bmu=B0m;
    	float psigma;						/*Mean and standard deviation for log-10 of initial period.  From Perna et al 2008.*/
    	float pmu;
    	psigma=P0s;
    	pmu=P0m;
    	long double k;						/*Constant -  16pi^2R^6B^2/3Ic^3*/						
    	long double Pt;						/*Time-dependent Period*/	
    	long double Pdot;					/*Time derivative of period*/	 	
    	long double o;						/*Omega, the angular velocity at time t*/		
    	long double odot;					/*Time derivative of o*/
    	long double E;						/*Rotational Kinetic Energy*/				
    	long double Edot;					/*Time derivative of Rotational Kinetic Energy*/
    	long double Lx;						/*X-Ray Luminosity in erg/s*/	
    	long double muLx;					/*Log-base-10 of Lx*/
    	float sigmaLx;
    	float sigmaLxa;
    	float sigmaLxb;
    	float lEdot;
    	float lLx;
    	
    	srand ( time(NULL) );					/*Initializes random seed for t:*/
    	do	
    		{						/*Loop runs until i=N*/	
    		if(FILETHERE==0)
    			{
    			t=random() % 2428272000u;
    			}
    		if(FILETHERE==1)
    			{
    			t=agesinyears[i]*31536000u;
    			}
    		b=(gasdev(&seed)*bsigma)+bmu;			/*B0 and P0 along normal distribution.*/
    		B0=pow(10.0, b);				
    		p=(gasdev(&seed)*psigma)+pmu;
    		P0=pow(10.0, p);
    
    		y=t/31536000u;					/*Age in years (for testing purposes)*/			
    		/*X-Ray Luminosity calculation*/
    		k=((16*pi*pi*R*R*R*R*R*R*B0*B0) /( 3*I*c*c*c));	/*Constant*/
    		Pt= sqrt((P0*P0)+(k*t));			/*Time-dependent Period Function*/
    		Pdot= 0.5*k/sqrt(P0*P0+k*t);			/*Time derivative of period*/
    		o=(2*pi)/Pt;					/*Angular velocity*/
    		odot=(-2*pi*Pdot)/(Pt*Pt);			/*Time derivative of o*/
    		E=(I*o*o)/2;					/*Rotational kinetic energy*/
    		/*Edot=I*o*odot;*/
    		Edot=-(B0*B0*o*o*o*o*R*R*R*R*R*R)/(6*c*c*c);	/*Time derivative of E*/
    		Lx=pow(10, (1.34*log10(-Edot)-15.34));		/*X-Ray Luminosity function*/
    		muLx = log10(Lx);
    
    
    		sigmaLxa=0.03;
    		sigmaLxb=1.11;
    			
    		lEdot= log10(-Edot);		
    		sigmaLx = sqrt(((pow(sigmaLxa,2))*(pow(lEdot,2)))+(pow(sigmaLxb,2)));
    		
    		lLx=(gasdev(&seed)*sigmaLx)+muLx;
    
    		/*
    		printf("%Lf\n",muLx);
    		printf("%f\n",sigmaLx);
    		printf("%f\n",sigmaLxa);
    		printf("%f\n",sigmaLxb);
    		printf("%f\n",lEdot);		
    		*/
    		
    		/*Prints*/
    		printf("%g\n", lLx ); 
    		i++;						/*Updates pulsar count*/
    		}while(i<N+1);					/*Loop conditions*/	
    	
    	return 0;	
    
    					
    	}							/*End of Main*/

    Thanks for reading!

  2. #2
    Registered User
    Join Date
    May 2011
    Location
    Around 8.3 light-minutes from the Sun
    Posts
    1,866
    Ok, this is where you set your seed in your code:
    Code:
    long utime;
    utime=(long)time(NULL);					/*Unix time*/
    long seed=utime;
    To return the num seconds since midnight....of Jan 1 1970:
    Code:
    time_t utime;
    time(&utime);
    long seed = (long)utime;
    Quote Originally Posted by anduril462 View Post
    Now, please, for the love of all things good and holy, think about what you're doing! Don't just run around willy-nilly, coding like a drunk two-year-old....
    Quote Originally Posted by quzah View Post
    ..... Just don't be surprised when I say you aren't using standard C anymore, and as such,are off in your own little universe that I will completely disregard.
    Warning: Some or all of my posted code may be non-standard and as such should not be used and in no case looked at.

  3. #3
    Registered User
    Join Date
    Jun 2011
    Posts
    2
    I just did that, and as it turns out I also had to make the seed negative in order for it to work. Thanks a billion!

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Map of sets?
    By Orange Lozenge in forum C++ Programming
    Replies: 6
    Last Post: 01-23-2011, 10:16 AM
  2. Replies: 3
    Last Post: 01-22-2011, 12:22 PM
  3. Sets Problem
    By inevitable in forum C++ Programming
    Replies: 5
    Last Post: 04-03-2008, 02:45 PM
  4. ATD: C++ Sets
    By stalker in forum C++ Programming
    Replies: 9
    Last Post: 11-09-2004, 05:40 AM
  5. standard language, standard compiler?
    By doubleanti in forum A Brief History of Cprogramming.com
    Replies: 8
    Last Post: 09-03-2001, 04:21 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21