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!
Originally Posted by anduril462
