def reference_metric(SymPySimplifyExpressions=True):
    global f0_of_xx0_funcform, f1_of_xx1_funcform, f2_of_xx0_xx1_funcform, f3_of_xx0_funcform, f4_of_xx2_funcform
    global f0_of_xx0, f1_of_xx1, f2_of_xx1, f2_of_xx0_xx1, f3_of_xx0, f4_of_xx2
    f0_of_xx0_funcform = sp.Function('f0_of_xx0_funcform')(xx[0])
    f1_of_xx1_funcform = sp.Function('f1_of_xx1_funcform')(xx[1])
    f2_of_xx0_xx1_funcform = sp.Function('f2_of_xx0_xx1_funcform')(xx[0], xx[1])
    f3_of_xx0_funcform = sp.Function('f3_of_xx0_funcform')(xx[0])
    f4_of_xx2_funcform = sp.Function('f4_of_xx2_funcform')(xx[2])
    (f0_of_xx0, f1_of_xx1, f2_of_xx1, f2_of_xx0_xx1, f3_of_xx0, f4_of_xx2) = par.Cparameters('REAL', thismodule, ['f0_of_xx0', 'f1_of_xx1', 'f2_of_xx1', 'f2_of_xx0_xx1', 'f3_of_xx0', 'f4_of_xx2'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f0_of_xx0__D0', 'f0_of_xx0__DD00', 'f0_of_xx0__DDD000'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f1_of_xx1__D1', 'f1_of_xx1__DD11', 'f1_of_xx1__DDD111'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f2_of_xx1__D1', 'f2_of_xx1__DD11', 'f2_of_xx1__DDD111'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f2_of_xx0_xx1__D0', 'f2_of_xx0_xx1__D1', 'f2_of_xx0_xx1__DD00', 'f2_of_xx0_xx1__DD11'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f3_of_xx0__D0', 'f3_of_xx0__DD00'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f4_of_xx2__D2', 'f4_of_xx2__DD22'], 1e+300)
    global have_already_called_reference_metric_function
    have_already_called_reference_metric_function = True
    CoordSystem = par.parval_from_str('reference_metric::CoordSystem')
    (M_PI, M_SQRT1_2) = par.Cparameters('#define', thismodule, ['M_PI', 'M_SQRT1_2'], '')
    global xxmin
    global xxmax
    global UnitVectors
    UnitVectors = ixp.zerorank2(DIM=3)
    if (CoordSystem in ('Spherical', 'SinhSpherical', 'SinhSphericalv2')):
        xx[0] = sp.symbols('xx0', real=True)
        xx[1] = sp.symbols('xx1', real=True)
        if (CoordSystem == 'Spherical'):
            RMAX = par.Cparameters('REAL', thismodule, ['RMAX'], 10.0)
            xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
            xxmax = [RMAX, M_PI, M_PI]
            r = xx[0]
            th = xx[1]
            ph = xx[2]
            Cart_to_xx[0] = sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))
            Cart_to_xx[1] = sp.acos((Cartz / Cart_to_xx[0]))
            Cart_to_xx[2] = sp.atan2(Carty, Cartx)
        elif (CoordSystem == 'SinhSpherical'):
            xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
            xxmax = [sp.sympify(1), M_PI, M_PI]
            (AMPL, SINHW) = par.Cparameters('REAL', thismodule, ['AMPL', 'SINHW'], [10.0, 0.2])
            r = ((AMPL * (sp.exp((xx[0] / SINHW)) - sp.exp(((- xx[0]) / SINHW)))) / (sp.exp((1 / SINHW)) - sp.exp(((- 1) / SINHW))))
            th = xx[1]
            ph = xx[2]
            Cart_to_xx[0] = (SINHW * sp.asinh(((sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2))) * sp.sinh((1 / SINHW))) / AMPL)))
            Cart_to_xx[1] = sp.acos((Cartz / sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))))
            Cart_to_xx[2] = sp.atan2(Carty, Cartx)
        elif (CoordSystem == 'SinhSphericalv2'):
            xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
            xxmax = [sp.sympify(1), M_PI, M_PI]
            (AMPL, SINHW) = par.Cparameters('REAL', thismodule, ['AMPL', 'SINHW'], [10.0, 0.2])
            const_dr = par.Cparameters('REAL', thismodule, ['const_dr'], 0.0625)
            r = (AMPL * ((const_dr * xx[0]) + ((sp.exp((xx[0] / SINHW)) - sp.exp(((- xx[0]) / SINHW))) / (sp.exp((1 / SINHW)) - sp.exp(((- 1) / SINHW))))))
            th = xx[1]
            ph = xx[2]
        xxSph[0] = r
        xxSph[1] = th
        xxSph[2] = ph
        xx_to_Cart[0] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.cos(xxSph[2]))
        xx_to_Cart[1] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.sin(xxSph[2]))
        xx_to_Cart[2] = (xxSph[0] * sp.cos(xxSph[1]))
        scalefactor_orthog[0] = sp.diff(xxSph[0], xx[0])
        scalefactor_orthog[1] = xxSph[0]
        scalefactor_orthog[2] = (xxSph[0] * sp.sin(xxSph[1]))
        f0_of_xx0 = xxSph[0]
        f1_of_xx1 = sp.sin(xxSph[1])
        scalefactor_orthog_funcform[0] = sp.diff(f0_of_xx0_funcform, xx[0])
        scalefactor_orthog_funcform[1] = f0_of_xx0_funcform
        scalefactor_orthog_funcform[2] = (f0_of_xx0_funcform * f1_of_xx1_funcform)
        UnitVectors = [[(sp.sin(xxSph[1]) * sp.cos(xxSph[2])), (sp.sin(xxSph[1]) * sp.sin(xxSph[2])), sp.cos(xxSph[1])], [(sp.cos(xxSph[1]) * sp.cos(xxSph[2])), (sp.cos(xxSph[1]) * sp.sin(xxSph[2])), (- sp.sin(xxSph[1]))], [(- sp.sin(xxSph[2])), sp.cos(xxSph[2]), sp.sympify(0)]]
    elif (CoordSystem in ('NobleSphericalThetaOptionOne', 'NobleSphericalThetaOptionTwo')):
        (Rin, R0) = par.Cparameters('REAL', thismodule, ['Rin', 'R0'], [1.08986052555408, 0.0])
        x0beg = sp.log((Rin - R0))
        xx[0] = sp.symbols('xx0', real=True)
        r = (R0 + sp.exp((x0beg + xx[0])))
        (th_c, xi, x1beg) = par.Cparameters('REAL', thismodule, ['th_c', 'xi', 'x1beg'], [0.053407075111026485, 0.25, 0.0])
        xx[1] = sp.symbols('xx1', real=True)
        x1j = (x1beg + xx[1])
        if (CoordSystem == 'NobleSphericalThetaOptionOne'):
            th = ((th_c + ((M_PI - (2 * th_c)) * x1j)) + (xi * sp.sin(((2 * M_PI) * x1j))))
        elif (CoordSystem == 'NobleSphericalThetaOptionTwo'):
            x1_n_exponent = par.Cparameters('REAL', thismodule, ['x1_n_exponent'], 9.0)
            th = ((M_PI / 2) * ((1 + ((1 - xi) * ((2 * x1j) - 1))) + ((xi - ((2 * th_c) / M_PI)) * (((2 * x1j) - 1) ** x1_n_exponent))))
        xx[2] = sp.symbols('xx2', real=True)
        ph = xx[2]
        xxSph[0] = r
        xxSph[1] = th
        xxSph[2] = ph
        xx_to_Cart[0] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.cos(xxSph[2]))
        xx_to_Cart[1] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.sin(xxSph[2]))
        xx_to_Cart[2] = (xxSph[0] * sp.cos(xxSph[1]))
        scalefactor_orthog[0] = sp.diff(xxSph[0], xx[0])
        scalefactor_orthog[1] = xxSph[0]
        scalefactor_orthog[2] = (xxSph[0] * sp.sin(xxSph[1]))
        UnitVectors = [[(sp.sin(xxSph[1]) * sp.cos(xxSph[2])), (sp.sin(xxSph[1]) * sp.sin(xxSph[2])), sp.cos(xxSph[1])], [(sp.cos(xxSph[1]) * sp.cos(xxSph[2])), (sp.cos(xxSph[1]) * sp.sin(xxSph[2])), (- sp.sin(xxSph[1]))], [(- sp.sin(xxSph[2])), sp.cos(xxSph[2]), sp.sympify(0)]]
    elif (CoordSystem in ('Cylindrical', 'SinhCylindrical', 'SinhCylindricalv2')):
        xx[0] = sp.symbols('xx0', real=True)
        if (CoordSystem == 'Cylindrical'):
            (RHOMAX, ZMIN, ZMAX) = par.Cparameters('REAL', thismodule, ['RHOMAX', 'ZMIN', 'ZMAX'], [10.0, (- 10.0), 10.0])
            xxmin = [sp.sympify(0), (- M_PI), ZMIN]
            xxmax = [RHOMAX, M_PI, ZMAX]
            RHOCYL = xx[0]
            PHICYL = xx[1]
            ZCYL = xx[2]
            Cart_to_xx[0] = sp.sqrt(((Cartx ** 2) + (Carty ** 2)))
            Cart_to_xx[1] = sp.atan2(Carty, Cartx)
            Cart_to_xx[2] = Cartz
        elif (CoordSystem == 'SinhCylindrical'):
            xxmin = [sp.sympify(0), (- M_PI), sp.sympify((- 1))]
            xxmax = [sp.sympify(1), M_PI, sp.sympify((+ 1))]
            (AMPLRHO, SINHWRHO, AMPLZ, SINHWZ) = par.Cparameters('REAL', thismodule, ['AMPLRHO', 'SINHWRHO', 'AMPLZ', 'SINHWZ'], [10.0, 0.2, 10.0, 0.2])
            RHOCYL = ((AMPLRHO * (sp.exp((xx[0] / SINHWRHO)) - sp.exp(((- xx[0]) / SINHWRHO)))) / (sp.exp((1 / SINHWRHO)) - sp.exp(((- 1) / SINHWRHO))))
            PHICYL = xx[1]
            ZCYL = ((AMPLZ * (sp.exp((xx[2] / SINHWZ)) - sp.exp(((- xx[2]) / SINHWZ)))) / (sp.exp((1 / SINHWZ)) - sp.exp(((- 1) / SINHWZ))))
            Cart_to_xx[0] = (SINHWRHO * sp.asinh(((sp.sqrt(((Cartx ** 2) + (Carty ** 2))) * sp.sinh((1 / SINHWRHO))) / AMPLRHO)))
            Cart_to_xx[1] = sp.atan2(Carty, Cartx)
            Cart_to_xx[2] = (SINHWZ * sp.asinh(((Cartz * sp.sinh((1 / SINHWZ))) / AMPLZ)))
        elif (CoordSystem == 'SinhCylindricalv2'):
            xxmin = [sp.sympify(0), (- M_PI), sp.sympify((- 1))]
            xxmax = [sp.sympify(1), M_PI, sp.sympify((+ 1))]
            (AMPLRHO, SINHWRHO, AMPLZ, SINHWZ) = par.Cparameters('REAL', thismodule, ['AMPLRHO', 'SINHWRHO', 'AMPLZ', 'SINHWZ'], [10.0, 0.2, 10.0, 0.2])
            (const_drho, const_dz) = par.Cparameters('REAL', thismodule, ['const_drho', 'const_dz'], [0.0625, 0.0625])
            RHOCYL = (AMPLRHO * ((const_drho * xx[0]) + ((sp.exp((xx[0] / SINHWRHO)) - sp.exp(((- xx[0]) / SINHWRHO))) / (sp.exp((1 / SINHWRHO)) - sp.exp(((- 1) / SINHWRHO))))))
            PHICYL = xx[1]
            ZCYL = (AMPLZ * ((const_dz * xx[2]) + ((sp.exp((xx[2] / SINHWZ)) - sp.exp(((- xx[2]) / SINHWZ))) / (sp.exp((1 / SINHWZ)) - sp.exp(((- 1) / SINHWZ))))))
        xx_to_Cart[0] = (RHOCYL * sp.cos(PHICYL))
        xx_to_Cart[1] = (RHOCYL * sp.sin(PHICYL))
        xx_to_Cart[2] = ZCYL
        xxSph[0] = sp.sqrt(((RHOCYL ** 2) + (ZCYL ** 2)))
        xxSph[1] = sp.acos((ZCYL / xxSph[0]))
        xxSph[2] = PHICYL
        scalefactor_orthog[0] = sp.diff(RHOCYL, xx[0])
        scalefactor_orthog[1] = RHOCYL
        scalefactor_orthog[2] = sp.diff(ZCYL, xx[2])
        f0_of_xx0 = RHOCYL
        f4_of_xx2 = sp.diff(ZCYL, xx[2])
        scalefactor_orthog_funcform[0] = sp.diff(f0_of_xx0_funcform, xx[0])
        scalefactor_orthog_funcform[1] = f0_of_xx0_funcform
        scalefactor_orthog_funcform[2] = f4_of_xx2_funcform
        UnitVectors = [[sp.cos(PHICYL), sp.sin(PHICYL), sp.sympify(0)], [(- sp.sin(PHICYL)), sp.cos(PHICYL), sp.sympify(0)], [sp.sympify(0), sp.sympify(0), sp.sympify(1)]]
    elif (CoordSystem in ('SymTP', 'SinhSymTP')):
        (bScale, SINHWAA, AMAX) = par.Cparameters('REAL', thismodule, ['bScale', 'SINHWAA', 'AMAX'], [0.5, 0.2, 10.0])
        (xx[0], xx[1]) = sp.symbols('xx0 xx1', real=True)
        xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
        xxmax = [AMAX, M_PI, M_PI]
        AA = xx[0]
        if (CoordSystem == 'SinhSymTP'):
            xxmax[0] = sp.sympify(1)
            AA = ((AMAX * (sp.exp((xx[0] / SINHWAA)) - sp.exp(((- xx[0]) / SINHWAA)))) / (sp.exp((1 / SINHWAA)) - sp.exp(((- 1) / SINHWAA))))
        var1 = sp.sqrt(((AA ** 2) + ((bScale * sp.sin(xx[1])) ** 2)))
        var2 = sp.sqrt(((AA ** 2) + (bScale ** 2)))
        RHOSYMTP = (AA * sp.sin(xx[1]))
        PHSYMTP = xx[2]
        ZSYMTP = (var2 * sp.cos(xx[1]))
        xx_to_Cart[0] = ((AA * sp.sin(xx[1])) * sp.cos(xx[2]))
        xx_to_Cart[1] = ((AA * sp.sin(xx[1])) * sp.sin(xx[2]))
        xx_to_Cart[2] = ZSYMTP
        xxSph[0] = sp.sqrt(((RHOSYMTP ** 2) + (ZSYMTP ** 2)))
        xxSph[1] = sp.acos((ZSYMTP / xxSph[0]))
        xxSph[2] = PHSYMTP
        if (CoordSystem == 'SymTP'):
            rSph = sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))
            thSph = sp.acos((Cartz / rSph))
            phSph = sp.atan2(Carty, Cartx)
            Cart_to_xx[0] = (sp.sqrt((((- (bScale ** 2)) + (rSph ** 2)) + sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))))) * M_SQRT1_2)
            Cart_to_xx[1] = sp.acos((sp.sign(Cartz) * (sp.sqrt(((1 + ((rSph ** 2) / (bScale ** 2))) - (sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))) / (bScale ** 2)))) * M_SQRT1_2)))
            Cart_to_xx[2] = phSph
        elif (CoordSystem == 'SinhSymTP'):
            rSph = sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))
            thSph = sp.acos((Cartz / rSph))
            phSph = sp.atan2(Carty, Cartx)
            Cart_to_xx[0] = (sp.sqrt((((- (bScale ** 2)) + (rSph ** 2)) + sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))))) * M_SQRT1_2)
            Cart_to_xx[1] = sp.acos((sp.sign(Cartz) * (sp.sqrt(((1 + ((rSph ** 2) / (bScale ** 2))) - (sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))) / (bScale ** 2)))) * M_SQRT1_2)))
            Cart_to_xx[2] = phSph
        scalefactor_orthog[0] = ((sp.diff(AA, xx[0]) * var1) / var2)
        scalefactor_orthog[1] = var1
        scalefactor_orthog[2] = (AA * sp.sin(xx[1]))
        f0_of_xx0 = AA
        f1_of_xx1 = sp.sin(xx[1])
        f2_of_xx0_xx1 = var1
        f3_of_xx0 = var2
        scalefactor_orthog_funcform[0] = ((sp.diff(f0_of_xx0_funcform, xx[0]) * f2_of_xx0_xx1_funcform) / f3_of_xx0_funcform)
        scalefactor_orthog_funcform[1] = f2_of_xx0_xx1_funcform
        scalefactor_orthog_funcform[2] = (f0_of_xx0_funcform * f1_of_xx1_funcform)
        UnitVectors = [[(((sp.sin(xx[1]) * sp.cos(xx[2])) * var2) / var1), (((sp.sin(xx[1]) * sp.sin(xx[2])) * var2) / var1), ((AA * sp.cos(xx[1])) / var1)], [(((AA * sp.cos(xx[1])) * sp.cos(xx[2])) / var1), (((AA * sp.cos(xx[1])) * sp.sin(xx[2])) / var1), (((- sp.sin(xx[1])) * var2) / var1)], [(- sp.sin(xx[2])), sp.cos(xx[2]), sp.sympify(0)]]
    elif (CoordSystem == 'Cartesian'):
        par.Cparameters('REAL', thismodule, ['xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax'], [(- 10.0), 10.0, (- 10.0), 10.0, (- 10.0), 10.0])
        xxmin = ['xmin', 'ymin', 'zmin']
        xxmax = ['xmax', 'ymax', 'zmax']
        xx_to_Cart[0] = xx[0]
        xx_to_Cart[1] = xx[1]
        xx_to_Cart[2] = xx[2]
        xxSph[0] = sp.sqrt((((xx[0] ** 2) + (xx[1] ** 2)) + (xx[2] ** 2)))
        xxSph[1] = sp.acos((xx[2] / xxSph[0]))
        xxSph[2] = sp.atan2(xx[1], xx[0])
        Cart_to_xx[0] = Cartx
        Cart_to_xx[1] = Carty
        Cart_to_xx[2] = Cartz
        scalefactor_orthog[0] = sp.sympify(1)
        scalefactor_orthog[1] = sp.sympify(1)
        scalefactor_orthog[2] = sp.sympify(1)
        scalefactor_orthog_funcform[0] = sp.sympify(1)
        scalefactor_orthog_funcform[1] = sp.sympify(1)
        scalefactor_orthog_funcform[2] = sp.sympify(1)
        UnitVectors = [[sp.sympify(1), sp.sympify(0), sp.sympify(0)], [sp.sympify(0), sp.sympify(1), sp.sympify(0)], [sp.sympify(0), sp.sympify(0), sp.sympify(1)]]
    elif (CoordSystem == 'SinhCartesian'):
        xxmin = [sp.sympify((- 1)), sp.sympify((- 1)), sp.sympify((- 1))]
        xxmax = [sp.sympify((+ 1)), sp.sympify((+ 1)), sp.sympify((+ 1))]
        (AMPLXYZ, SINHWXYZ) = par.Cparameters('REAL', thismodule, ['AMPLXYZ', 'SINHWXYZ'], [10.0, 0.2])
        for ii in [0, 1, 2]:
            xx_to_Cart[ii] = ((AMPLXYZ * (sp.exp((xx[ii] / SINHWXYZ)) - sp.exp(((- xx[ii]) / SINHWXYZ)))) / (sp.exp((1 / SINHWXYZ)) - sp.exp(((- 1) / SINHWXYZ))))
        xxSph[0] = sp.sqrt((((xx_to_Cart[0] ** 2) + (xx_to_Cart[1] ** 2)) + (xx_to_Cart[2] ** 2)))
        xxSph[1] = sp.acos((xx_to_Cart[2] / xxSph[0]))
        xxSph[2] = sp.atan2(xx_to_Cart[1], xx_to_Cart[0])
        Cart_to_xx[0] = (SINHWXYZ * sp.asinh(((Cartx * sp.sinh((1 / SINHWXYZ))) / AMPLXYZ)))
        Cart_to_xx[1] = (SINHWXYZ * sp.asinh(((Carty * sp.sinh((1 / SINHWXYZ))) / AMPLXYZ)))
        Cart_to_xx[2] = (SINHWXYZ * sp.asinh(((Cartz * sp.sinh((1 / SINHWXYZ))) / AMPLXYZ)))
        scalefactor_orthog[0] = sp.diff(xx_to_Cart[0], xx[0])
        scalefactor_orthog[1] = sp.diff(xx_to_Cart[1], xx[1])
        scalefactor_orthog[2] = sp.diff(xx_to_Cart[2], xx[2])
        f0_of_xx0 = sp.diff(xx_to_Cart[0], xx[0])
        f1_of_xx1 = sp.diff(xx_to_Cart[1], xx[1])
        f4_of_xx2 = sp.diff(xx_to_Cart[2], xx[2])
        scalefactor_orthog_funcform[0] = f0_of_xx0_funcform
        scalefactor_orthog_funcform[1] = f1_of_xx1_funcform
        scalefactor_orthog_funcform[2] = f4_of_xx2_funcform
        UnitVectors = [[sp.sympify(1), sp.sympify(0), sp.sympify(0)], [sp.sympify(0), sp.sympify(1), sp.sympify(0)], [sp.sympify(0), sp.sympify(0), sp.sympify(1)]]
    else:
        print((('CoordSystem == ' + CoordSystem) + ' is not supported.'))
        sys.exit(1)
    ref_metric__hatted_quantities(SymPySimplifyExpressions)
def reference_metric(SymPySimplifyExpressions=True):
    global f0_of_xx0_funcform, f1_of_xx1_funcform, f2_of_xx0_xx1_funcform, f3_of_xx0_funcform, f4_of_xx2_funcform
    global f0_of_xx0, f1_of_xx1, f2_of_xx1, f2_of_xx0_xx1, f3_of_xx0, f4_of_xx2
    f0_of_xx0_funcform = sp.Function('f0_of_xx0_funcform')(xx[0])
    f1_of_xx1_funcform = sp.Function('f1_of_xx1_funcform')(xx[1])
    f2_of_xx0_xx1_funcform = sp.Function('f2_of_xx0_xx1_funcform')(xx[0], xx[1])
    f3_of_xx0_funcform = sp.Function('f3_of_xx0_funcform')(xx[0])
    f4_of_xx2_funcform = sp.Function('f4_of_xx2_funcform')(xx[2])
    (f0_of_xx0, f1_of_xx1, f2_of_xx1, f2_of_xx0_xx1, f3_of_xx0, f4_of_xx2) = par.Cparameters('REAL', thismodule, ['f0_of_xx0', 'f1_of_xx1', 'f2_of_xx1', 'f2_of_xx0_xx1', 'f3_of_xx0', 'f4_of_xx2'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f0_of_xx0__D0', 'f0_of_xx0__DD00', 'f0_of_xx0__DDD000'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f1_of_xx1__D1', 'f1_of_xx1__DD11', 'f1_of_xx1__DDD111'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f2_of_xx1__D1', 'f2_of_xx1__DD11', 'f2_of_xx1__DDD111'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f2_of_xx0_xx1__D0', 'f2_of_xx0_xx1__D1', 'f2_of_xx0_xx1__DD00', 'f2_of_xx0_xx1__DD11'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f3_of_xx0__D0', 'f3_of_xx0__DD00'], 1e+300)
    par.Cparameters('REAL', thismodule, ['f4_of_xx2__D2', 'f4_of_xx2__DD22'], 1e+300)
    global have_already_called_reference_metric_function
    have_already_called_reference_metric_function = True
    CoordSystem = par.parval_from_str('reference_metric::CoordSystem')
    (M_PI, M_SQRT1_2) = par.Cparameters('#define', thismodule, ['M_PI', 'M_SQRT1_2'], )
    global xxmin
    global xxmax
    global UnitVectors
    UnitVectors = ixp.zerorank2(DIM=3)
    if (CoordSystem in ('Spherical', 'SinhSpherical', 'SinhSphericalv2')):
        xx[0] = sp.symbols('xx0', real=True)
        xx[1] = sp.symbols('xx1', real=True)
        if (CoordSystem == 'Spherical'):
            RMAX = par.Cparameters('REAL', thismodule, ['RMAX'], 10.0)
            xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
            xxmax = [RMAX, M_PI, M_PI]
            r = xx[0]
            th = xx[1]
            ph = xx[2]
            Cart_to_xx[0] = sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))
            Cart_to_xx[1] = sp.acos((Cartz / Cart_to_xx[0]))
            Cart_to_xx[2] = sp.atan2(Carty, Cartx)
        elif (CoordSystem == 'SinhSpherical'):
            xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
            xxmax = [sp.sympify(1), M_PI, M_PI]
            (AMPL, SINHW) = par.Cparameters('REAL', thismodule, ['AMPL', 'SINHW'], [10.0, 0.2])
            r = ((AMPL * (sp.exp((xx[0] / SINHW)) - sp.exp(((- xx[0]) / SINHW)))) / (sp.exp((1 / SINHW)) - sp.exp(((- 1) / SINHW))))
            th = xx[1]
            ph = xx[2]
            Cart_to_xx[0] = (SINHW * sp.asinh(((sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2))) * sp.sinh((1 / SINHW))) / AMPL)))
            Cart_to_xx[1] = sp.acos((Cartz / sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))))
            Cart_to_xx[2] = sp.atan2(Carty, Cartx)
        elif (CoordSystem == 'SinhSphericalv2'):
            xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
            xxmax = [sp.sympify(1), M_PI, M_PI]
            (AMPL, SINHW) = par.Cparameters('REAL', thismodule, ['AMPL', 'SINHW'], [10.0, 0.2])
            const_dr = par.Cparameters('REAL', thismodule, ['const_dr'], 0.0625)
            r = (AMPL * ((const_dr * xx[0]) + ((sp.exp((xx[0] / SINHW)) - sp.exp(((- xx[0]) / SINHW))) / (sp.exp((1 / SINHW)) - sp.exp(((- 1) / SINHW))))))
            th = xx[1]
            ph = xx[2]
        xxSph[0] = r
        xxSph[1] = th
        xxSph[2] = ph
        xx_to_Cart[0] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.cos(xxSph[2]))
        xx_to_Cart[1] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.sin(xxSph[2]))
        xx_to_Cart[2] = (xxSph[0] * sp.cos(xxSph[1]))
        scalefactor_orthog[0] = sp.diff(xxSph[0], xx[0])
        scalefactor_orthog[1] = xxSph[0]
        scalefactor_orthog[2] = (xxSph[0] * sp.sin(xxSph[1]))
        f0_of_xx0 = xxSph[0]
        f1_of_xx1 = sp.sin(xxSph[1])
        scalefactor_orthog_funcform[0] = sp.diff(f0_of_xx0_funcform, xx[0])
        scalefactor_orthog_funcform[1] = f0_of_xx0_funcform
        scalefactor_orthog_funcform[2] = (f0_of_xx0_funcform * f1_of_xx1_funcform)
        UnitVectors = [[(sp.sin(xxSph[1]) * sp.cos(xxSph[2])), (sp.sin(xxSph[1]) * sp.sin(xxSph[2])), sp.cos(xxSph[1])], [(sp.cos(xxSph[1]) * sp.cos(xxSph[2])), (sp.cos(xxSph[1]) * sp.sin(xxSph[2])), (- sp.sin(xxSph[1]))], [(- sp.sin(xxSph[2])), sp.cos(xxSph[2]), sp.sympify(0)]]
    elif (CoordSystem in ('NobleSphericalThetaOptionOne', 'NobleSphericalThetaOptionTwo')):
        (Rin, R0) = par.Cparameters('REAL', thismodule, ['Rin', 'R0'], [1.08986052555408, 0.0])
        x0beg = sp.log((Rin - R0))
        xx[0] = sp.symbols('xx0', real=True)
        r = (R0 + sp.exp((x0beg + xx[0])))
        (th_c, xi, x1beg) = par.Cparameters('REAL', thismodule, ['th_c', 'xi', 'x1beg'], [0.053407075111026485, 0.25, 0.0])
        xx[1] = sp.symbols('xx1', real=True)
        x1j = (x1beg + xx[1])
        if (CoordSystem == 'NobleSphericalThetaOptionOne'):
            th = ((th_c + ((M_PI - (2 * th_c)) * x1j)) + (xi * sp.sin(((2 * M_PI) * x1j))))
        elif (CoordSystem == 'NobleSphericalThetaOptionTwo'):
            x1_n_exponent = par.Cparameters('REAL', thismodule, ['x1_n_exponent'], 9.0)
            th = ((M_PI / 2) * ((1 + ((1 - xi) * ((2 * x1j) - 1))) + ((xi - ((2 * th_c) / M_PI)) * (((2 * x1j) - 1) ** x1_n_exponent))))
        xx[2] = sp.symbols('xx2', real=True)
        ph = xx[2]
        xxSph[0] = r
        xxSph[1] = th
        xxSph[2] = ph
        xx_to_Cart[0] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.cos(xxSph[2]))
        xx_to_Cart[1] = ((xxSph[0] * sp.sin(xxSph[1])) * sp.sin(xxSph[2]))
        xx_to_Cart[2] = (xxSph[0] * sp.cos(xxSph[1]))
        scalefactor_orthog[0] = sp.diff(xxSph[0], xx[0])
        scalefactor_orthog[1] = xxSph[0]
        scalefactor_orthog[2] = (xxSph[0] * sp.sin(xxSph[1]))
        UnitVectors = [[(sp.sin(xxSph[1]) * sp.cos(xxSph[2])), (sp.sin(xxSph[1]) * sp.sin(xxSph[2])), sp.cos(xxSph[1])], [(sp.cos(xxSph[1]) * sp.cos(xxSph[2])), (sp.cos(xxSph[1]) * sp.sin(xxSph[2])), (- sp.sin(xxSph[1]))], [(- sp.sin(xxSph[2])), sp.cos(xxSph[2]), sp.sympify(0)]]
    elif (CoordSystem in ('Cylindrical', 'SinhCylindrical', 'SinhCylindricalv2')):
        xx[0] = sp.symbols('xx0', real=True)
        if (CoordSystem == 'Cylindrical'):
            (RHOMAX, ZMIN, ZMAX) = par.Cparameters('REAL', thismodule, ['RHOMAX', 'ZMIN', 'ZMAX'], [10.0, (- 10.0), 10.0])
            xxmin = [sp.sympify(0), (- M_PI), ZMIN]
            xxmax = [RHOMAX, M_PI, ZMAX]
            RHOCYL = xx[0]
            PHICYL = xx[1]
            ZCYL = xx[2]
            Cart_to_xx[0] = sp.sqrt(((Cartx ** 2) + (Carty ** 2)))
            Cart_to_xx[1] = sp.atan2(Carty, Cartx)
            Cart_to_xx[2] = Cartz
        elif (CoordSystem == 'SinhCylindrical'):
            xxmin = [sp.sympify(0), (- M_PI), sp.sympify((- 1))]
            xxmax = [sp.sympify(1), M_PI, sp.sympify((+ 1))]
            (AMPLRHO, SINHWRHO, AMPLZ, SINHWZ) = par.Cparameters('REAL', thismodule, ['AMPLRHO', 'SINHWRHO', 'AMPLZ', 'SINHWZ'], [10.0, 0.2, 10.0, 0.2])
            RHOCYL = ((AMPLRHO * (sp.exp((xx[0] / SINHWRHO)) - sp.exp(((- xx[0]) / SINHWRHO)))) / (sp.exp((1 / SINHWRHO)) - sp.exp(((- 1) / SINHWRHO))))
            PHICYL = xx[1]
            ZCYL = ((AMPLZ * (sp.exp((xx[2] / SINHWZ)) - sp.exp(((- xx[2]) / SINHWZ)))) / (sp.exp((1 / SINHWZ)) - sp.exp(((- 1) / SINHWZ))))
            Cart_to_xx[0] = (SINHWRHO * sp.asinh(((sp.sqrt(((Cartx ** 2) + (Carty ** 2))) * sp.sinh((1 / SINHWRHO))) / AMPLRHO)))
            Cart_to_xx[1] = sp.atan2(Carty, Cartx)
            Cart_to_xx[2] = (SINHWZ * sp.asinh(((Cartz * sp.sinh((1 / SINHWZ))) / AMPLZ)))
        elif (CoordSystem == 'SinhCylindricalv2'):
            xxmin = [sp.sympify(0), (- M_PI), sp.sympify((- 1))]
            xxmax = [sp.sympify(1), M_PI, sp.sympify((+ 1))]
            (AMPLRHO, SINHWRHO, AMPLZ, SINHWZ) = par.Cparameters('REAL', thismodule, ['AMPLRHO', 'SINHWRHO', 'AMPLZ', 'SINHWZ'], [10.0, 0.2, 10.0, 0.2])
            (const_drho, const_dz) = par.Cparameters('REAL', thismodule, ['const_drho', 'const_dz'], [0.0625, 0.0625])
            RHOCYL = (AMPLRHO * ((const_drho * xx[0]) + ((sp.exp((xx[0] / SINHWRHO)) - sp.exp(((- xx[0]) / SINHWRHO))) / (sp.exp((1 / SINHWRHO)) - sp.exp(((- 1) / SINHWRHO))))))
            PHICYL = xx[1]
            ZCYL = (AMPLZ * ((const_dz * xx[2]) + ((sp.exp((xx[2] / SINHWZ)) - sp.exp(((- xx[2]) / SINHWZ))) / (sp.exp((1 / SINHWZ)) - sp.exp(((- 1) / SINHWZ))))))
        xx_to_Cart[0] = (RHOCYL * sp.cos(PHICYL))
        xx_to_Cart[1] = (RHOCYL * sp.sin(PHICYL))
        xx_to_Cart[2] = ZCYL
        xxSph[0] = sp.sqrt(((RHOCYL ** 2) + (ZCYL ** 2)))
        xxSph[1] = sp.acos((ZCYL / xxSph[0]))
        xxSph[2] = PHICYL
        scalefactor_orthog[0] = sp.diff(RHOCYL, xx[0])
        scalefactor_orthog[1] = RHOCYL
        scalefactor_orthog[2] = sp.diff(ZCYL, xx[2])
        f0_of_xx0 = RHOCYL
        f4_of_xx2 = sp.diff(ZCYL, xx[2])
        scalefactor_orthog_funcform[0] = sp.diff(f0_of_xx0_funcform, xx[0])
        scalefactor_orthog_funcform[1] = f0_of_xx0_funcform
        scalefactor_orthog_funcform[2] = f4_of_xx2_funcform
        UnitVectors = [[sp.cos(PHICYL), sp.sin(PHICYL), sp.sympify(0)], [(- sp.sin(PHICYL)), sp.cos(PHICYL), sp.sympify(0)], [sp.sympify(0), sp.sympify(0), sp.sympify(1)]]
    elif (CoordSystem in ('SymTP', 'SinhSymTP')):
        (bScale, SINHWAA, AMAX) = par.Cparameters('REAL', thismodule, ['bScale', 'SINHWAA', 'AMAX'], [0.5, 0.2, 10.0])
        (xx[0], xx[1]) = sp.symbols('xx0 xx1', real=True)
        xxmin = [sp.sympify(0), sp.sympify(0), (- M_PI)]
        xxmax = [AMAX, M_PI, M_PI]
        AA = xx[0]
        if (CoordSystem == 'SinhSymTP'):
            xxmax[0] = sp.sympify(1)
            AA = ((AMAX * (sp.exp((xx[0] / SINHWAA)) - sp.exp(((- xx[0]) / SINHWAA)))) / (sp.exp((1 / SINHWAA)) - sp.exp(((- 1) / SINHWAA))))
        var1 = sp.sqrt(((AA ** 2) + ((bScale * sp.sin(xx[1])) ** 2)))
        var2 = sp.sqrt(((AA ** 2) + (bScale ** 2)))
        RHOSYMTP = (AA * sp.sin(xx[1]))
        PHSYMTP = xx[2]
        ZSYMTP = (var2 * sp.cos(xx[1]))
        xx_to_Cart[0] = ((AA * sp.sin(xx[1])) * sp.cos(xx[2]))
        xx_to_Cart[1] = ((AA * sp.sin(xx[1])) * sp.sin(xx[2]))
        xx_to_Cart[2] = ZSYMTP
        xxSph[0] = sp.sqrt(((RHOSYMTP ** 2) + (ZSYMTP ** 2)))
        xxSph[1] = sp.acos((ZSYMTP / xxSph[0]))
        xxSph[2] = PHSYMTP
        if (CoordSystem == 'SymTP'):
            rSph = sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))
            thSph = sp.acos((Cartz / rSph))
            phSph = sp.atan2(Carty, Cartx)
            Cart_to_xx[0] = (sp.sqrt((((- (bScale ** 2)) + (rSph ** 2)) + sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))))) * M_SQRT1_2)
            Cart_to_xx[1] = sp.acos((sp.sign(Cartz) * (sp.sqrt(((1 + ((rSph ** 2) / (bScale ** 2))) - (sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))) / (bScale ** 2)))) * M_SQRT1_2)))
            Cart_to_xx[2] = phSph
        elif (CoordSystem == 'SinhSymTP'):
            rSph = sp.sqrt((((Cartx ** 2) + (Carty ** 2)) + (Cartz ** 2)))
            thSph = sp.acos((Cartz / rSph))
            phSph = sp.atan2(Carty, Cartx)
            Cart_to_xx[0] = (sp.sqrt((((- (bScale ** 2)) + (rSph ** 2)) + sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))))) * M_SQRT1_2)
            Cart_to_xx[1] = sp.acos((sp.sign(Cartz) * (sp.sqrt(((1 + ((rSph ** 2) / (bScale ** 2))) - (sp.sqrt(((((bScale ** 4) + ((2 * (bScale ** 2)) * (rSph ** 2))) + (rSph ** 4)) - (((4 * (bScale ** 2)) * (rSph ** 2)) * (sp.cos(thSph) ** 2)))) / (bScale ** 2)))) * M_SQRT1_2)))
            Cart_to_xx[2] = phSph
        scalefactor_orthog[0] = ((sp.diff(AA, xx[0]) * var1) / var2)
        scalefactor_orthog[1] = var1
        scalefactor_orthog[2] = (AA * sp.sin(xx[1]))
        f0_of_xx0 = AA
        f1_of_xx1 = sp.sin(xx[1])
        f2_of_xx0_xx1 = var1
        f3_of_xx0 = var2
        scalefactor_orthog_funcform[0] = ((sp.diff(f0_of_xx0_funcform, xx[0]) * f2_of_xx0_xx1_funcform) / f3_of_xx0_funcform)
        scalefactor_orthog_funcform[1] = f2_of_xx0_xx1_funcform
        scalefactor_orthog_funcform[2] = (f0_of_xx0_funcform * f1_of_xx1_funcform)
        UnitVectors = [[(((sp.sin(xx[1]) * sp.cos(xx[2])) * var2) / var1), (((sp.sin(xx[1]) * sp.sin(xx[2])) * var2) / var1), ((AA * sp.cos(xx[1])) / var1)], [(((AA * sp.cos(xx[1])) * sp.cos(xx[2])) / var1), (((AA * sp.cos(xx[1])) * sp.sin(xx[2])) / var1), (((- sp.sin(xx[1])) * var2) / var1)], [(- sp.sin(xx[2])), sp.cos(xx[2]), sp.sympify(0)]]
    elif (CoordSystem == 'Cartesian'):
        par.Cparameters('REAL', thismodule, ['xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax'], [(- 10.0), 10.0, (- 10.0), 10.0, (- 10.0), 10.0])
        xxmin = ['xmin', 'ymin', 'zmin']
        xxmax = ['xmax', 'ymax', 'zmax']
        xx_to_Cart[0] = xx[0]
        xx_to_Cart[1] = xx[1]
        xx_to_Cart[2] = xx[2]
        xxSph[0] = sp.sqrt((((xx[0] ** 2) + (xx[1] ** 2)) + (xx[2] ** 2)))
        xxSph[1] = sp.acos((xx[2] / xxSph[0]))
        xxSph[2] = sp.atan2(xx[1], xx[0])
        Cart_to_xx[0] = Cartx
        Cart_to_xx[1] = Carty
        Cart_to_xx[2] = Cartz
        scalefactor_orthog[0] = sp.sympify(1)
        scalefactor_orthog[1] = sp.sympify(1)
        scalefactor_orthog[2] = sp.sympify(1)
        scalefactor_orthog_funcform[0] = sp.sympify(1)
        scalefactor_orthog_funcform[1] = sp.sympify(1)
        scalefactor_orthog_funcform[2] = sp.sympify(1)
        UnitVectors = [[sp.sympify(1), sp.sympify(0), sp.sympify(0)], [sp.sympify(0), sp.sympify(1), sp.sympify(0)], [sp.sympify(0), sp.sympify(0), sp.sympify(1)]]
    elif (CoordSystem == 'SinhCartesian'):
        xxmin = [sp.sympify((- 1)), sp.sympify((- 1)), sp.sympify((- 1))]
        xxmax = [sp.sympify((+ 1)), sp.sympify((+ 1)), sp.sympify((+ 1))]
        (AMPLXYZ, SINHWXYZ) = par.Cparameters('REAL', thismodule, ['AMPLXYZ', 'SINHWXYZ'], [10.0, 0.2])
        for ii in [0, 1, 2]:
            xx_to_Cart[ii] = ((AMPLXYZ * (sp.exp((xx[ii] / SINHWXYZ)) - sp.exp(((- xx[ii]) / SINHWXYZ)))) / (sp.exp((1 / SINHWXYZ)) - sp.exp(((- 1) / SINHWXYZ))))
        xxSph[0] = sp.sqrt((((xx_to_Cart[0] ** 2) + (xx_to_Cart[1] ** 2)) + (xx_to_Cart[2] ** 2)))
        xxSph[1] = sp.acos((xx_to_Cart[2] / xxSph[0]))
        xxSph[2] = sp.atan2(xx_to_Cart[1], xx_to_Cart[0])
        Cart_to_xx[0] = (SINHWXYZ * sp.asinh(((Cartx * sp.sinh((1 / SINHWXYZ))) / AMPLXYZ)))
        Cart_to_xx[1] = (SINHWXYZ * sp.asinh(((Carty * sp.sinh((1 / SINHWXYZ))) / AMPLXYZ)))
        Cart_to_xx[2] = (SINHWXYZ * sp.asinh(((Cartz * sp.sinh((1 / SINHWXYZ))) / AMPLXYZ)))
        scalefactor_orthog[0] = sp.diff(xx_to_Cart[0], xx[0])
        scalefactor_orthog[1] = sp.diff(xx_to_Cart[1], xx[1])
        scalefactor_orthog[2] = sp.diff(xx_to_Cart[2], xx[2])
        f0_of_xx0 = sp.diff(xx_to_Cart[0], xx[0])
        f1_of_xx1 = sp.diff(xx_to_Cart[1], xx[1])
        f4_of_xx2 = sp.diff(xx_to_Cart[2], xx[2])
        scalefactor_orthog_funcform[0] = f0_of_xx0_funcform
        scalefactor_orthog_funcform[1] = f1_of_xx1_funcform
        scalefactor_orthog_funcform[2] = f4_of_xx2_funcform
        UnitVectors = [[sp.sympify(1), sp.sympify(0), sp.sympify(0)], [sp.sympify(0), sp.sympify(1), sp.sympify(0)], [sp.sympify(0), sp.sympify(0), sp.sympify(1)]]
    else:
        print((('CoordSystem == ' + CoordSystem) + ' is not supported.'))
        sys.exit(1)
    ref_metric__hatted_quantities(SymPySimplifyExpressions)