StatMech
Loading...
Searching...
No Matches
maximum.c File Reference

Functions

int main (int argc, char **argv)
 

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)
16 {
17 if (argc < 3) {
18 printf("Usage: 1.This 2.n 3.ndata\n");
19 exit(EXIT_FAILURE);
20 }
21 Integer_t i, j, count;
22 Integer_t const n = atoi(argv[1]);
23 Integer_t const ndata = atoi(argv[2]);
24 Integer_t seed=0;
25 double max, maxData;
26 double rand, rand1, rand2;
27
28 Integer_t subMax;
29 Integer_t *histogram;
30 double mean, stddev;
31 double binwidth=0.05;
32
33 char str[100];
34 char *filename = str;
35
36 FILE *Infp;
37 FILE *Outfp;
38 sprintf(str, "Dist_maximum_%lld_raw.txt", n);
39 Outfp = fopen(filename, "w");
40 if (Outfp == NULL) {
41 printf("Can't open a file %s\n", filename);
42 exit(EXIT_FAILURE);
43 }
44 fprintf(Outfp, "#(n=%lld, seed=%lld) Distribution of maximum of |x| \n\n", n, seed);
45
46 init_genrand(seed);
47 maxData = NAN;
48 mean = 0;
49 stddev = 0;
50 for(i = 0;i < ndata; ++i) {
51 max = NAN;
52 for(j = 0;j < n; ++j) {
53 rand = genrand_real3();
54 max = fmax(max, rand);
55 }
56 maxData = fmax(maxData, max);
57 mean += max;
58 stddev += max*max;
59 fprintf(Outfp, "%lld %le\n", i, max);
60 }
61 mean /= (double)ndata;
62 stddev /= (double)ndata;
63 stddev = sqrt( stddev -mean*mean );
64 binwidth *= stddev;
65 subMax = (Integer_t)(maxData/binwidth) +1;
66 histogram = alloc_ivector(subMax); for(i = 0;i < subMax; ++i) histogram[i] = 0;
67
68 sprintf(str, "Dist_maximum_%lld_raw.txt", n);
69 Infp = fopen(filename, "r");
70 if (Infp == NULL) {
71 printf("Can't open a file %s\n", filename);
72 exit(EXIT_FAILURE);
73 }
74 skip_header(Infp);
75
76 count=0;
77 while ( fscanf(Infp, "%lld %le", &i, &max) != EOF ) {
78 histogram[ (Integer_t)(max/binwidth) ] += 1;
79 count += 1;
80 }
81
82 sprintf(str, "Dist_maximum_%lld.txt", n);
83 Outfp = fopen(filename, "w");
84 if (Outfp == NULL) {
85 printf("Can't open a file %s\n", filename);
86 exit(EXIT_FAILURE);
87 }
88 fprintf(Outfp, "# mean=%le, stddev=%le\n\n", mean, stddev);
89 fprintf(Outfp, "# 1.Maximum value 2.Probability density\n\n");
90 for(i = 0;i < subMax; ++i) {
91 fprintf(Outfp, "%le %le\n", binwidth*(i+0.5), histogram[i]/(binwidth*count));
92 }
93
94 return 0;
95}
MKL_INT Integer_t
Definition mytypes.hpp:359