subroutine fh4(x,yt,xt,y,dydx,indh)
*
* Lagrange interpolating polynomial
*
      implicit double precision(a-h,o-z)

      real*8 yt(4),xt(4),xmt(4)

      a1=yt(1)/((xt(1)-xt(2))*(xt(1)-xt(3))*(xt(1)-xt(4)))
      a2=yt(2)/((xt(2)-xt(1))*(xt(2)-xt(3))*(xt(2)-xt(4)))
      a3=yt(3)/((xt(3)-xt(1))*(xt(3)-xt(2))*(xt(3)-xt(4)))
      a4=yt(4)/((xt(4)-xt(1))*(xt(4)-xt(2))*(xt(4)-xt(3)))

      do 100 k=1,4
100   xmt(k)=x-xt(k)
      y=a1*xmt(2)*xmt(3)*xmt(4)+a2*xmt(1)*xmt(3)*xmt(4) +
     1 a3*xmt(1)*xmt(2)*xmt(4)+a4*xmt(1)*xmt(2)*xmt(3)
      if( indh .eq. 0) return 
      fr=xmt(2)/(xt(3)-xt(2)) 
      if(xmt(4)*xmt(3) .lt. 0.) go to 150
      dydx1=(xt(1)-xt(4))*a1*(xmt(2)+xmt(3))
     1  +(xt(2)-xt(4))*a2*(xmt(1)+xmt(3))
     2 +(xt(3)-xt(4))*a3*(xmt(2)+xmt(1))
      go to 200
150   dydx1=(yt(4)-yt(3))/(xt(4)-xt(3)) 
      fr=-xmt(4)/(xt(4)-xt(3))
200   if(xmt(1)*xmt(2) .lt. 0.) go to 250
      dydx2 = (xt(2)-xt(1))*a2*(xmt(3)+xmt(4)) +
     1  (xt(3)-xt(1))*a3*(xmt(2)+xmt(4)) +
     2  (xt(4)-xt(1))*a4*(xmt(3)+xmt(2))
      go to 300
250   dydx2=(yt(2)-yt(1))/(xt(2)-xt(1)) 
      fr=-xmt(2)/(xt(2)-xt(1))
300   dydx=dydx2*fr+dydx1*(1.-fr)

      return 
      end

************************************************************************