1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
#include <stdio.h>
#include <errno.h>
#include <string.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
static int initialized = 0;
static double *x, *y;
static size_t size;
gsl_interp_accel *acc;
gsl_spline *spline;
static int init()
{
int i;
char line[256];
char *str;
double value;
int n;
FILE *f = fopen("pmt_pcath_response.dat", "r");
if (!f) {
fprintf(stderr, "failed to open pmt_pcath_response.dat: %s", strerror(errno));
return -1;
}
i = 0;
n = 0;
/* For the first pass, we just count how many values there are. */
while (fgets(line, sizeof(line), f)) {
size_t len = strlen(line);
if (len && (line[len-1] != '\n')) {
fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
goto err;
}
i += 1;
/* Skip the first 32 lines since it's just a header. */
if (i < 32) continue;
if (!len) continue;
else if (line[0] == '#') continue;
str = strtok(line," \n");
while (str) {
value = strtod(str, NULL);
n += 1;
str = strtok(NULL," \n");
}
}
x = malloc(sizeof(double)*(n+2));
y = malloc(sizeof(double)*(n+2));
size = n + 2;
/* Make sure we extrapolate to 0. */
x[0] = 0;
y[0] = 0;
x[n+1] = 1000;
y[n+1] = 0;
i = 0;
n = 0;
/* Now, we actually store the values. */
rewind(f);
while (fgets(line, sizeof(line), f)) {
size_t len = strlen(line);
if (len && (line[len-1] != '\n')) {
fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
goto err;
}
i += 1;
/* Skip the first 32 lines since it's just a header. */
if (i < 32) continue;
if (!len) continue;
else if (line[0] == '#') continue;
str = strtok(line," \n");
while (str) {
value = strtod(str, NULL);
/* According to the file, the values are stored for wavelengths
* between 230 and 700 in 1 nm increments. */
x[n+1] = 230 + n*1.0;
y[n+1] = value;
n += 1;
str = strtok(NULL," \n");
}
}
fclose(f);
acc = gsl_interp_accel_alloc();
spline = gsl_spline_alloc(gsl_interp_linear, size);
gsl_spline_init(spline, x, y, size);
initialized = 1;
return 0;
err:
fclose(f);
return -1;
}
double get_quantum_efficiency(double wavelength)
{
if (!initialized) {
if (init()) {
exit(1);
}
}
return gsl_spline_eval(spline, wavelength, acc);
}
|
