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 ************************************************************************