* This program takes a DIFFaX ".spc" file (or a suite of files)
* and re-applies any peak broadening to the unbroadened data
* (which is the second column of the ".spc" file).
* The new data is saved as ".blur" file.
* The purpose of this program is to save the user from having
* to re-run DIFFaX simply to modify the broadening component
* of the output.
*
* This program has been compiled and run using Language System's
* Fortran on the Macintosh, and has not been tested on any other
* computer system. Basically, it is nothing more than a useful hack.
* You may wish to tailor it to your needs,
* or force it to comply with your compiler's foibles.
* mmjt 6/89
* mmjt 6/98: Updated. Bug in PV() fixed
************************************************************************
      program Broaden
      implicit none

      integer*4 MAX_SP
      parameter (MAX_SP=20001)
      real*8 spec(MAX_SP), brd_spc(MAX_SP)
      real*8 theta(MAX_SP)
      integer*4 i, j, index, no_records, length, file_cnt
      real*8 theta2_low, theta2_min, theta2_max, d_theta, deg_to_rad
      character*32 filename, savename, rootname, temproot
      character*1 YN, tab
      real*8 pv_u, pv_v, pv_w, pv_gamma, std_dev, h_wdth
      integer*4 in, out, file_io
      integer*4 p_vgt
      logical*4 file_exists, multiple, ok, TRMSPC

      common /spec/ spec, brd_spc
      common /params/ theta2_min, theta2_max, d_theta
      common /PV_blurr/ pv_u, pv_v, pv_w, pv_gamma
      common /G_blurr/ std_dev
      common /L_blurr/ h_wdth
      common /in_out/ in, out
      
	tab = char(9)
	
      in = 6
      out = 6
      file_io = 10
      
      deg_to_rad = datan(1.0D0) / 45.0D0
      
      YN = 'y'
      multiple = .true.
      write(out,*)'Do you wish to blur multiple files? (y/n)'
      read(in,*) YN
      if(YN.ne.'y' .and. YN.ne.'Y') multiple = .false.
      if(multiple) then
        write(out,*)'Enter root name'
        write(out,*)'(ie. if files are ''Fau01.spc'''
        write(out,*)'to ''Fau99.spc'', then enter ''Fau'')'
        read(in,*) rootname
      else
        write(out,*)'Enter filename'
        read(in,*) filename
      endif

    5 write(out,*) 'What type of instrumental broadening?'
      write(out,*) '1: pseudo-Voigt,  2: Gaussian,  3: Lorentzian'
      read(in,*) p_vgt
*
      if(p_vgt .eq. 1) then
        write(out,*) 'Enter pseudo-Voigt u, v, w, gamma'
        pv_u = 0.2D0
        pv_v = -0.18D0
        pv_w  = 0.50000D-01
        pv_gamma  = 0.6D0
        write(out,*) 'Enter u (def = 0.2)'
        read(in,*) pv_u      
        write(out,*) 'Enter v (def = -0.18)'
        read(in,*) pv_v      
        write(out,*) 'Enter w (def = 0.05)'
        read(in,*) pv_w      
        write(out,*) 'Enter gamma (def = 0.6)'
        read(in,*) pv_gamma      
        write(out,'(1x,a,f8.5)') 'u = ',pv_u
        write(out,'(1x,a,f8.5)') 'v = ',pv_v
        write(out,'(1x,a,f8.5)') 'w = ',pv_w
        write(out,'(1x,a,f8.5)') 'gamma = ',pv_gamma      
      else if(p_vgt .eq. 2) then
        write(out,*) 'Enter Gaussian standard deviation'
        read(in,*) std_dev
        write(out,*) 'Gaussian standard deviation = ',std_dev
      else if(p_vgt .eq. 3) then
        write(out,*) 'Enter Lorentzian halfwidth'
        read(in,*) h_wdth
        write(out,*) 'Lorentzian halfwidth = ',h_wdth
      else
        write(out,*) 'Illegal broadening type'
        goto 5
      endif      
      write(out,*) 'okay? (y/n)'
      read(in,*) YN
      if(YN.ne.'y' .and. YN.ne.'Y') then
        write(out,*) 'Re-enter braodening data...'
        goto 5
      endif
*        
*      write(out,*)'multiple = ',multiple

      file_cnt = 0
      if(.not.multiple) then
        inquire(file=filename,exist=file_exists)
         if(.not.file_exists) goto 80
        file_cnt =file_cnt + 1
        open(unit = file_io,file=filename,status = 'old')
        i = 0
10        i = i + 1
        read(file_io,*,end=20) theta(i), spec(i)
        goto 10
20        continue
        close (unit = file_io)
        no_records = i - 1

        write(out,*) ' '
        write(out,*) 'File ''',trim(filename),''' read'
        write(out,*) 'File contains ',no_records,' records'
*
        if(no_records.gt.1) then
          d_theta = (theta(no_records) - theta(1)) / no_records
        elseif(no_records.eq.1) then
          d_theta = theta(no_records)
        else
          write(out,*) 'Insufficient number of records'
          stop
        end if
        d_theta = 0.5D0 * d_theta * deg_to_rad
        theta2_min = theta(1) * deg_to_rad
        theta2_max = theta(no_records) * deg_to_rad
        ok = TRMSPC(theta2_low)
*        write(out,'(1x,a13,f9.5,a7)') 
*     |        'theta2_low = ',theta2_low,' radian'
*        write(out,'(1x,a13,f9.5,a7)') 
*     |        'theta2_min = ',theta2_min,' radian'
*        write(out,'(1x,a13,f9.5,a7)') 
*     |        'theta2_max = ',theta2_max,' radian'
*        write(out,'(1x,a10,f10.7,a7)') 'd_theta = ',d_theta,' radian'
        if(ok) then
          if(p_vgt .eq. 1) call PV(theta2_low)
          if(p_vgt .eq. 2) call GAUSSN(theta2_low)
          if(p_vgt .eq. 3) call LORNZN(theta2_low)
          filename  = trim(filename) //'.blur'
          write(out,*) 'new filename is ',filename
          open(unit = file_io, file = filename, status = 'new')
          do 30 j = 1, no_records
            write(file_io,100) theta(j), tab, brd_spc(j)
30          continue
          close (unit = file_io)
        endif
      endif
      
      if(multiple) then
        do 70 index = 1, 999
          temproot = rootname
          length = len(trim(rootname))
          if(index .lt. 10) then
            write(temproot(length+1:length+2),'(a)') '00'
            write(temproot(length+3:length+3),'(i1)') index
          else if(index .ge. 10 .and. index .lt. 100) then
            write(temproot(length+1:length+1),'(a)') '0'
            write(temproot(length+2:length+3),'(i2)') index
          else
            write(temproot(length+1:length+3),'(i3)') index
          endif
          filename = trim(temproot) //'.spc'
          inquire(file=filename,exist=file_exists)
c          if(.not.file_exists) cycle
          if(.not.file_exists) goto 70
        write(out,*)'filename = ',filename
c        pause
        open(unit = file_io,file=filename,status = 'old')
        file_cnt = file_cnt + 1
        i = 0
40        i = i + 1
        read(file_io,*,end=50) theta(i), spec(i)
        goto 40
50        continue
        close (unit = file_io)
        no_records = i - 1

        if(no_records.gt.1) then
          d_theta = (theta(no_records) - theta(1)) / no_records
        elseif(no_records.eq.1) then
          d_theta = theta(no_records)
        else
          write(out,*) 'Insufficient number of records'
          stop
        end if
        d_theta = 0.5D0 * d_theta * deg_to_rad
        theta2_min = theta(1) * deg_to_rad
        theta2_max = theta(no_records) * deg_to_rad
        ok = TRMSPC(theta2_low)
*        write(out,'(1x,a13,f9.5,a7)') 
*     |        'theta2_low = ',theta2_low,' radian'
*        write(out,'(1x,a13,f9.5,a7)') 
*     |        'theta2_min = ',theta2_min,' radian'
*        write(out,'(1x,a13,f9.5,a7)') 
*     |        'theta2_max = ',theta2_max,' radian'
*        write(out,'(1x,a10,f10.7,a7)') 'd_theta = ',d_theta,' radian'
        if(ok) then
          if(p_vgt .eq. 1) call PV(theta2_low)
          if(p_vgt .eq. 2) call GAUSSN(theta2_low)
          if(p_vgt .eq. 3) call LORNZN(theta2_low)
        
          savename  = trim(temproot) //'.blur'
          write(out,*) 'new filename is ',savename
          open(unit = file_io, file = savename, status = 'new')
          do 60 j = 1, no_records
            write(file_io,100) theta(j), tab, brd_spc(j)
60          continue
*
          close (unit = file_io)
        endif
70      continue
      endif
*
80    if(file_cnt.eq.0) then
        write(out,*) 'Sorry. No files were found.'
      else
        write(out,*) file_cnt, ' file(s) were read.'
      endif
*
  100 format(1x, e12.5, a, g13.6)
      end
*
*********************************
*
      subroutine PV(theta2_low)
      implicit none
*
      real*8 theta2_low
*
      integer*4 MAX_SP
      parameter (MAX_SP=20001)
      real*8 spec(MAX_SP), brd_spc(MAX_SP)
      integer*4 i, j, n_low, n_high
      real*8 PI, th_rng, tn_th, c00, hk_inv
      real*8 k1, k2, k3, k4, k5, pVoigt, const, tmp, speci
      real*8 theta2_min, theta2_max, d_theta, th0
      real*8 pv_u, pv_v, pv_w, pv_gamma
      integer*4 in, out

      common /spec/ spec, brd_spc
      common /params/ theta2_min, theta2_max, d_theta
      common /PV_blurr/ pv_u, pv_v, pv_w, pv_gamma
      common /in_out/ in, out

      write(out,*)'Adding pseudo-Voigt instrumental broadening...'
      PI = 4.0D0 * atan(1.0D0)
      n_low  = int(0.5D0*theta2_low/d_theta) + 1
	
	if(n_low.lt.1) then
        write(out,*) 'PV: n_low = ', n_low
	endif
	
      n_high = int(0.5D0*(theta2_max-theta2_min)/d_theta) + 1
      c00 = 4.0D0 * log(2.0D0)
      const = 360.0D0 * d_theta / PI
      do 10 i = 1, n_high
        brd_spc(i) = 0.0D0
10    continue
      k1 = pv_gamma*2.0D0/PI
      k2 = (1.0D0 - pv_gamma)*sqrt(c00/PI)
      k3 = -c00
      th0 = 0.5D0*theta2_min
      do 20 i = n_low, n_high
C get tan((2theta)/2)
        tn_th = tan(i*d_theta + th0)
        hk_inv = 1.0D0/sqrt((pv_u*tn_th + pv_v)*tn_th + pv_w)
        tmp = (const*hk_inv)**2
        k4 = k1*hk_inv*const
        k5 = k2*hk_inv*const
	  speci = spec(i)
        do 30 j = n_low - i, n_high - i
          th_rng = tmp*j*j
          pVoigt = k4/(1.0D0 + 4.0D0*th_rng) + k5*exp(k3*th_rng)
          brd_spc(i+j) = brd_spc(i+j) + pVoigt*speci
30      continue
20    continue
      return
      end
*
*********************************
*
      subroutine GAUSSN(theta2_low)
      implicit none
*
      real*8 theta2_low
*
      integer*4 MAX_SP
      parameter (MAX_SP=20001)
      real*8 spec(MAX_SP), brd_spc(MAX_SP)
      integer*4 i, j, n_low, n_high, m
      real*8 k1, k2, k3, const, gss, tmp, tmp1, tmp2, PI
      real*8 theta2_min, theta2_max, d_theta, std_dev
      integer*4 in, out
*
      common /spec/ spec, brd_spc
      common /params/ theta2_min, theta2_max, d_theta
      common /G_blurr/ std_dev
      common /in_out/ in, out
*
      PI = 4.0D0 * atan(1.0D0)
*
      if(std_dev.le. 0.0D0) goto 999
      write(out,100)'Adding Gaussian instrumental broadening...'
*
      n_low  = int(0.5D0*theta2_low/d_theta) + 1
      n_high = int(0.5D0*(theta2_max-theta2_min)/d_theta) + 1
      const = 360.0D0 * d_theta / PI
      k1 = const / (sqrt(2.0D0*PI) * std_dev)
      k2 = -0.5D0 * (const / std_dev)**2
*
      do 10 i = 1, n_high
   10   brd_spc(i) = 0.0D0
*
* go out to 40 standard deviations, or to the end of the spectrum
      m = nint(40.0D0 * std_dev / const)
      if(m.gt.n_high) m = n_high
      do 20 i = 0, m-1
        k3 = k2*i*i
        gss = k1*exp(k3)
        do 20 j = n_low+1, n_high
          tmp1 = 0.0D0
          tmp2 = 0.0D0
          if((j-i).gt.n_low)  tmp1 = spec(j-i)
          if((j+i).le.n_high) tmp2 = spec(j+i)
          tmp = tmp1 + tmp2
          if(i.eq.0) tmp = 0.5D0 * tmp
   20     brd_spc(j) = brd_spc(j) + gss * tmp
      return
  999 write(out,100)'Illegal halfwidth in GAUSSN()'
      write(out,100)'Gaussian instrumental broadening not added'
      return
  100 format(1X, A)
      end
*
*********************************
*
      subroutine LORNZN(theta2_low)
      implicit none
*
      real*8 theta2_low
*
      integer*4 MAX_SP
      parameter (MAX_SP=20001)
      real*8 spec(MAX_SP), brd_spc(MAX_SP)
*
      integer*4 i, j, n_low, n_high, m
      real*8 k1, k2, k3, const, lrnz, tmp, tmp1, tmp2, PI
      real*8 theta2_min, theta2_max, d_theta, h_wdth
      integer*4 in, out
*
      common /spec/ spec, brd_spc
      common /params/ theta2_min, theta2_max, d_theta
      common /L_blurr/ h_wdth
      common /in_out/ in, out
*
      PI = 4.0D0 * atan(1.0D0)
*
      if(h_wdth.le. 0.0D0) goto 999
      write(out,100)'Adding Lorentzian instrumental broadening...'
*
      n_low  = int(0.5D0*theta2_low/d_theta) + 1
      n_high = int(0.5D0*(theta2_max-theta2_min)/d_theta) + 1
      const = 360.0D0 * d_theta / PI
      k1 = const / (PI * h_wdth)
      k2 = (const / h_wdth)**2
*
      do 10 i = 1, n_high
   10   brd_spc(i) = 0.0D0
*
* go out to 50 halfwidths, or to the end of the spectrum
      m = nint(50.0D0 * h_wdth / const)
      if(m.gt.n_high) m = n_high
      do 20 i = 0, m-1
        k3 = 1.0D0 + k2*i*i
        lrnz = k1/k3
        do 20 j = n_low+1, n_high
          tmp1 = 0.0D0
          tmp2 = 0.0D0
          if((j-i).gt.n_low)  tmp1 = spec(j-i)
          if((j+i).le.n_high) tmp2 = spec(j+i)
          tmp = tmp1 + tmp2
          if(i.eq.0) tmp = 0.5D0 * tmp
   20     brd_spc(j) = brd_spc(j) + lrnz * tmp
      return
  999 write(out,100)'Illegal halfwidth in LORNZN()'
      write(out,100)'Lorentzian instrumental broadening not added'
      return
  100 format(1X, A)
      end
*
*********************************
*
      logical function TRMSPC(theta2_low)
      implicit none
*
      real*8 theta2_low
*
      integer*4 i, i_min, i_max
      integer*4 MAX_SP
      parameter (MAX_SP=20001)
      real*8 spec(MAX_SP), brd_spc(MAX_SP)
      real*8 theta2_min, theta2_max, d_theta
      integer*4 in, out
*
      common /spec/ spec, brd_spc
      common /params/ theta2_min, theta2_max, d_theta
      common /in_out/ in, out
*
      i_max = int(0.5D0*(theta2_max-theta2_min) / d_theta) + 1
* spec(1) corresponds to the intensity at the origin and is always zero.
      i = 2
   20 i = i + 1
        if(i.ge.i_max+1) then
          write(out,100) 'No peaks were found in spectrum.'
          goto 900
        endif
* locate the first minimum after the huge peak at the origin
        if(spec(i).le.spec(i-1))
     |goto 20
      i_min = i - 1
*
* NOTE: The absolute angle is th2_low + th2_min
      theta2_low = i_min * d_theta
*
  900 TRMSPC = .true.
*
      return
  100 format(1x, a)
      end
*