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
|
#include "Main.h"
#include <stdio.h>
#include <stdlib.h>
#include "Defines.h"
#include "Error.h"
#include "HpFile.h"
#include "Utilities.h"
/* own stuff */
#include "AreaBelow.h"
/*
* Return the area enclosed by all of the curves. The algorithm
* used is the same as the trapizoidal rule for integration.
*/
floatish
AreaBelow(void)
{
intish i;
intish j;
intish bucket;
floatish value;
struct chunk *ch;
floatish area;
floatish trap;
floatish base;
floatish *maxima;
maxima = (floatish *) xmalloc(nsamples * sizeof(floatish));
for (i = 0; i < nsamples; i++) {
maxima[i] = 0.0;
}
for (i = 0; i < nidents; i++) {
for (ch = identtable[i]->chk; ch; ch = ch->next) {
for (j = 0; j < ch->nd; j++) {
bucket = ch->d[j].bucket;
value = ch->d[j].value;
if (bucket >= nsamples)
Disaster("bucket out of range");
maxima[ bucket ] += value;
}
}
}
area = 0.0;
for (i = 1; i < nsamples; i++) {
base = samplemap[i] - samplemap[i-1];
if (maxima[i] > maxima[i-1]) {
trap = base * maxima[i-1] + ((base * (maxima[i] - maxima[i-1]))/ 2.0);
} else {
trap = base * maxima[i] + ((base * (maxima[i-1] - maxima[i]))/ 2.0);
}
area += trap;
}
free(maxima);
return area;
}
|
