
	subroutine givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	implicit NONE
	double precision xa1,xa2,xa3,ya1,ya2,ya3,a,b,c
	double precision den
	
	den=(xa1-xa2)*(xa2-xa3)*(xa3-xa1)
	a = -ya1*(xa2-xa3)/den
	a = a-ya2*(xa3-xa1)/den
	a = a-ya3*(xa1-xa2)/den
      
	b =  ya1*(xa2*xa2-xa3*xa3)/den
	b = b+ya2*(xa3*xa3-xa1*xa1)/den
	b = b+ya3*(xa1*xa1-xa2*xa2)/den
      
	c = -ya1*(xa2-xa3)*xa2*xa3/den
	c = c-ya2*(xa3-xa1)*xa1*xa3/den
	c = c-ya3*(xa1-xa2)*xa1*xa2/den
	return
	end 




	double precision function fitted(x)
	implicit NONE
	integer Np,Np2
	parameter(Np=258,Np2=516)
	double precision x
	double precision mass(Np), deltaew(Np), yy(Np2)
	double precision a,b,c,xa1,xa2,xa3,ya1,ya2,ya3
	integer i, n
	double precision ris
	double precision diff, diffL, diffR
	integer position, i1, i2, i3
	external givecoeff
#include "./electroweak.h"

	do i=1,Np
	   mass(i) = yy(2*i -1)
	   deltaew(i) = yy(2*i)
	enddo

	if ((x.lt.100d0).or.(x.gt.500d0)) then 
	   print*, "EWK correction not available: set to zero"
	   fitted = 0d0
	   return
	endif

	position = Np
	do i=1,Np
	   diff = x-mass(i)
	   if (diff.lt.0d0) then 
	      position = i-1
	      goto 10
	   endif
	enddo
 10	continue

	if (position.eq.1) then 
	   i1 = position
	   i2 = position+1
	   i3 = position+2
	elseif (position.eq.Np) then 
	   i1 = position
	   i2 = position+1
	   i3 = position+2
	else
	   diffL = x-mass(position)
	   diffR = mass(position+1)-x
	   if (diffL.lt.diffR) then 
	      i1 = position-1
	      i2 = position
	      i3 = position+1
	   else
	      i1 = position
	      i2 = position+1
	      i3 = position+2
	   endif
	endif

	xa1=mass(i1)
	xa2=mass(i2)
	xa3=mass(i3)
	ya1=deltaew(i1)
	ya2=deltaew(i2)
	ya3=deltaew(i3)
	
C	print*, "mass =", x
C	print*, xa1, xa2, xa3
C	print*, ya1, ya2, ya3

	call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)

	ris = c + b*x +a*x*x

	fitted = ris
	return
	end
