program answer

*** compute a robust position

      implicit double precision(a-h,o-z)
      parameter(mx=90000)
      dimension dlas(mx),dlos(mx),dhts(mx)
      character*80 fname

      write(*,*) 'answer -- average values -- 25mar02'

      factor=2.5d0

      write(*,'(a$)') ' Enter input file name = '
      read(*,'(a80)') fname

      lin=1
      open(lin,file=fname,form='formatted',status='old')

      n=1
    1 read(lin,*,err=1,end=9) 
     *  t,dlas(n),dlos(n),dhts(n)
      n=n+1
      if(n.gt.mx) go to 9
      go to 1
    9 nall=n-1

*** first rejection pass

      aas=0.d0
      aos=0.d0
      ahs=0.d0

      tol1=9.9d20
      tol2=9.9d20
      tol3=9.9d20

*** iteration

      do 100 iloop=1,10

      das=0.d0
      dos=0.d0
      dhs=0.d0
      da2=0.d0
      do2=0.d0
      dh2=0.d0
      n=0
      do 10 i=1,nall
        if(dabs(dlas(i)-aas).le.tol1.and.
     *     dabs(dlos(i)-aos).le.tol2.and.
     *     dabs(dhts(i)-ahs).le.tol3) then
          das=das+dlas(i)
          dos=dos+dlos(i)
          dhs=dhs+dhts(i)
          da2=da2+dlas(i)*dlas(i)
          do2=do2+dlos(i)*dlos(i)
          dh2=dh2+dhts(i)*dhts(i)
          n=n+1
        endif
   10 continue
      en=dble(n)
      aas=das/en
      aos=dos/en
      ahs=dhs/en
      ras=dsqrt((da2-das*das/en)/en)
      ros=dsqrt((do2-dos*dos/en)/en)
      rhs=dsqrt((dh2-dhs*dhs/en)/en)
      write(*,'(2f15.9,f11.4,i7)') aas,aos,ahs,n
      write(*,'(2f15.9,f11.4   )') ras,ros,rhs

      tol1=factor*ras
      tol2=factor*ros
      tol3=factor*rhs

  100 continue

      end