      program apply_hvca
      implicit none
c---- constants
      integer   maxOverflow
      parameter (maxOverflow=10000)
c---- functions
      integer   length
      integer   hkl
c---- local vars
      character filename*80
      integer   lf
      character string*80
      integer   unit
      logical   todo(MaxOverflow)
      integer   overflows(MaxOverflow)
      integer   t
      integer   iostat
      integer   nr,h,k,l,code
      real      a,p
      real      ri,rs
      integer   nin
      integer   nout
      character infile*80
      integer   li
      character outfile*80
      integer   lo
      integer   error
      integer   out,in
      integer   nd
      integer   ier
      logical   found
      integer   i,j,m
      integer   nf
      real      minlow,maxlow
      character line*80
      integer   ntoo_strong
      real      ymsobs,zmsobs,rotobs
      character formatstring*80
      real      intensity,sigma
      character message*256
 
c
      call finish(-6)
      call str_real_len(4)
      call get_commandline
      formatstring = '!! FORMAT h k l Int SigInt Code Yms Zms Rot'
      error = 1
 
      call get_inp_default(outfile,nd,'Output file name [FINAL.FAST]',
     1 *900)
      call build_file(outfile,'FINAL','.FAST',outfile,lo)
      call create_file(outfile(1:lo),out,message,iostat)
      call str_write(6,message)
      if(iostat.ne.0) goto 900
 
      call build_file_type(outfile,'HVCA','.FAC',filename,lf)
      call open_file(filename(1:lf),unit,message,iostat)
      call str_write(6,message)
      if(iostat.ne.0) then
        call str_write(6,'Assuming no attenuation data available')
        a = -1.0
        p = -1.0
      else
        read (unit,*,iostat=iostat) a,p
        if(iostat.ne.0) then
          call str_write(6,'Empty file: '//filename)
          goto 900
        endif
        close (unit)
      endif
      call str_2real_write(6,a,p,'I(High) = # * I(Low)^#')
 
      call get_inp_default(infile,nd,'High FAST file [HIGH]',*900)
      call build_file(infile,'HIGH','.FAST',infile,Li)
      call open_file(infile(1:li),in,message,iostat)
      call str_write(6,message)
      if(iostat.ne.0) goto 900
      call str_write(6,'Opened '//infile(1:Li))
      read (in,'(a)',iostat=iostat) line
      if(line(1:8).ne.'!!_BEGIN') then
        call str_write(6,'Expected !!_BEGIN')
        call str_write(6,'Found '//line)
        ier = 1
        goto 900
      endif
      call begin(out)
      if(iostat.eq.0) then
        read (in,'(a)',iostat=iostat) line
        do while(line(1:6).ne.'!!_END')
          if(line(1:9).ne.'!! FORMAT') call str_write(out,line)
          read (in,'(a)',iostat=iostat) line
        enddo
        call str_write(out,formatstring)
        call str_write(out,line)
      endif
 
      if(iostat.eq.0) read(in,*,iostat=iostat) nr,h,k,l,intensity,sigma,
     1 code,ymsobs,zmsobs,rotobs
      if(iostat.ne.0) then
        call str_write(6,'Empty file:'//infile(1:Li))
        goto 900
      endif
 
      t = 0
      nout = 0
      nin = 0
      ntoo_strong = 0
      call str_write(6,'Reading '//infile(1:li)//'....')
      do while(iostat.eq.0)
        nin = nin + 1
        if(code.ge.32768) then
          if(a.gt.0.0) then
            t = t + 1
            if(t.gt.MaxOverflow) then
              call str_write(6,
     1         'Too many overflows. Increase MaxOverflow')
              goto 900
            endif
            todo(t) = .true.
            overflows(t) = hkl(h,k,l)
          else
            ntoo_strong = ntoo_strong + 1
            nout = nout + 1
            write (out,300) h,k,l,intensity,sigma,code,
     1       ymsobs,zmsobs,rotobs
          endif
        else
          nout = nout + 1
          write (out,300) h,k,l,intensity,sigma,code,
     1     ymsobs,zmsobs,rotobs
300       format(3i5,2(1x,f8.2),1x,i6,2(1x,f6.1),1x,f9.3)
        endif
        read(in,*,iostat=iostat) nr,h,k,l,intensity,sigma,code,
     1   ymsobs,zmsobs,rotobs
      enddo
      close(in)
      call str_num(nin,'# reflections from '//infile(1:li),string)
      call str_write(6,string)
      call log_add(string)
      call str_num(nout,'# reflections to '//outfile(1:lo),string)
      call str_write(6,string)
      call log_add(string)
      if(a.gt.0.0) then
        call str_num(t,'# overflow reflection^',string)
        call str_write(6,string)
        call log_add(string)
      else
        call str_num(ntoo_strong,'# too strong reflection^',string)
        call str_write(6,string)
        call log_add(string)
      endif
 
      if(t.gt.0.and.a.gt.0.0) then
        minlow = 999999.0
        maxlow =      0.0
        call get_inp_default(infile,nd,'Low FAST file [LOW]',*900)
        call build_file(infile,'LOW','.FAST',infile,li)
        call open_file(infile(1:li),in,message,iostat)
        call str_write(6,message)
        if(iostat.ne.0) goto 900
 
        if(iostat.eq.0) then
          read (in,'(a)',iostat=iostat) line
          do while(line(1:6).ne.'!!_END')
            if(line(1:9).ne.'!! FORMAT') call str_write(out,line)
            read (in,'(a)',iostat=iostat) line
          enddo
          call str_2real_write(out,a,p,'!! I(HIGH) = # * I(Low)^#')
          call str_write(out,formatstring)
          call str_write(out,line)
        endif
 
        if(iostat.eq.0) read(in,*,iostat=iostat) nr,h,k,l,intensity,
     1   sigma,code,ymsobs,zmsobs,rotobs
        if(iostat.ne.0) then
          call str_write(6,'Empty file:'//infile(1:li))
          goto 900
        endif
        nin = 0
        nf = 0
        call str_write(6,'Reading '//infile(1:li)//'....')
        do while(iostat.eq.0)
          nin = nin + 1
          m = hkl(h,k,l)
          found = .false.
          j = 0
50        continue
            j = j + 1
            if (m.eq.overflows(j).and.todo(j)) found = .true.
          if(.not.found.and.j.lt.t) goto 50
          if(found) then
            todo(j) = .false.
            minlow = min(minlow,intensity)
            maxlow = max(maxlow,intensity)
            nf = nf + 1
            if(intensity.lt.0.0) then
              intensity=-intensity
              ri = -a*(intensity**p)
              intensity=-intensity
            else
              ri = a*(intensity**p)
            endif
            rs = a*(sigma**p)
            write (out,300) h,k,l,ri,rs,code,
     1       ymsobs,zmsobs,rotobs
            if(code.ge.32768) ntoo_strong = ntoo_strong + 1
            nout = nout + 1
          endif ! found
          read(in,*,iostat=iostat) nr,h,k,l,intensity,sigma,code,
     1     ymsobs,zmsobs,rotobs
        enddo
        close(in)
        call str_num(nin,'# reflection^ from '//infile(1:li),string)
        call str_write(6,string)
        call log_add(string)
        call str_num(nf,'# attenuated reflection^ used',string)
        call str_write(6,string)
        call log_add(string)
        call str_num(ntoo_strong,'# too strong reflection^',string)
        call str_write(6,string)
        call log_add(string)
        call str_num(nout,'Total:# reflections written to '//
     1   outfile(1:lo),string)
        call str_write(6,string)
        call log_add(string)
        call str_2real_write(6,minlow,maxlow,
     1   'Intensity range for low values of attenuated refls=')
        j = 0
        do i=1,t
          if(todo(i)) then
            j = j + 1
            if(j.lt.20) then
              call extract(overflows(i),h,k,l)
              call str_3integer_write(6,h,k,l,'Low hkl missing:')
            else if(j.eq.20) then
              call str_write(6,'Missing List truncated...')
            endif
          endif
        enddo
        call str_num(j,'# reflection^ missing',string)
        call str_write(6,string)
        call log_add(string)
      endif ! t
c      write(out,120)
c120   format(1x)
      close(out)
      error = 0
900   continue
      call log_close
      call finish(6)
      if(error.eq.0) then
        call exit(1)
      else
        call exit(2)
      endif
      end
 
      integer function hkl(h,k,l)
      implicit none
      integer   h,k,l
c---- constants
      integer   maxindex
      parameter (maxindex = 100)
      integer   factor
      parameter (factor=2*maxindex)
c
      hkl = (h+maxindex)*factor**2 +
     1      (k+maxindex)*factor +
     2      (l+maxindex)
      return
      end
 
      subroutine extract(n,h,k,l)
      integer n
      integer h,k,l
c---- constants
      integer   maxindex
      parameter (maxindex = 100)
      integer   factor
      parameter (factor=2*maxindex)
c---- local vars
      integer j
c
      j = n
      l = mod(j,factor) - maxindex
      j = j/factor
      k = mod(j,factor) - maxindex
      j = j/factor
      h = j - maxindex
      return
      end
 
