function phi(x)
common B, C, ieq

x2 =x*x

if(ieq ==1) then

phi = log(4*x2+3)

elseif(ieq ==2) then

phi = -1.2*x2 +1.9*log(4*x2+3)

else
x4 =x2*x2

phi = ( 45*x4 -109.5*x2 +117.125*log(4*x2+3) )/35
endif

end


! ---------------
real function k(x)
common B, C, ieq

x2 =x*x
x3 =x2*x

if(ieq ==1) then
k = 4*x/(4*x2+3)

elseif(ieq ==2) then

k = (-4.8*x3 +4*x)/(4*x2+3)

else	! ieq =3
x5 =x3*x2;

k =(360*x5 -168*x3 +140*x)/(4*x2+3)/35
endif

end
