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
| /* input_matrix() reads a gnuplot-compatible binary file, i.e. single precision *
* matrix stored in the following format: *
* *
* <N+1> <y0> <y1> <y2> ... <yN> *
* <x0> <z0,0> <z0,1> <z0,2> ... <z0,N> *
* <x1> <z1,0> <z1,1> <z1,2> ... <z1,N> *
* *
* For example: input_matrix(T,fp,N,X0,Y0,L0); *
* Reads a square matrix of size N from file descriptor "fp" into scalar field *
* T, defined at points (X,Y). *
*/
struct InputMatrix {
// compulsory
scalar s;
FILE * fp;
// optional
int n;
double ox, oy, width, oz;
double arz;
};
void input_matrix (struct InputMatrix p) {
scalar s = p.s;
if (p.width == 0.) p.width = L0;
float width=0;
fread(&width, sizeof(float), 1, p.fp);
if (p.n != width) p.n = (int) width;
float yp[p.n], xp[p.n], v[p.n][p.n];
fread(&yp, sizeof(float), p.n, p.fp);
for (int i = 0; i < p.n; i++) {
fread(&xp[i], sizeof(float), 1, p.fp);
for (int j = 0; j < p.n; j++) {
fread(&v[i][j], sizeof(float), 1, p.fp);
}
}
foreach() {
int i = (x - p.ox)*width/p.width, j = (y - p.oy)*width/p.width;
if (i >= 0 && i < width && j >= 0 && j < width)
s[] = v[i][j];
else
s[] = 0.;
}
}
|
