from Numeric import *
import LinearAlgebra

par = [4.88120712, -0.073295799999999994, -0.088838760000000003, 0.30230843000000002, -0.99391039999999997, 0.99581405999999995, -0.64791631000000005, 0.39092521000000002, -0.37218462000000002, 0.22443650000000001, -0.11762777000000001, 0.089297920000000003, -0.035368660000000003, 0.01604531, -0.027709569999999999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.33862102999999999, -0.31581143, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

def Expm(alpha,x):
    rCut = 4.85
    withincutoffradius = less_equal(x,rCut)
    return (exp(-alpha*x)-(alpha*(rCut-x)+1)*exp(-alpha*rCut))*withincutoffradius

def Density(lat,alpha):
    return 8*Expm(alpha,lat*sqrt(3.0)/2.0)+6*Expm(alpha,lat)+12*Expm(alpha,lat*sqrt(2.0))

def LatticeConstant(rho,c):
    return c[0] + c[1]*log(rho)+c[2]*rho*log(rho)+c[3]*rho**2*log(rho)+c[4]*rho**3*log(rho) + rho*c[5] + 1/rho*c[6] + 1/rho**2*c[7]
    
alpha = 4.88120712
Y = []
A = []
norm = Density(3.1872,alpha)
for lat in arange(2.9,3.51,0.1):
    Y.append(lat)
    rho = Density(lat,alpha)/norm
    A.append([])
    for i in range(8):
        c = zeros((8,))
        c[i] = 1
        A[-1].append(LatticeConstant(rho,c))
for lat in arange(3.12,3.23,0.01):
    Y.append(lat)
    rho = Density(lat,alpha)/norm
    A.append([])
    for i in range(8):
        c = zeros((8,))
        c[i] = 1
        A[-1].append(LatticeConstant(rho,c))
Y = array(Y)
A = array(A)

print Y
print A
U,w,V = LinearAlgebra.singular_value_decomposition(A)