#!/usr/bin/python
# Quick/rough Python code to check (numerically) the derivations...
import numpy as np
for i in range(10):
x = np.random.uniform(-10,10)
z = np.random.uniform(-10,10)
d = np.random.uniform(10,1000)
phi = np.random.uniform(-np.pi, +np.pi)
theta = np.random.uniform(0, +np.pi)
print( 'x{:+05.1f} d{:+06.1f} z{:+05.1f} phi{:+06.1f} theta{:+06.1f}'.format(x, d, z, np.degrees(theta), np.degrees(phi)) )
#k = 1.0
v2y = d*np.cos(theta) - z*np.sin(theta)
dprime = np.sqrt( np.square(x) + np.square(d) + np.square(z) )
qx = x*np.cos(phi) - np.sin(phi)*v2y
qy = x*np.sin(phi) + np.cos(phi)*v2y - dprime
qz = d*np.sin(theta) + z*np.cos(theta)
q = np.sqrt( np.square(qx) + np.square(qy) + np.square(qz) )
print( ' pieces qx{:g} qy{:g} qz{:g} q{:g}'.format(qx, qy, qz, q) )
#print( '{}{:g}'.format( ' '*50, q) )
if False:
l1 = np.square(qx)
l2 = np.square(qy)
l3 = np.square(qz)
qalt2 = l1 + l2 + l3
if False:
l1 = np.square(x*np.cos(phi) - np.sin(phi)*(d*np.cos(theta) - z*np.sin(theta)))
l2 = np.square(x*np.sin(phi) + np.cos(phi)*(d*np.cos(theta) - z*np.sin(theta)) - dprime)
l3 = np.square(d*np.sin(theta) + z*np.cos(theta))
qalt2 = l1 + l2 + l3
if True:
l1 = np.square(x)*np.square(np.cos(phi)) - 2*x*np.cos(phi)*np.sin(phi)*v2y + np.square(np.sin(phi))*np.square(v2y)
l2 = np.square(x)*np.square(np.sin(phi)) + x*np.sin(phi)*np.cos(phi)*v2y - dprime*x*np.sin(phi)
l3 = x*np.sin(phi)*np.cos(phi)*v2y + np.square(np.cos(phi))*np.square(v2y) - dprime*np.cos(phi)*v2y
l4 = -1*dprime*x*np.sin(phi) - dprime*np.cos(phi)*v2y + np.square(dprime)
l5 = np.square(d)*np.square(np.sin(theta)) + 2*d*z*np.sin(theta)*np.cos(theta) + np.square(z)*np.square(np.cos(theta))
qalt2 = l1 + l2 + l3 + l4 + l5
if True:
qalt2 = 2*dprime*( dprime - x*np.sin(phi) - np.cos(phi)*v2y )
qalt = np.sqrt(qalt2)
print( ' qalt {:g}'.format(qalt) )
print( '{}{:g}'.format( ' '*50, q-qalt) )