funcion.py

from sympy import sqrt,sin,atan2,cos
mu = 398600441500000.000000
R = 6378136.3            
J2 = 1.75553*(10**10)*1000**5
J3 = -2.61913*(10**11)*1000**6
J4 = -0.16193312050719E-05/(R**4 * mu)
J5 = -0.22771610163688E-06/(R**5 * mu)
J6 = 0.53964849049834E-06/(R**6 * mu)
J7 = -0.35136844210318E-06/(R**7 * mu)
J8 = -0.20251871520885E-06/(R**8 * mu)

def Fg(x,y,z):
    r2 =  x**2 + y**2 + z**2
    Fx = x*mu*(J2*R**2*(3*z**2/(2*(x**2 + y**2 + z**2)) - 0.5)/(x**2 + y**2 + z**2) + J3*R**3*(5*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 3*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(3/2) + J4*R**4*(35*z**4/(8*(x**2 + y**2 + z**2)**2) - 15*z**2/(4*(x**2 + y**2 + z**2)) + 0.375)/(x**2 + y**2 + z**2)**2 + J5*R**5*(63*z**5/(8*(x**2 + y**2 + z**2)**(5/2)) - 35*z**3/(4*(x**2 + y**2 + z**2)**(3/2)) + 15*z/(8*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) + 945*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - R**5*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(x**2 + y**2 + z**2)**(5/2) + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 945*R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 105*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**4*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**2 + 105*R**4*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + R**4*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + 15*R**3*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - R**3*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(x**2 + y**2 + z**2)**(3/2) - 15*R**3*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(3/2) - 3*R**2*z*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(3/2) + R**2*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(3/2) - mu*(-3*J2*R**2*x*z**2/(x**2 + y**2 + z**2)**3 - 2*J2*R**2*x*(3*z**2/(2*(x**2 + y**2 + z**2)) - 0.5)/(x**2 + y**2 + z**2)**2 - 3*J3*R**3*x*(5*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 3*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) + J3*R**3*(-15*x*z**3/(2*(x**2 + y**2 + z**2)**(5/2)) + 3*x*z/(2*(x**2 + y**2 + z**2)**(3/2)))/(x**2 + y**2 + z**2)**(3/2) - 4*J4*R**4*x*(35*z**4/(8*(x**2 + y**2 + z**2)**2) - 15*z**2/(4*(x**2 + y**2 + z**2)) + 0.375)/(x**2 + y**2 + z**2)**3 + J4*R**4*(-35*x*z**4/(2*(x**2 + y**2 + z**2)**3) + 15*x*z**2/(2*(x**2 + y**2 + z**2)**2))/(x**2 + y**2 + z**2)**2 - 5*J5*R**5*x*(63*z**5/(8*(x**2 + y**2 + z**2)**(5/2)) - 35*z**3/(4*(x**2 + y**2 + z**2)**(3/2)) + 15*z/(8*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(7/2) + J5*R**5*(-315*x*z**5/(8*(x**2 + y**2 + z**2)**(7/2)) + 105*x*z**3/(4*(x**2 + y**2 + z**2)**(5/2)) - 15*x*z/(8*(x**2 + y**2 + z**2)**(3/2)))/(x**2 + y**2 + z**2)**(5/2) + 3780*R**5*x*z**3*(-z**2/(x**2 + y**2 + z**2) + 1)*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**5 - R**5*x*z**2*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(9/2)) - 3.0*R**5*x*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) + 945*R**5*x*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) - 4725.0*R**5*x*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) + 2*R**5*x*z**2*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) - 5670*R**5*x*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 5*R**5*x*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(x**2 + y**2 + z**2)**(7/2) - 5*R**5*x*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 5*R**5*x*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 4725*R**5*x*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 945*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(2.18467428144253e-57*y*sin(4*atan2(y, x))/(x**2 + y**2) + 3.68214853051585e-58*y*cos(4*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**3 - R**5*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-315*x*z**4/(2*(x**2 + y**2 + z**2)**3) + 105*x*z**2/(2*(x**2 + y**2 + z**2)**2))/(x**2 + y**2 + z**2)**(5/2) - R**5*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-1.27673579762669e-56*x**2/(x**2 + y**2)**(3/2) - 1.91721188645097e-56*x*y/(x**2 + y**2)**(3/2) + 1.27673579762669e-56/sqrt(x**2 + y**2))*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(x**2 + y**2 + z**2)**(5/2) + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(-5.01935870174171e-56*y*sin(2*atan2(y, x))/(x**2 + y**2) - 2.48740662616188e-56*y*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(-945*x*z**3/(2*(x**2 + y**2 + z**2)**(5/2)) + 105*x*z/(2*(x**2 + y**2 + z**2)**(3/2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.06429596208709e-56*y*sin(3*atan2(y, x))/(x**2 + y**2) - 5.06321650926682e-57*y*cos(3*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) - 945*R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(-5.1157928486588e-58*y*sin(5*atan2(y, x))/(x**2 + y**2) - 1.95872830662192e-57*y*cos(5*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) - 315.0*R**4*x*z**3*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) - R**4*x*z**2*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**4) - 105*R**4*x*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 420*R**4*x*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 2*R**4*x*z**2*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 525*R**4*x*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 4*R**4*x*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**3 - 420*R**4*x*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - 4*R**4*x*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - 105*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-2.69305659624137e-49*y*sin(3*atan2(y, x))/(x**2 + y**2) - 5.46258257352023e-50*y*cos(3*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) - R**4*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-105*x*z**3/(2*(x**2 + y**2 + z**2)**(5/2)) + 15*x*z/(2*(x**2 + y**2 + z**2)**(3/2)))/(x**2 + y**2 + z**2)**2 - R**4*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))*(-7.71205966695584e-49*x**2/(x**2 + y**2)**(3/2) - 6.81362183555786e-49*x*y/(x**2 + y**2)**(3/2) + 7.71205966695584e-49/sqrt(x**2 + y**2))/(x**2 + y**2 + z**2)**2 + 105*R**4*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(2.41485716138944e-50*y*sin(4*atan2(y, x))/(x**2 + y**2) + 3.95701656819374e-50*y*cos(4*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**2 + R**4*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.3773922698255e-49*y*sin(2*atan2(y, x))/(x**2 + y**2) + 4.49194768085293e-49*y*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**2 + 30*R**3*x*z**3*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 15*R**3*x*z**2*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(7/2) - R**3*x*z**2*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(7/2)) - 45.0*R**3*x*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) - 60*R**3*x*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 + 3*R**3*x*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(x**2 + y**2 + z**2)**(5/2) + 45*R**3*x*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) + 15*R**3*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(-5.97574231983965e-42*y*sin(2*atan2(y, x))/(x**2 + y**2) - 4.08809276302108e-42*y*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**2 - R**3*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)*(2.1202136313377e-41*x**2/(x**2 + y**2)**(3/2) + 2.5914026863665e-42*x*y/(x**2 + y**2)**(3/2) - 2.1202136313377e-41/sqrt(x**2 + y**2))/(x**2 + y**2 + z**2)**(3/2) - 15*R**3*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-2.9169059202925e-42*y*sin(3*atan2(y, x))/(x**2 + y**2) + 5.7202092780239e-42*y*cos(3*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(3/2) - 3*R**2*x*z**3*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(7/2)) + 6*R**2*x*z**2*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 + 9*R**2*x*z*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(5/2) - 2*R**2*x*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - 3*R**2*z*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-1.48871617866734e-38*x**2/(x**2 + y**2)**(3/2) + 9.51631253760479e-38*x*y/(x**2 + y**2)**(3/2) + 1.48871617866734e-38/sqrt(x**2 + y**2))/(x**2 + y**2 + z**2)**(3/2) + R**2*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(-1.94203581611747e-34*y*sin(2*atan2(y, x))/(x**2 + y**2) - 1.11483263842116e-34*y*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2))/sqrt(x**2 + y**2 + z**2)
    Fy = y*mu*(J2*R**2*(3*z**2/(2*(x**2 + y**2 + z**2)) - 0.5)/(x**2 + y**2 + z**2) + J3*R**3*(5*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 3*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(3/2) + J4*R**4*(35*z**4/(8*(x**2 + y**2 + z**2)**2) - 15*z**2/(4*(x**2 + y**2 + z**2)) + 0.375)/(x**2 + y**2 + z**2)**2 + J5*R**5*(63*z**5/(8*(x**2 + y**2 + z**2)**(5/2)) - 35*z**3/(4*(x**2 + y**2 + z**2)**(3/2)) + 15*z/(8*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) + 945*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - R**5*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(x**2 + y**2 + z**2)**(5/2) + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 945*R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 105*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**4*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**2 + 105*R**4*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + R**4*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + 15*R**3*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - R**3*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(x**2 + y**2 + z**2)**(3/2) - 15*R**3*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(3/2) - 3*R**2*z*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(3/2) + R**2*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(3/2) - mu*(-3*J2*R**2*y*z**2/(x**2 + y**2 + z**2)**3 - 2*J2*R**2*y*(3*z**2/(2*(x**2 + y**2 + z**2)) - 0.5)/(x**2 + y**2 + z**2)**2 - 3*J3*R**3*y*(5*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 3*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) + J3*R**3*(-15*y*z**3/(2*(x**2 + y**2 + z**2)**(5/2)) + 3*y*z/(2*(x**2 + y**2 + z**2)**(3/2)))/(x**2 + y**2 + z**2)**(3/2) - 4*J4*R**4*y*(35*z**4/(8*(x**2 + y**2 + z**2)**2) - 15*z**2/(4*(x**2 + y**2 + z**2)) + 0.375)/(x**2 + y**2 + z**2)**3 + J4*R**4*(-35*y*z**4/(2*(x**2 + y**2 + z**2)**3) + 15*y*z**2/(2*(x**2 + y**2 + z**2)**2))/(x**2 + y**2 + z**2)**2 - 5*J5*R**5*y*(63*z**5/(8*(x**2 + y**2 + z**2)**(5/2)) - 35*z**3/(4*(x**2 + y**2 + z**2)**(3/2)) + 15*z/(8*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(7/2) + J5*R**5*(-315*y*z**5/(8*(x**2 + y**2 + z**2)**(7/2)) + 105*y*z**3/(4*(x**2 + y**2 + z**2)**(5/2)) - 15*y*z/(8*(x**2 + y**2 + z**2)**(3/2)))/(x**2 + y**2 + z**2)**(5/2) + 3780*R**5*y*z**3*(-z**2/(x**2 + y**2 + z**2) + 1)*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**5 - R**5*y*z**2*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(9/2)) - 3.0*R**5*y*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) + 945*R**5*y*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) - 4725.0*R**5*y*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) + 2*R**5*y*z**2*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) - 5670*R**5*y*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 5*R**5*y*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(x**2 + y**2 + z**2)**(7/2) - 5*R**5*y*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 5*R**5*y*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 4725*R**5*y*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 945*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-2.18467428144253e-57*x*sin(4*atan2(y, x))/(x**2 + y**2) - 3.68214853051585e-58*x*cos(4*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**3 - R**5*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-315*y*z**4/(2*(x**2 + y**2 + z**2)**3) + 105*y*z**2/(2*(x**2 + y**2 + z**2)**2))/(x**2 + y**2 + z**2)**(5/2) - R**5*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)*(-1.27673579762669e-56*x*y/(x**2 + y**2)**(3/2) - 1.91721188645097e-56*y**2/(x**2 + y**2)**(3/2) + 1.91721188645097e-56/sqrt(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(5.01935870174171e-56*x*sin(2*atan2(y, x))/(x**2 + y**2) + 2.48740662616188e-56*x*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(-945*y*z**3/(2*(x**2 + y**2 + z**2)**(5/2)) + 105*y*z/(2*(x**2 + y**2 + z**2)**(3/2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(-1.06429596208709e-56*x*sin(3*atan2(y, x))/(x**2 + y**2) + 5.06321650926682e-57*x*cos(3*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) - 945*R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(5.1157928486588e-58*x*sin(5*atan2(y, x))/(x**2 + y**2) + 1.95872830662192e-57*x*cos(5*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) - 315.0*R**4*y*z**3*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(9/2) - R**4*y*z**2*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**4) - 105*R**4*y*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 420*R**4*y*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 2*R**4*y*z**2*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 525*R**4*y*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 4*R**4*y*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**3 - 420*R**4*y*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - 4*R**4*y*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - 105*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(2.69305659624137e-49*x*sin(3*atan2(y, x))/(x**2 + y**2) + 5.46258257352023e-50*x*cos(3*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(5/2) - R**4*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-105*y*z**3/(2*(x**2 + y**2 + z**2)**(5/2)) + 15*y*z/(2*(x**2 + y**2 + z**2)**(3/2)))/(x**2 + y**2 + z**2)**2 - R**4*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))*(-7.71205966695584e-49*x*y/(x**2 + y**2)**(3/2) - 6.81362183555786e-49*y**2/(x**2 + y**2)**(3/2) + 6.81362183555786e-49/sqrt(x**2 + y**2))/(x**2 + y**2 + z**2)**2 + 105*R**4*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-2.41485716138944e-50*x*sin(4*atan2(y, x))/(x**2 + y**2) - 3.95701656819374e-50*x*cos(4*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**2 + R**4*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(2.3773922698255e-49*x*sin(2*atan2(y, x))/(x**2 + y**2) - 4.49194768085293e-49*x*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**2 + 30*R**3*y*z**3*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 15*R**3*y*z**2*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(7/2) - R**3*y*z**2*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(7/2)) - 45.0*R**3*y*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) - 60*R**3*y*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 + 3*R**3*y*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(x**2 + y**2 + z**2)**(5/2) + 45*R**3*y*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) + 15*R**3*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(5.97574231983965e-42*x*sin(2*atan2(y, x))/(x**2 + y**2) + 4.08809276302108e-42*x*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**2 - R**3*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)*(2.1202136313377e-41*x*y/(x**2 + y**2)**(3/2) + 2.5914026863665e-42*y**2/(x**2 + y**2)**(3/2) - 2.5914026863665e-42/sqrt(x**2 + y**2))/(x**2 + y**2 + z**2)**(3/2) - 15*R**3*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(2.9169059202925e-42*x*sin(3*atan2(y, x))/(x**2 + y**2) - 5.7202092780239e-42*x*cos(3*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2)**(3/2) - 3*R**2*y*z**3*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(7/2)) + 6*R**2*y*z**2*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 + 9*R**2*y*z*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(5/2) - 2*R**2*y*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - 3*R**2*z*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-1.48871617866734e-38*x*y/(x**2 + y**2)**(3/2) + 9.51631253760479e-38*y**2/(x**2 + y**2)**(3/2) - 9.51631253760479e-38/sqrt(x**2 + y**2))/(x**2 + y**2 + z**2)**(3/2) + R**2*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(1.94203581611747e-34*x*sin(2*atan2(y, x))/(x**2 + y**2) + 1.11483263842116e-34*x*cos(2*atan2(y, x))/(x**2 + y**2))/(x**2 + y**2 + z**2))/sqrt(x**2 + y**2 + z**2)
    Fz = z*mu*(J2*R**2*(3*z**2/(2*(x**2 + y**2 + z**2)) - 0.5)/(x**2 + y**2 + z**2) + J3*R**3*(5*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 3*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(3/2) + J4*R**4*(35*z**4/(8*(x**2 + y**2 + z**2)**2) - 15*z**2/(4*(x**2 + y**2 + z**2)) + 0.375)/(x**2 + y**2 + z**2)**2 + J5*R**5*(63*z**5/(8*(x**2 + y**2 + z**2)**(5/2)) - 35*z**3/(4*(x**2 + y**2 + z**2)**(3/2)) + 15*z/(8*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) + 945*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - R**5*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(x**2 + y**2 + z**2)**(5/2) + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 945*R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 105*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - R**4*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**2 + 105*R**4*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + R**4*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + 15*R**3*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - R**3*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(x**2 + y**2 + z**2)**(3/2) - 15*R**3*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(3/2) - 3*R**2*z*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(3/2) + R**2*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(3/2) - mu*(-2*J2*R**2*z*(3*z**2/(2*(x**2 + y**2 + z**2)) - 0.5)/(x**2 + y**2 + z**2)**2 + J2*R**2*(-3*z**3/(x**2 + y**2 + z**2)**2 + 3*z/(x**2 + y**2 + z**2))/(x**2 + y**2 + z**2) - 3*J3*R**3*z*(5*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 3*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) + J3*R**3*(-15*z**4/(2*(x**2 + y**2 + z**2)**(5/2)) + 9*z**2/(x**2 + y**2 + z**2)**(3/2) - 3/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(3/2) - 4*J4*R**4*z*(35*z**4/(8*(x**2 + y**2 + z**2)**2) - 15*z**2/(4*(x**2 + y**2 + z**2)) + 0.375)/(x**2 + y**2 + z**2)**3 + J4*R**4*(-35*z**5/(2*(x**2 + y**2 + z**2)**3) + 25*z**3/(x**2 + y**2 + z**2)**2 - 15*z/(2*(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**2 - 5*J5*R**5*z*(63*z**5/(8*(x**2 + y**2 + z**2)**(5/2)) - 35*z**3/(4*(x**2 + y**2 + z**2)**(3/2)) + 15*z/(8*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(7/2) + J5*R**5*(-315*z**6/(8*(x**2 + y**2 + z**2)**(7/2)) + 525*z**4/(8*(x**2 + y**2 + z**2)**(5/2)) - 225*z**2/(8*(x**2 + y**2 + z**2)**(3/2)) + 15/(8*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) - 5670*R**5*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**4 + 5*R**5*z*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(x**2 + y**2 + z**2)**(7/2) + 945*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(4*z**3/(x**2 + y**2 + z**2)**2 - 4*z/(x**2 + y**2 + z**2))*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - 5*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 5*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 4725*R**5*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2.5*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) - R**5*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-315*z**5/(2*(x**2 + y**2 + z**2)**3) + 210*z**3/(x**2 + y**2 + z**2)**2 - 105*z/(2*(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) - R**5*(1.27673579762669e-56*x/sqrt(x**2 + y**2) + 1.91721188645097e-56*y/sqrt(x**2 + y**2))*(z**3/(x**2 + y**2 + z**2)**2 - z/(x**2 + y**2 + z**2))*(315*z**4/(8*(x**2 + y**2 + z**2)**2) - 105*z**2/(4*(x**2 + y**2 + z**2)) + 1.875)/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(5/2)) - R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(945*z**2/(2*(x**2 + y**2 + z**2)) - 52.5)*(3.0*z**3/(x**2 + y**2 + z**2)**2 - 3.0*z/(x**2 + y**2 + z**2))*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) + 945*R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.20537132628962e-59*sin(4*atan2(y, x)) + 5.46168570360633e-58*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 + R**5*(-z**2/(x**2 + y**2 + z**2) + 1)*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))*(-945*z**4/(2*(x**2 + y**2 + z**2)**(5/2)) + 525*z**2/(x**2 + y**2 + z**2)**(3/2) - 105/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**(5/2) - R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-945*z**3/(x**2 + y**2 + z**2)**2 + 945*z/(x**2 + y**2 + z**2))*(1.68773883642227e-57*sin(3*atan2(y, x)) + 3.54765320695698e-57*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 945*R**5*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(5.0*z**3/(x**2 + y**2 + z**2)**2 - 5.0*z/(x**2 + y**2 + z**2))*(3.91745661324383e-58*sin(5*atan2(y, x)) - 1.02315856973176e-58*cos(5*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) + R**5*(2*z**3/(x**2 + y**2 + z**2)**2 - 2*z/(x**2 + y**2 + z**2))*(315*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 105*z/(2*sqrt(x**2 + y**2 + z**2)))*(1.24370331308094e-56*sin(2*atan2(y, x)) - 2.50967935087085e-56*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) + 525*R**4*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(7/2) + 4*R**4*z*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**3 - 105*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(3.0*z**3/(x**2 + y**2 + z**2)**2 - 3.0*z/(x**2 + y**2 + z**2))*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) - 420*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)**2*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - 4*R**4*z*(-z**2/(x**2 + y**2 + z**2) + 1)*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 - R**4*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-105*z**4/(2*(x**2 + y**2 + z**2)**(5/2)) + 60*z**2/(x**2 + y**2 + z**2)**(3/2) - 15/(2*sqrt(x**2 + y**2 + z**2)))/(x**2 + y**2 + z**2)**2 - R**4*(7.71205966695584e-49*x/sqrt(x**2 + y**2) + 6.81362183555786e-49*y/sqrt(x**2 + y**2))*(z**3/(x**2 + y**2 + z**2)**2 - z/(x**2 + y**2 + z**2))*(35*z**3/(2*(x**2 + y**2 + z**2)**(3/2)) - 15*z/(2*sqrt(x**2 + y**2 + z**2)))/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**2) + R**4*(-z**2/(x**2 + y**2 + z**2) + 1)*(-105*z**3/(x**2 + y**2 + z**2)**2 + 105*z/(x**2 + y**2 + z**2))*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + 105*R**4*(-z**2/(x**2 + y**2 + z**2) + 1)*(4*z**3/(x**2 + y**2 + z**2)**2 - 4*z/(x**2 + y**2 + z**2))*(-9.89254142048434e-51*sin(4*atan2(y, x)) + 6.03714290347361e-51*cos(4*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - 105*R**4*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(1.82086085784008e-50*sin(3*atan2(y, x)) - 8.97685532080458e-50*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) + R**4*(105*z**2/(2*(x**2 + y**2 + z**2)) - 7.5)*(2*z**3/(x**2 + y**2 + z**2)**2 - 2*z/(x**2 + y**2 + z**2))*(-2.24597384042647e-49*sin(2*atan2(y, x)) - 1.18869613491275e-49*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - 60*R**3*z**2*(-z**2/(x**2 + y**2 + z**2) + 1)*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**3 + 3*R**3*z*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)/(x**2 + y**2 + z**2)**(5/2) + 45*R**3*z*(-z**2/(x**2 + y**2 + z**2) + 1)**1.5*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(5/2) + 15*R**3*z*(2*z**3/(x**2 + y**2 + z**2)**2 - 2*z/(x**2 + y**2 + z**2))*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - R**3*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(-15*z**3/(x**2 + y**2 + z**2)**2 + 15*z/(x**2 + y**2 + z**2))/(x**2 + y**2 + z**2)**(3/2) - R**3*(-2.1202136313377e-41*x/sqrt(x**2 + y**2) - 2.5914026863665e-42*y/sqrt(x**2 + y**2))*(15*z**2/(2*(x**2 + y**2 + z**2)) - 1.5)*(z**3/(x**2 + y**2 + z**2)**2 - z/(x**2 + y**2 + z**2))/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(3/2)) - 15*R**3*(-z**2/(x**2 + y**2 + z**2) + 1)**0.5*(3.0*z**3/(x**2 + y**2 + z**2)**2 - 3.0*z/(x**2 + y**2 + z**2))*(-1.90673642600797e-42*sin(3*atan2(y, x)) - 9.72301973430834e-43*cos(3*atan2(y, x)))/(x**2 + y**2 + z**2)**(3/2) + 15*R**3*(-z**2/(x**2 + y**2 + z**2) + 1)*(2.04404638151054e-42*sin(2*atan2(y, x)) - 2.98787115991983e-42*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 + 9*R**2*z**2*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(5/2) - 3*R**2*z*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*(z**3/(x**2 + y**2 + z**2)**2 - z/(x**2 + y**2 + z**2))/(sqrt(-z**2/(x**2 + y**2 + z**2) + 1)*(x**2 + y**2 + z**2)**(3/2)) - 2*R**2*z*(-3*z**2/(x**2 + y**2 + z**2) + 3)*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2)**2 - 3*R**2*(1.48871617866734e-38*x/sqrt(x**2 + y**2) - 9.51631253760479e-38*y/sqrt(x**2 + y**2))*sqrt(-z**2/(x**2 + y**2 + z**2) + 1)/(x**2 + y**2 + z**2)**(3/2) + R**2*(6*z**3/(x**2 + y**2 + z**2)**2 - 6*z/(x**2 + y**2 + z**2))*(5.57416319210581e-35*sin(2*atan2(y, x)) - 9.71017908058733e-35*cos(2*atan2(y, x)))/(x**2 + y**2 + z**2))/sqrt(x**2 + y**2 + z**2)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
    return Fx,Fy,Fz