40 use iso_fortran_env,
only: real64
44 real(real64),
parameter :: sigma = 1e-3
45 integer,
parameter :: niter = 5
47 real(real64),
allocatable :: dips(:, :)
48 character(len=1024) :: rawname, smoothname
49 integer :: iarg, narg, dot
51 narg = command_argument_count()
54 print
'(A)',
'Just give me a list of files to smooth...'
61 call get_command_argument(iarg, rawname)
65 dot = scan(rawname,
'.', back = .true.)
67 smoothname = smoothname(1:dot-1) //
'-smooth' // smoothname(dot:)
70 print
'(3A)', trim(rawname),
' -> ', trim(smoothname)
83 character(len=*),
intent(in) :: filename
84 real(real64),
allocatable,
intent(inout) :: dips(:, :)
86 integer :: i, u, ierr, n
89 open (newunit = u, file = filename, action =
'read', form =
'formatted')
95 read (u, *, iostat = ierr) x, y
102 if (
allocated(dips))
then
103 if (
size(dips, 2) /= n)
then
107 if (.not.
allocated(dips))
then
108 allocate (dips(2, n))
114 read (u, *) dips(:, i)
124 character(len=*),
intent(in) :: filename
125 real(real64),
allocatable,
intent(in) :: dips(:, :)
129 open (newunit = u, file = filename, action =
'write', form =
'formatted')
131 write (u,
'(2E25.15E3)') dips
140 real(real64),
allocatable,
intent(inout) :: dips(:, :)
141 real(real64),
allocatable :: weights(:), dips0(:, :), dips1(:, :)
142 real(real64) :: weight, weight_sum, re_d, im_d
143 integer :: i, j, iter, n0, n, d
147 allocate (weights(n), dips0(2, n), dips1(2, n))
151 weights(i) = exp(-sigma*(i - 1)**2)
158 if (any(dips(:, i) /= 0))
then
174 re_d = dips(1, j) - dips0(1, j)
175 im_d = dips(2, j) - dips0(2, j)
176 weight = weights(d) / sqrt(1 + re_d*re_d + im_d*im_d)
177 dips1(:, i) = dips1(:, i) + dips(:, j) * weight
178 weight_sum = weight_sum + weight
180 dips1(:, i) = dips1(:, i) / weight_sum
subroutine read_dipoles(filename, dips)
program multidip_smooth_rawdips
Smoothing of MULTIDIP raw transition dipoles.
subroutine smooth_dipoles(dips)
subroutine write_dipoles(filename, dips)