import numpy as np
import matplotlib.pyplot as plt

C0_const =  (np.log(np.pi*2)-np.euler_gamma)/2.
low_array = np.array( [ 0., np.pi**2/72.,  -np.pi**4/10800., np.pi**6/238140., -np.pi**8/2268000., np.pi**10/12349260. , -np.pi**12*477481/20917681785000., np.pi**14/109459350.] )
lowfit = np.poly1d( np.flip(low_array) )
high_array = np.array( [ 7.88123841835372274931e-05, -6.421972828775741795e-04, 2.311480349378722439e-03, -4.682647339174743698e-03, 5.328447750470300531e-03, -1.867534587731306452e-03, -4.086589757373656692e-03, 6.619996012020041599e-03, -1.014525569624682121e-03, -8.748354666170923832e-03, -4.958538288302578824e-05, 1.370838127635092829e-01, -4.277093729386082283e-07, 1.378955328194695811e-08 ] )

highfit = np.poly1d( high_array) 

def Sofa(a):
    if(a < 0.):
        return Sofa(-a)
    if( a > 1.):
        return  a*Sofa(1./a) + (a+1.)/2.*np.log(a) + C0_const * (1-a)
    if a < 0.2 :
        asq = a ** 2
        return lowfit(asq)
    else:
        return highfit(a)

    
#a = float(input("Enter a = \n"))
#print("S(a) = ", Sofa(a))
alist = np.linspace(0., 2.)
slist = np.array([ Sofa(a) for a in alist ])
plt.plot(alist, slist)
plt.show()