| 1 |
gezelter |
1468 |
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 |
|
|
|
