program gauss

implicit real*8 (a - h, o - z)
parameter (maxkpno = 20001)
dimension pos (maxkpno), spec (maxkpno), specmod (maxkpno), slit (maxkpno)
character*40 input, output
write (*, '(5x, a)') 'enter filespec of input file.'
read (*, '(a)') input
write (*, '(5x, a)') 'enter filespec of output file.'
read (*, '(a)') output
open (unit = 1, file = input, status = 'old')
open (unit = 2, file = output, status = 'unknown')

write (*, '(5x, a)') 'enter hw1/e.'
read (*, *) hw1e

i = 1
20 read (1, *, end = 30) pos (i), spec (i)
  i = i + 1
  go to 20
30 npoints = i - 1
write (*, *) ' npoints = ', npoints

call gaussian (pos, spec, specmod, npoints, hw1e)
do i = 1, npoints
  write (2, '(f10.5, 1pe13.5)') pos (i), specmod (i)
end do

close (unit = 1)
close (unit = 2)
stop
end

subroutine gaussian (pos, spec, specmod, npoints, hw1e)

! convolves input spectrum with gaussian slit function of specified hw1e
! (half-width at 1/e intensity). Assumes input spectrum has constant spacing.

implicit real*8 (a - h, o - z)
parameter (maxpts = 20001)
dimension pos (maxpts), spec (maxpts)
dimension specmod (maxpts), slit (maxpts)

if (hw1e .eq. 0) then
  write (*, *) ' no gaussian convolution applied.'
  do i = 1, npoints
    specmod (i) = spec (i)
  end do
  return
end if

emult = - 1.d0 / (hw1e**2)
delpos = pos (2) - pos (1)

! apply slit function convolution.

! calculate slit function values out to 0.001 times x0 value, normalize so that
! sum = 1.
slitsum = 1.d0
slit0 = 1.d0
do i = 1, 1000
  slit (i) = exp (emult * (delpos * i)**2)
  slitsum = slitsum + 2.d0 * slit (i)
  if (slit (i) / slit0 .le. 0.001) then
    nslit = i
    write (*, *) ' nslit = ', nslit
    go to 40
  end if
end do
40 continue
slit0 = slit0 / slitsum
do i = 1, nslit
  slit (i) = slit (i) / slitsum
end do

! convolve spectrum. don't reflect at endpoints.
do i = 1, npoints
  specmod (i) = slit0 * spec (i)
  do j = 1, nslit
    nlo = i - j
    nhi = i + j
    if (nlo .ge. 1) specmod (i) = specmod (i) + slit (j) * spec (nlo)
    if (nhi .le. npoints) specmod (i) = specmod (i) + slit (j) * spec (nhi)
  end do
end do

! Convolve spectrum.  reflect at endpoints.
! do i = 1, npoints
!   specmod (i) = slit0 * spec (i)
!   do j = 1, nslit
!     nlo = i - j
!     nhi = i + j
!     if (nlo .lt. 0) nlo = - nlo
!     if (nhi .gt. npoints) nhi = 2 * npoints - nhi
!     if (nlo .ne. 0) then
!       specmod (i) = specmod (i) + slit (j) * (spec (nlo) + spec (nhi))
!     else
!       specmod (i) = specmod (i) + slit (j) * spec (nhi)
!     end if
!   end do
! end do

return
end