#!/usr/bin/python

#from Numeric import *
#from LinearAlgebra import *
from numpy.oldnumeric import *
from numpy.oldnumeric.linear_algebra import *
import sys

def read_file(name):
  old=open(name,'r')
  xvect=[]
  yvect=[]
  for line in old.readlines():
    prefield=line.split()
    xvect.append(float(prefield[0]))
    yvect.append(float(prefield[1]))
  old.close()
  return xvect,yvect

def energy_set(x,e0,v0,b0,db0):
  energy=zeros(len(x),Float)
  for i in range(len(x)):
    energy[i]=energy_point(x[i],e0,v0,b0,db0)
  return energy
  
def energy_point(x,e0,v0,b0,db0):
  energy=e0+b0*x/db0*(((v0/x)**db0)/(db0-1.)+1.)-b0*v0/(db0-1.)
  return energy

def pressure_point(x,v0,b0,db0):
  pressure=b0/db0*((v0/x)**db0-1.)
  return pressure

  
def residual_vector(x,y,e0,v0,b0,db0):
  y2=energy_set(x,e0,v0,b0,db0)
  res=y-y2
  return res

def jacobi_mat(x,y,e0,v0,b0,db):
  bmat=zeros((len(x),4),Float)
  for i in range(len(x)):
    dere=deriv_e(x[i],e0,v0,b0,db0)
    derv0=deriv_v0(x[i],e0,v0,b0,db0)
    derb0=deriv_b0(x[i],e0,v0,b0,db0)
    derdb0=deriv_db0(x[i],e0,v0,b0,db0)
    bmat[i,0]=dere
    bmat[i,1]=derv0
    bmat[i,2]=derb0
    bmat[i,3]=derdb0    
  return bmat
    
def jacobi_mat_num(x,y,e0,v0,b0,db):
  step=0.001
  bmatP=zeros((len(x),4),Float)
  bmatM=zeros((len(x),4),Float)
  #c forward
  for i in range(len(x)):
    bmatP[i,0]=energy_point(x[i],e0+step,v0,b0,db0)
    bmatP[i,1]=energy_point(x[i],e0,v0+step,b0,db0)
    bmatP[i,2]=energy_point(x[i],e0,v0,b0+step,db0)
    bmatP[i,3]=energy_point(x[i],e0,v0,b0,db0+step)
    
  #c backward
  for i in range(len(x)):
    bmatM[i,0]=energy_point(x[i],e0-step,v0,b0,db0)
    bmatM[i,1]=energy_point(x[i],e0,v0-step,b0,db0)
    bmatM[i,2]=energy_point(x[i],e0,v0,b0-step,db0)
    bmatM[i,3]=energy_point(x[i],e0,v0,b0,db0-step)
  
  bmat=(bmatP-bmatM)/2/step  
  return bmat 
    
def deriv_e(x,e0,v0,b0,db0):
  dere=1.
  return dere
 
def deriv_b0(x,e0,v0,b0,db0): 
  derb0=x/db0*(((v0/x)**db0)/(db0-1.)+1.)-v0/(db0-1.)
  return derb0
  
def deriv_v0(x,e0,v0,b0,db0): 
  derv0=b0/(db0-1.)*((v0/x)**(db0-1.)-1.)
  return derv0
  
def deriv_db0(x,e0,v0,b0,db0): 
  derdb0=b0*v0/(db0-1.)**2
  derdb0+=-b0*x/db0**2*(((v0/x)**db0)/(db0-1.)+1.)
  derdb0+=b0*x/db0*((log(v0/x)*(v0/x)**db0)/(db0-1.)-((v0/x)**db0)/(db0-1.)**2)
  return derdb0


#c setting 
step=0.5
nmax=200
#tolerance=1.e-5


name=sys.argv[1]
xvect,yvect=read_file(name)

#c initial estimate of fitting parameters 
e0=min(yvect)
indx=yvect.index(e0)
v0=xvect[indx]
b0=0.1
db0=3.5

xvect=array(xvect)
yvect=array(yvect)

ksi=array([e0,v0,b0,db0])

res=residual_vector(xvect,yvect,ksi[0],ksi[1],ksi[2],ksi[3])

i=1
while i<nmax:
  #bmat=jacobi_mat_num(xvect,yvect,ksi[0],ksi[1],ksi[2],ksi[3])
  bmat=jacobi_mat(xvect,yvect,ksi[0],ksi[1],ksi[2],ksi[3])
  gmat=dot(transpose(bmat),bmat)
  ginv=inverse(gmat)
  amat=dot(ginv,transpose(bmat))
  dksi=dot(amat,res)
  ksi+=step*dksi
  res=residual_vector(xvect,yvect,ksi[0],ksi[1],ksi[2],ksi[3])
  err=sum(res**2)
  i=i+1

for i in range(len(xvect)):
  pressure=pressure_point(xvect[i],ksi[1],ksi[2],ksi[3])
  energy=energy_point(xvect[i],ksi[0],ksi[1],ksi[2],ksi[3])
  #print xvect[i],energy,pressure*160.217657*10

print nmax,"iterations performed"
print "Ch-square:",err
print '---------'

print ""
print "\t","E0(eV):","\t",ksi[0]
print "\t","d0(A):","\t",ksi[1]