#!/usr/bin/python3

import itertools
import sympy

from sympy import sin, cos, acos, atan2

# coords
x, y, z, theta, phi = sympy.symbols("x y z theta phi", real = True)
r = sympy.symbols("r", positive = True)

coords_cart      = [x, y, z]
coords_spherical = [r, theta, phi]

X = r * sin(theta) * cos(phi)
Y = r * sin(theta) * sin(phi)
Z = r * cos(theta)

coords_cart_of_spherical = [X, Y, Z]
dcart_dspherical = sympy.derive_by_array(coords_cart_of_spherical, coords_spherical)

R     = (x**2 + y**2 + z**2) ** (1./2)
Theta = acos(z / R)
Phi   = atan2(y, x)

coords_spherical_of_cart = [R, Theta, Phi]
dspherical_dcart = sympy.derive_by_array(coords_spherical_of_cart, coords_cart)

# metric
metric_cart   = sympy.MutableDenseNDimArray.zeros(3, 3)
metric_cart_u = sympy.MutableDenseNDimArray.zeros(3, 3)

for i, j in itertools.product(range(3), repeat = 2):
    metric_cart[i, j]   = sympy.Function('g_%s%s'  % (str(coords_cart[i]), coords_cart[j]))(x, y, z)
    metric_cart_u[i, j] = sympy.Function('gu_%s%s' % (str(coords_cart[i]), coords_cart[j]))(x, y, z)

metric_spherical   = sympy.MutableDenseNDimArray.zeros(3, 3)
metric_spherical_u = sympy.MutableDenseNDimArray.zeros(3, 3)
for i, j in itertools.product(range(3), repeat = 2):
    val   = 0
    val_u = 0
    for k, l in itertools.product(range(3), repeat = 2):
        val   += dcart_dspherical[i, k] * dcart_dspherical[j, l] * metric_cart[k, l]
        val_u += dspherical_dcart[k, i] * dcart_dspherical[l, j] * metric_cart_u[k, l]

    metric_spherical[i, j]   = val
    metric_spherical_u[i, j] = val_u

print(metric_spherical)
#dmetric_spherical =

Gamma = sympy.MutableDenseNDimArray.zeros(3, 3, 3)
