Talk:Geometry:WAXS 3D

From GISAXS
Jump to: navigation, search
#!/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) )