program main

  use accuracy
  use xy3_alpha
  implicit none

  integer(ik), parameter :: natoms=4
  integer(ik) :: iounit=10, info, iline, iterm, i, j, k, l, m, iopt, iatom
  real(ark) :: xyz(natoms,3), xyz_angs(natoms,3), f(maxnelem_sym), enr, dipole(3), polariz(3,3), hyperpolariz(3,3,3), &!
               second_hyperpolariz(3,3,3,3), polariz_lin(maxnelem_sym), hyperpolariz_lin(maxnelem_sym), &!
               second_hyperpolariz_lin(maxnelem_sym), angs=0.529177209_ark
  character(cl) :: param_fname, data_fname, lab
  type(ext_func) :: extf

  ! read name of file with tensor function parameters
  ! and file with data points (molecular Cartesian coordinates)

  call getarg(1, param_fname)
  call getarg(2, data_fname)


  ! read file with electric tensor function parameters

  open(iounit,file=param_fname,form='formatted',status='old',iostat=info)
  if (info/=0) then
    write(out, '(/a,1x,a,1x,a)') 'error: file', trim(param_fname), 'is not found'
    stop
  endif

  iline = 0

  read(iounit,*,iostat=info) extf%rank
  iline = iline + 1
  if (info/=0) then
    write(out, '(/a,1x,a,1x,a,1x,i6)') 'error while reading file', trim(param_fname), ', line =', iline
    stop
  endif

  extf%nelem = 3**extf%rank
  select case(extf%rank)
  case default
    write(out, '(/a,1x,i3,1x,a)') 'error: illegal tensor rank =', extf%rank, ', expected rank=2'
    stop
  case(2) ! polarizability
    extf%nelem_sym = 6
  end select

  read(iounit,*,iostat=info) extf%nterms(1:extf%nelem_sym)
  iline = iline + 1
  if (info/=0) then
    write(out, '(/a,1x,a,1x,a,1x,i6)') 'error while reading file', trim(param_fname), ', line =', iline
    stop
  endif

  do i=1, extf%nelem_sym
    do iterm=1, extf%nterms(i)
      read(iounit,*,iostat=info) lab, iopt, extf%coef(iterm,i)
      iline = iline + 1
      if (info/=0) then
        write(out, '(/a,1x,a,1x,a,1x,i6)') 'error while reading parameters from file', trim(param_fname), ', line =', iline
        stop
      endif
    enddo
  enddo

  close(iounit)


  ! read file with data points and compute electric-field tensor
  ! compare results with reference ab initio values

  open(iounit,file=data_fname,form='formatted',status='old',iostat=info)
  if (info/=0) then
    write(out, '(/a,1x,a,1x,a)') 'error: file', trim(data_fname), 'is not found'
    stop
  endif

  iline = 0

  select case(extf%rank)
  case default
    write(out, '(/a,1x,i3,1x,a)') 'error: illegal rank =', extf%rank, ', expected rank=1..4'
    stop

  case(2) ! polarizability

    do while(.true.)
      read(iounit,*,iostat=info) (xyz(iatom,1:3), iatom=1, natoms), (polariz(i,1:3), i=1,3), enr
      iline = iline + 1
      if (info>0) then
        write(out, '(/a,1x,a,1x,a,1x,i6)') 'error while reading molecular cooridnates from file', trim(data_fname), ', line =', iline
        stop
      elseif (info<0) then
        exit
      endif
      xyz_angs = xyz !!* angs
      call xy3_polariz_symadap(extf%nelem_sym, extf, natoms, xyz_angs, f(1:extf%nelem_sym))
      k = 0
      do i=1, 3
        do j=i, 3
          k = k + 1
          polariz_lin(k) = polariz(i,j)
        enddo
      enddo
      write(out, '(1x,i6,3x,4(3x,3(1x,f12.8)),3x,6(1x,es16.8),3x,6(1x,f10.6),3x,f12.4)') &!
      iline, (xyz(iatom,1:3), iatom=1, natoms), f(1:extf%nelem_sym), (f(i)-polariz_lin(i), i=1, extf%nelem_sym), enr
    enddo

  end select

  close(iounit)

end program main
