      program hvca_factor
      implicit none
c---- constants
      logical   debug
      parameter (debug=.false.)
c---- functions
      integer   length
c---- local vars
      character string*80
      character d*80
      integer   nd
      character infile1*80
      integer   li1
      character infile2*80
      integer   li2
      character outfile*80
      integer   lo
      character filename*80
      integer   lf
      integer   iostat
      integer   m1,H1,K1,L1,code1,nr1
      real      I1,S1
      integer   m2,H2,K2,L2,code2,nr2
      real      I2,S2
      real      i3
      integer   error
      integer   in1,in2,out,out2
      integer   nin1,nin2,nout,npair
      real      det,p,a
      integer   nx
      real      x,y
      real      sx,sy,sxx,sxy
      integer   unit
      character message*256
c
      call finish(-6)
      call str_unit(6)
      call str_real_len(2)
      error = 1
      call get_commandline
      call get_inp_default(d,nd,'High HKL filename [HIGH.SCR]',*900)
      call build_file(d,'HIGH','.SCR',infile1,li1)
      call open_file(infile1(1:li1),in1,message,iostat)
      call str_write(6,message)
      if(iostat.ne.0) goto 900
      read(in1,*,iostat=iostat) M1,H1,K1,L1,I1,S1,code1,nr1
      if(iostat.ne.0) then
        call str_write(6,'Empty file:'//infile1(1:li1))
        goto 900
      endif
      nin1 = 1
 
      call get_inp_default(d,nd,'Low HKL filename [LOW.SCR]',*900)
      call build_file(d,'LOW','.SCR',infile2,li2)
      call open_file(infile2(1:li2),in2,message,iostat)
      call str_write(6,message)
      if(iostat.ne.0) goto 900
      read (in2,*,iostat=iostat) M2,H2,K2,L2,I2,S2,code2,nr2
      if(iostat.ne.0) then
        call str_write(6,'Empty file '//infile2(1:li2))
        goto 900
      endif
      nin2 = 1
 
      call get_inp_default(outfile,nd,
     1 'Combined HKL filename [TEST.HKL2]',*900)
      call build_file(outfile,'HVCA','.HKL2',outfile,lo)
      call create_file(outfile(1:lo),out,message,iostat)
      call str_write(6,message)
      if(iostat.ne.0) goto 900
      nout = 0
 
      nx = 0
      sx = 0.0
      sy = 0.0
      sxx = 0.0
      sxy = 0.0
      iostat = 0
      npair  = 0
      do while(iostat.eq.0)
        if(debug) write (6,600) m1,h1,k1,l1,m2,h2,k2,l2,code2,nr2
600     format(1x,i10,3i5,5x,i10,3i5,i6,i6)
        if(m1.eq.m2) then
          npair = npair + 1
          if(i1.gt.1.0.and.i2.gt.1.0) then
c-     y = A * x^p   y = HVCAHIGH  x = HVCALOW
            nx = nx + 1
            x = log(i2)
            y = log(i1)
            sx = sx + x
            sxx = sxx + x*x
            sy = sy + y
            sxy = sxy + x*y
            write (out,300) h1,k2,l1,i1,s1,i2,s2,
     1       1000.0*s1/i1,1000.0*s2/i2,code1,code2,nr1,nr2
300         format(3i4,4f10.2,2f8.2,2i6,2i6)
            nout = nout + 1
          endif
          read(in1,*,iostat=iostat) M1,H1,K1,L1,I1,S1,code1,nr1
          if(iostat.eq.0) then
            nin1 = nin1 + 1
            read(in2,*,iostat=iostat) M2,H2,K2,L2,I2,S2,code2,nr2
            if(iostat.eq.0) nin2 = nin2 + 1
          endif
        else if(m1.lt.m2) then
          read(in1,*,iostat=iostat) M1,H1,K1,L1,I1,S1,code1,nr1
          if(iostat.eq.0) nin1 = nin1 + 1
        else
          read (in2,*,iostat=iostat) M2,H2,K2,L2,I2,S2,code2,nr2
          if(iostat.eq.0) nin2 = nin2 + 1
        endif
      enddo ! while iostat
      close(in1)
      close(in2)
c      close(out)
      call str_num(nin1,'# reflection^ from '//infile1(1:li1),string)
      call str_write(6,string)
      call log_add(string)
      call str_num(nin2,'# reflection^ from '//infile2(1:li2),string)
      call str_write(6,string)
      call log_add(string)
      call str_num(npair,'# pair^',string)
      call str_write(6,string)
      call log_add(string)
      call str_num(nx,'# pair^ used',string)
      call str_write(6,string)
      call log_add(string)
      call str_num(nout,'# reflection^ to '//outfile(1:lo),string)
      call str_write(6,string)
      call log_add(string)
 
      det = (sxx*nx-sx*sx)
      p = (sxy*nx-sx*sy)/det
      A = (sxx*sy-sxy*sx)/det
      a = exp(a)
      call str_real_len(4)
      call str_2real(a,p,'I(high) = # * I(low)^##',string)
      call str_out(string)
      call log_add(string)
      call str_flush
 
      call build_file_type(outfile,'HVCA','.FAC',filename,lf)
      call create_file(filename(1:lf),unit,message,iostat)
      call str_write(6,message)
 
      if(iostat.ne.0) goto 900
      call str_2real_write(unit,a,p,'# #')
      call str_2real(a,p,'I(high) = # * I(low)^##',string)
      close(unit)
c
c      rewind out
c      call build_file_type(outfile,'HVCA','.HKL4',outfile,lo)
c      call create_file(outfile(1:lo),out2,message,iostat)
c      call str_write(6,message)
c      if(iostat.ne.0) goto 900
c      read (out,300,iostat=iostat) h1,k1,l1,i1,s1,i2,s2
c      do while(iostat.eq.0)
c        i3 = a*(i2**p)
c        write(out2,310) h1,k1,l1,i1,i2,i3
c310     format(3i4,4f10.2)
c        read (out,300,iostat=iostat) h1,k1,l1,i1,s1,i2,s2
c      enddo
c      close(out2)
 
      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
