MODULE interpol
	IMPLICIT NONE
	CONTAINS
	FUNCTION interp(n,x,y,x1) result(y1)
		INTEGER,INTENT(in)::n
		REAL(8),INTENT(in)::x(n),y(n),x1
		REAL(8)::y1
		INTEGER::i,i1,i2
		IF(x1<=x(1))THEN
			y1=y(1)
			RETURN
		ELSE IF(x1>=x(n))THEN
			y1=y(n)
			RETURN
		ELSE IF(x1<=x(2))THEN
			y1=lin_interp(x(1),y(1),x(2),y(2),x1)
			RETURN
		ELSE IF(x1>=x(n-1))THEN
			y1=lin_interp(x(n-1),y(n-1),x(n),y(n),x1)
		END IF
		i1=1
		i2=n
		DO while(i2-i1 > 1)                  					      					 	!Lion Hunt
			i=(i1+i2)/2
			IF(x1 > x(i))THEN
				i1=i
			ELSE
				i2=i
			END IF
		END DO
		y1=lag3(x(i1-1:i2+1),y(i1-1:i2+1),x1)
		RETURN
	END FUNCTION interp

	FUNCTION lag3(x,y,x1) result(y1)                                                       	!Third order Lagrange interpolation
		REAL(8),INTENT(in)::x(4),y(4),x1
		REAL(8):: y1,f1,f2,f3,f4
		f1=((x1-x(2))*(x1-x(3))*(x1-x(4)))/((x(1)-x(2))*(x(1)-x(3))*(x(1)-x(4)))
		f2=((x1-x(1))*(x1-x(3))*(x1-x(4)))/((x(2)-x(1))*(x(2)-x(3))*(x(2)-x(4)))
		f3=((x1-x(1))*(x1-x(2))*(x1-x(4)))/((x(3)-x(1))*(x(3)-x(2))*(x(3)-x(4)))
		f4=((x1-x(1))*(x1-x(2))*(x1-x(3)))/((x(4)-x(1))*(x(4)-x(2))*(x(4)-x(3)))
		y1=y(1)*f1+y(2)*f2+y(3)*f3+y(4)*f4
		RETURN
	END FUNCTION lag3

	FUNCTION hermit3(x,y,x1) result(y1)                                                    	!Third order Hermit interpolation
		REAL(8),INTENT(in)::x(4),y(4),x1
		REAL(8)::y1,dy1,dy2,g1,g2,f1,f2
		dy1=(y(3)-y(1))/(2*(x(3)-x(1)))
		dy2=(y(4)-y(2))/(2*(x(4)-x(2)))
		g1=((x1-x(2))*((x1-x(3))*(x1-x(3))))/((x(2)-x(3))*(x(2)-x(3)))
		g2=(((x1-x(2))*(x1-x(2)))*(x1-x(3)))/((x(2)-x(3))*(x(2)-x(3)))
		f1=(((x1-x(3))*(x1-x(3)))/((x(2)-x(3))*(x(2)-x(3)))-(2*g1)/(x(2)-x(3)))
		f2=(((x1-x(2))*(x1-x(2)))/((x(3)-x(2))*(x(3)-x(2)))-(2*g2)/(x(3)-x(2)))
		y1=dy1*g1+dy2*g2+y(2)*f1+y(3)*f2
		RETURN
	END FUNCTION hermit3

	FUNCTION lin_interp(x1,y1,x2,y2,x) result(y)                                            !linear interpolation
		REAL(8),INTENT(in):: x1,y1,x2,y2,x
		REAL(8):: y
		y=y1+(y2-y1)/(x2-x1)*(x-x1)
		RETURN
	END FUNCTION lin_interp
END MODULE interpol
