1 |
chrisfen |
1845 |
program convolve |
2 |
|
|
character*80 junk, option |
3 |
|
|
integer nargs, iargc, npts, maxpts, i, j |
4 |
|
|
parameter (maxpts=10000) |
5 |
|
|
double precision width, x(maxpts), y(maxpts), myX, myY, dx, convol |
6 |
|
|
double precision instrument_gaussian, deltaX |
7 |
|
|
|
8 |
|
|
external iargc, instrument_gaussian |
9 |
|
|
nargs = iargc() |
10 |
|
|
|
11 |
|
|
if (nargs.ne.2) then |
12 |
|
|
write(0,*) 'This program is a filter for convoluting a data file by' |
13 |
|
|
write(0,*) 'a gaussian instrument function of a particular width.' |
14 |
|
|
write(0,*) |
15 |
|
|
write(0,*) 'Sample usage:' |
16 |
|
|
write(0,*) ' convolve -w 0.5 < file.dat > new.dat' |
17 |
|
|
write(0,*) |
18 |
|
|
write(0,*) 'The -w flag tells the program that the argument that' |
19 |
|
|
write(0,*) 'follows is the width of the gaussian instrument function' |
20 |
|
|
write(0,*) 'in the same units as the x column in the data file.' |
21 |
|
|
stop |
22 |
|
|
endif |
23 |
|
|
|
24 |
|
|
call getarg(1,junk) |
25 |
|
|
read(junk,*) option |
26 |
|
|
if (option.ne.'-w') then |
27 |
|
|
write(0,*) 'Your first argument is not a -w! What do you want to do?' |
28 |
|
|
stop |
29 |
|
|
endif |
30 |
|
|
call getarg(2,junk) |
31 |
|
|
read(junk,*) width |
32 |
|
|
|
33 |
|
|
npts = 1 |
34 |
|
|
do while (.true.) |
35 |
|
|
read(5,*,end=40) myX, myY |
36 |
|
|
x(npts) = myX |
37 |
|
|
y(npts) = myY |
38 |
|
|
npts = npts + 1 |
39 |
|
|
end do |
40 |
|
|
|
41 |
|
|
40 continue |
42 |
|
|
|
43 |
|
|
deltaX = x(2) - x(1) |
44 |
|
|
|
45 |
|
|
do i = 1, npts-1 |
46 |
|
|
convol = 0.0d0 |
47 |
|
|
do j = 1, npts-1 |
48 |
|
|
dx = x(j) - x(i) |
49 |
|
|
convol = convol + y(j)*instrument_gaussian(dx, width) * deltaX |
50 |
|
|
enddo |
51 |
|
|
convol = convol |
52 |
|
|
write(*,*) x(i), convol |
53 |
|
|
enddo |
54 |
|
|
|
55 |
|
|
end program convolve |
56 |
|
|
|
57 |
|
|
double precision function instrument_gaussian(x, sigma) |
58 |
|
|
|
59 |
|
|
double precision pi, sigma, x |
60 |
|
|
|
61 |
|
|
pi = 4.0d0 * datan(1.0d0) |
62 |
|
|
instrument_gaussian = exp(- x*x / (2.0d0 * sigma * sigma)) / (sigma * sqrt(2.0d0 * pi)) |
63 |
|
|
|
64 |
|
|
return |
65 |
|
|
end function instrument_gaussian |
66 |
|
|
|