test_leibnitz

using leibnitz

x = randn()
@test  cos''(x)== -cos(x)

function rotmat(v)
   a = sqrt(dot(v,v)) # norm() does not take Fado's
   if a < 1e-9
      return eye(3) + [0. -v[3] v[2];v[3] 0. -v[1];-v[2] v[1] 0.]
   else
      c = cos(a)
      s = sin(a)
      v = v/a
      return (1.-c)*v*v'+c*eye(3) + s*[0. -v[3] v[2];v[3] 0. -v[1];-v[2] v[1] 0.]
   end
end
v      = [1.,1.,1.]
m_v    = derive(rotmat,v)
m_v1   = reshape(m_v[1,:,:],3,3)
Dv     = [1e-9,0.,0.]
m_vcmp = (rotmat(v+Dv)-rotmat(v-Dv))/2e-9
@test maximum(abs(m_v1-m_vcmp)) < 1e-7

function U(x,t)
   L0 = 0.1
   EA = 2.1e11*0.01^2
   a  = 1e-4
   Dx = x[1:2]-x[3:4]
   L  = sqrt(dot(Dx,Dx))
   DL = L-(L0*(1+a*(t)))
   return 1/2*EA*DL^2
end
R(x,t) = -derive(z->U(z,t),x)
K(x,t) = -derive(z->R(z,t),x)
x      = [.5001/5,.0/5,0.,0.]
t      = 0.
Kfasit = [2.1e7 0.0 -2.1e7 0.0;0.0 4199.16 -0.0 -4199.16;-2.1e7 0.0 2.1e7 0.0;0.0 -4199.16 -0.0 4199.16]
@test_approx_eq_eps R(x,t)  [-420.0,-0.0,420.0,-0.0] 1e-7
@test_approx_eq_eps K(x,t)  Kfasit 1e0

a           = derivatives(sin,0.,5)
@test a == [0., 1., 0., -1., 0., 1.]

f(x,y)    = 3*x^3*y+y^2+5*x+2
f_x(x,y)  = 9*x^2*y+5
f_xx(x,y) = 18*x*y
f_y(x,y)  = 3*x^3+2y
f_xy(x,y) = 9*x^2
f_xxy(x,y)= 18*x
f_yy(x,y) = 2
f_xyy(x,y)= 0.
f_xxyy(x,y)= 0.
d = derivatives(f,2.,3.)
@test d[1,1] == f(2.,3.)
@test d[2,1] == f_x(2.,3.)
@test d[1,2] == f_y(2.,3.)
@test d[2,2] == f_xy(2.,3.)
d = derivatives(f,2.,3.,2,2)
@test d[1,1] == f(2.,3.)
@test d[2,1] == f_x(2.,3.)
@test d[3,1] == f_xx(2.,3.)
@test d[1,2] == f_y(2.,3.)
@test d[2,2] == f_xy(2.,3.)
@test d[3,2] == f_xxy(2.,3.)
@test d[1,3] == f_yy(2.,3.)
@test d[2,3] == f_xyy(2.,3.)
@test d[3,3] == f_xxyy(2.,3.)

t=Taylor(cos,0.,2)
@test t(1.)==0.5