subroutine rterr(x,y,sigmay,npts,p,dp,yfit)
c
c     This basic subroutine calculates the coordinates of the 
c     X^2=chilim point in a given direction from the X^2 minimum.
c
      real*8 x(*),y(*),sigmay(*),yfit(*)
      real*8 p(*),dp(*),dpstep,a(2)
      real*8 chisq,chisqo,chifr,chilim,chimin
      integer*4 maxitr
      parameter (maxitr=100)

      modef=0 ! Gaussian model in FCHISQ
      nfree=npts-2
      if(nfree.le.0)goto 100
      chifr=0.01D0
      dpstep=dp(1)/10.D0
      chilim=1.D0+(1.0D0/dble(nfree))
      chimin=1.D0
c
c     Evaluate start chi square at estimated X^2=chilim point
   11 call rtshift(p,dp,a)
      call rtfuncs(npts,x,yfit,a)
      call fchisq(y,sigmay,npts,nfree,yfit,chisqo,modef)
      if(chisqo.le.chimin)then
       chimin=1.1D0
       dp(1)=dp(1)+dpstep
       goto 11
      endif
      if(chisqo.gt.chilim)dpstep=-dpstep
      dp(1)=sqrt((chilim-1.D0)/((chisqo-1.D0)/(dp(1)**2)))-dpstep
c
      DO ITER=1,MAXITR
       if(iter.eq.20)dpstep=2.D0*dpstep
       if(iter.eq.30)dpstep=2.D0*dpstep
       if(iter.eq.40)dpstep=2.D0*dpstep
       if(iter.eq.50)dpstep=2.D0*dpstep
       if(iter.eq.60)dpstep=2.D0*dpstep
       if(iter.eq.70)dpstep=2.D0*dpstep
       dp(1)=dp(1)+dpstep
       call rtshift(p,dp,a)
       call rtfuncs(npts,x,yfit,a)
       call fchisq(y,sigmay,npts,nfree,yfit,chisq,modef)
       if(abs(1.D0-chisq/chilim).lt.chifr)goto 111
       if((chisqo-chilim)*(chisq-chilim).lt.0)then
        dpstep=-dpstep/2.D0
       else if(chisqo.ge.chilim.and.chisq.gt.chisqo)then
        dpstep=-dpstep
       else if(chisqo.le.chilim.and.chisq.lt.chisqo)then
        dpstep=-dpstep
       endif
        chisqo=chisq
      ENDDO
      write(6,*)'RTERR: had problems with one or more angles!'
      dp(1)=4.848D-9   ! set to some value, so that crash is avoided
      goto 111
c
  100 write(6,*)' Too few data points! Why?'
c
  111 return
      end