program Invinter
  ! Invinter -- Inverse interpolate in 1-D
  !
  ! <in.H Invinter o1= d1=  >out.H
  !  input- 
  ! in.H in(2,n2) = ordinate, coordinate pairs.
  ! n1 output model size.
  ! niter=n1 iteration count
  !  output-
  ! out.H inverse linear interpolated model mesh mm(m1)

  use sep
  use invint_mod

  implicit none
  integer                        :: n2, nd, m1, na, niter, i2, randseed, nframe
  real, dimension (:),   pointer :: mm, aa, coord, ord
  real, dimension (:,:), pointer :: pair
  real                           :: o1, d1, eps, randscale

  call sep_init ()

  ! create model
  call from_par ("n1", m1)
  allocate ( mm(m1) )
  call from_param("niter", niter, m1) ; niter = 1 + niter*3/2 
  call from_history ("o1", o1)
  call from_history ("d1", d1)
  call from_param( "nframe", nframe, 1)
  call from_param( "scale", randscale, 1.0)

  call to_history ("n2", nframe)

  call from_par ("eps",eps,0.)
  if (eps == 0.) then                       ! Epsilon is in the filter history file!
     ! read filter
     call from_aux ("filtin","epsilon",eps,1.) 
     call from_aux ("filtin","n1",na)
     allocate (aa (na)) 
     call sep_read (aa,"filtin")
  end if

  ! read data
  call from_history (nd, n2)
  if (n2 /= 2) call erexit ('bad n2')
  allocate( pair(nd,2), coord (nd), ord (nd))
  call sep_read (coord)
  call sep_read (ord)
  pair(:,1) = coord;  pair(:,2) = ord;
  call sep_close ()

  do i2=1,nframe
     randseed = 177 + 100*i2
     call init_random( .true., randseed, randscale) 
     call invint( niter, pair, nd, o1, d1, aa, mm, eps)
     call sep_write (mm)
  end do

  deallocate (mm, ord, coord, aa, pair)

  call exit (0) 
end program Invinter
