!-----------------------------------------------------------------------
!Main program
!-----------------------------------------------------------------------
program test
![specification part]
implicit none
integer, parameter :: rsp = kind(1.0)
integer, parameter :: rdp = selected_real_kind(2*precision(1.0_rsp))
real(rdp), allocatable, dimension(:) :: x,y
real(rdp) :: width
integer :: pulse_type, polarity
integer :: i, mult,n
real(rdp), parameter :: pi=3.14159265358979323846264338327950288419716939937510_rdp
real(rdp):: a,b,c,xi
open(20,file="sin.dat", status="unknown", action="write")
n=100000
mult=3.0_rdp
do i=1,n
xi=(mult*pi)*dble(i)/dble(n)
write(20,'(2E24.16)') xi,sin(xi)
enddo
close(20)
write(*,*) "exact width for sine=", 4.0_rdp*pi/6.0_rdp
call init_pulse ("sin.dat",x,y )
call fwhm(x,y, width, pulse_type, polarity)
write(*,*) "width=", width
write(*,*) "pulse_type=",pulse_type
write(*,*) "polarity=",polarity
open(20,file="gaussian.dat", status="unknown", action="write")
n=100000
mult=3.0_rdp
a=1.0_rdp/sqrt(2.0_rdp*pi)
b=0.d0
c=1
do i=-n,n
xi=mult*dble(i)/dble(n)
write(20,'(2E24.16)') xi ,a*exp(-((xi-b)**2)/(2.0_rdp*c**2))
enddo
close(20)
write(*,*) "exact width for gaussian=", 2.0_rdp*sqrt(2.0_rdp*log(2.0_rdp))*c
call init_pulse ("gaussian.dat",x,y )
call fwhm(x,y, width, pulse_type, polarity)
write(*,*) "width=", width
write(*,*) "pulse_type=",pulse_type
write(*,*) "polarity=",polarity
contains
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
subroutine init_pulse (fname,x,y )
implicit none
character(len=*)::fname
real(rdp), allocatable, dimension(:) :: x,y
integer, parameter :: find=20
integer :: nl,k,ios
real(rdp) ::a,b
open(find,file=fname, status="unknown", action="read")
nl=0 !number of lines
!detect number of lines in screen.dat
do
read(find,*, iostat=ios) a,b
if (ios /= 0) then
write(*,*) "Detected number of lines in file =",nl
exit
else
nl=nl+1
endif
end do
!go back to start
rewind(find,iostat=ios)
if (ios/=0) then !something is wrong
write(*,*) "cannot rewind file. STOP"
stop
endif
if (allocated(x)) deallocate(x)
if (allocated(y)) deallocate(y)
allocate(x(1:nl))
allocate(y(1:nl))
do k=1,nl
read(find,*) x(k),y(k)
enddo
close(find)
end subroutine init_pulse
!-----------------------------------------------------------------------
! full-width at half-maximum (fwhm) of the waveform y(x)
! and its polarity.
!
! polarity=+-1 polarity of the pulse
! pulse_type=+1 pulse is impulse or rectangular with 2 edges
! pulse_type=0 step-like pulse, no second edge
!-----------------------------------------------------------------------
subroutine fwhm(x,y, width, pulse_type, polarity)
implicit none
real(rdp), dimension(:), intent(in):: x,y
real(rdp) , intent(out):: width
integer , intent(out):: pulse_type, polarity
real(rdp) :: maxval_y,minval_y,interp,t_lead,t_trail
integer :: maxloc_y(1),minloc_y(1)
real(rdp), parameter :: lev_half=0.5_rdp
integer :: n,c_index,i
!initial bad values
polarity=0
pulse_type=0
width=0.0_rdp
! find index of center (max or min) of pulse
maxval_y=maxval(y)
maxloc_y=maxloc(y)
minval_y=minval(y)
minloc_y=minloc(y)
n=size(y)
!simple checks
if(n<3) then
write(*,*) "Warning from fwhm: dimension of signal <",n
return
endif
if(abs(maxval_y)<tiny(1.0_rdp) .and. abs(minval_y)<tiny(1.0_rdp)) then
write(*,*) "Warning from fwhm: no signal"
return
endif
i=1
if(y(i)/maxval_y<lev_half) then
c_index=maxloc_y(i)
polarity=1
else
c_index=minloc_y(i)
polarity=-1
endif
!position from the left is between y(i-1) & y(i)
do i=2,n-1
if (sign(1.0_rdp,y(i)/maxval_y-lev_half)/= sign(1.0_rdp,y(i-1)/maxval_y-lev_half)) then
write(*,*)" i in exit=",i
exit
endif
if(i==n-1) then
write(*,*) "Warning from fwhm: problem in signal, no edge"
return !problem, no edge
endif
end do
interp=(lev_half-y(i-1)/maxval_y)/(y(i)/maxval_y-y(i-1)/maxval_y)
t_lead=x(i-1)+interp*(x(i)-x(i-1))
! search for the next position starting from center of the pulse
do i=c_index+1,n-1
if (sign(1.0_rdp,y(i)/maxval_y-lev_half)/= sign(1.0_rdp,y(i-1)/maxval_y-lev_half)) then
exit
endif
if(i==n-1) then
write(*,*) "Warning from fwhm: problem in signal, step-like pulse, no second edge"
return !problem, step-like pulse, no second edge
endif
end do
!ok pulse is impulse or rectangular with 2 edges
pulse_type=1
interp=(lev_half-y(i-1)/maxval_y)/(y(i)/maxval_y-y(i-1)/maxval_y)
t_trail=x(i-1)+interp*(x(i)-x(i-1))
width=t_trail-t_lead
end subroutine fwhm
end program test
|
