#!/usr/bin/sage

# Schmidt Quasi-Normalized Legendre Functions for Spherical Harmonics

x = var('x');

def polypart(n,m):
    pol = (x^2-1)^n
    for I in range(1,n+1):
        pol = diff(pol,x)
    for I in range(1,m+1):
        pol = diff(pol,x)
    return pol

def normpart(n,m):
    if(m==0):
        nrm = 1/2^n/factorial(n)
    else:
        nrm = (1/2^n/factorial(n)) * sqrt(2*factorial(n-m)/factorial(n+m))
    return nrm

NN = 20


lns = []
coefs = []
nmaxcoefs = 0
for I in range(0,NN+1):
    for J in range(0,I+1):
        pol = polypart(I,J)
        nrm = normpart(I,J)
        pol2 = expand(pol)
        p2c = pol2.coefficients(sparse=False)

        pol3 = nrm*pol2
        p3c = pol3.coefficients(sparse=False)
        for K in range(0,len(p3c)):
            p3c[K] = p3c[K].n()

        coefs.append(p3c)
        if(nmaxcoefs<len(p3c)):
            nmaxcoefs = len(p3c)
        
        #print("{:d}\t{:d}\t{}\t\t{}".format(I,J,nrm,p2c))
        
        s = ""
        s = s + "{:d}\t{:d}\t".format(I,J)
        for K in range(0,len(p3c)):
            nn = p3c[K]
            if(abs(nn)<1E-14):
                nn = 0
            s = s + "{:1.8g}\t".format(nn)
        s = s + "\n"
        print(s)
        lns.append(s)

##C style coefficient output

fout = "c_legendrecoefs2.c"
fp = open(fout,"w+")


pref1 = "static const int amslegendre_ndegree = {:d};\n".format(NN)
pref2 = "static const int amslegendre_maxporder = {:d};\n".format(nmaxcoefs-1)
pref = pref1 + pref2 + "static const float amslegendrecoefs[] = { \n"

for I in range(0,NN+1):
    for J in range(0,I+1):
        ind = I*(I+1)/2+J
        for K in range(0,nmaxcoefs):
            if(K<len(coefs[ind])):
                nn = coefs[ind][K]
                if(abs(nn)<1E-14):
                    nn = 0
            else:
                nn = 0
            
            if(ind<len(coefs)-1):
                if(K<nmaxcoefs-1):
                    pref = pref + "{:1.8g},".format(nn)
                else:
                    pref = pref + "{:1.8g},\n".format(nn)
            else:
                if(K<nmaxcoefs-1):
                    pref = pref + "{:1.8g},".format(nn)
                else:
                    pref = pref + "{:1.8g}".format(nn) + "\n};\n\n"

fp.writelines(pref)
fp.close()