### # Copyright 2008-2011 Diamond Light Source Ltd. # This file is part of Diffcalc. # # Diffcalc is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # Diffcalc is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with Diffcalc. If not, see . ### from math import pi, sin, cos, acos, asin, atan2, sqrt from itertools import product try: from numpy import matrix from numpy.linalg import norm except ImportError: from numjy import matrix from numjy.linalg import norm from diffcalc.log import logging from diffcalc.hkl.calcbase import HklCalculatorBase from diffcalc.hkl.you.geometry import create_you_matrices, calcMU, calcPHI, \ calcCHI, calcETA from diffcalc.hkl.you.geometry import YouPosition from diffcalc.util import DiffcalcException, bound, angle_between_vectors,\ y_rotation from diffcalc.util import cross3, z_rotation, x_rotation from diffcalc.ub.calc import PaperSpecificUbCalcStrategy from diffcalc.hkl.you.constraints import NUNAME logger = logging.getLogger("diffcalc.hkl.you.calc") I = matrix('1 0 0; 0 1 0; 0 0 1') SMALL = 1e-6 TORAD = pi / 180 TODEG = 180 / pi PRINT_DEGENERATE = False def is_small(x): return abs(x) < SMALL def sign(x): if is_small(x): return 0 if x > 0: return 1 # x < 0 return -1 def normalised(vector): return vector * (1 / norm(vector)) def cut_at_minus_pi(value): if value < (-pi - SMALL): return value + 2 * pi if value >= pi + SMALL: return value - 2 * pi return value def _calc_N(Q, n): """Return N as described by Equation 31""" Q = normalised(Q) n = normalised(n) if is_small(angle_between_vectors(Q, n)): # Replace the reference vector with an alternative vector from Eq.(78) idx_min, _ = min(enumerate([abs(Q[0, 0]), abs(Q[1, 0]), abs(Q[2, 0])]), key=lambda v: v[1]) idx_1, idx_2 = [idx for idx in range(3) if idx != idx_min] qval = sqrt(Q[idx_1, 0] * Q[idx_1, 0] + Q[idx_2, 0] * Q[idx_2, 0]) n[idx_min, 0] = qval n[idx_1, 0] = -Q[idx_min, 0] * Q[idx_1, 0] / qval n[idx_2, 0] = -Q[idx_min, 0] * Q[idx_2, 0] / qval if is_small(norm(n)): n[idx_min, 0] = 0 n[idx_1, 0] = Q[idx_2, 0] / qval n[idx_2, 0] = -Q[idx_1, 0] / qval Qxn = cross3(Q, n) QxnxQ = cross3(Qxn, Q) QxnxQ = normalised(QxnxQ) Qxn = normalised(Qxn) return matrix([[Q[0, 0], QxnxQ[0, 0], Qxn[0, 0]], [Q[1, 0], QxnxQ[1, 0], Qxn[1, 0]], [Q[2, 0], QxnxQ[2, 0], Qxn[2, 0]]]) def _calc_angle_between_naz_and_qaz(theta, alpha, tau): # Equation 30: top = cos(tau) - sin(alpha) * sin(theta) bottom = cos(alpha) * cos(theta) if is_small(bottom): if is_small(cos(alpha)): raise ValueError('cos(alpha) is too small') if is_small(cos(theta)): raise ValueError('cos(theta) is too small') if is_small(sin(tau)): return 0. return acos(bound(top / bottom)) def youAnglesToHkl(pos, wavelength, UBmatrix): """Calculate miller indices from position in radians. """ [MU, DELTA, NU, ETA, CHI, PHI] = create_you_matrices(*pos.totuple()) q_lab = (NU * DELTA - I) * matrix([[0], [2 * pi / wavelength], [0]]) # 12 hkl = UBmatrix.I * PHI.I * CHI.I * ETA.I * MU.I * q_lab return hkl[0, 0], hkl[1, 0], hkl[2, 0] def _tidy_degenerate_solutions(pos, constraints): original = pos.inDegrees() detector_like_constraint = constraints.detector or constraints.naz nu_constrained_to_0 = is_small(pos.nu) and detector_like_constraint mu_constrained_to_0 = is_small(pos.mu) and 'mu' in constraints.sample delta_constrained_to_0 = is_small(pos.delta) and detector_like_constraint eta_constrained_to_0 = is_small(pos.eta) and 'eta' in constraints.sample phi_not_constrained = not 'phi' in constraints.sample if nu_constrained_to_0 and mu_constrained_to_0 and phi_not_constrained: # constrained to vertical 4-circle like mode if is_small(pos.chi): # phi || eta desired_eta = pos.delta / 2. eta_diff = desired_eta - pos.eta pos.eta = desired_eta pos.phi -= eta_diff if PRINT_DEGENERATE: print ('DEGENERATE: with chi=0, phi and eta are colinear:' 'choosing eta = delta/2 by adding % 7.3f to eta and ' 'removing it from phi. (mu=%s=0 only)' % (eta_diff * TODEG, NUNAME)) print ' original:', original elif delta_constrained_to_0 and eta_constrained_to_0 and phi_not_constrained: # constrained to horizontal 4-circle like mode if is_small(pos.chi - pi / 2): # phi || mu desired_mu = pos.nu / 2. mu_diff = desired_mu - pos.mu pos.mu = desired_mu pos.phi += mu_diff if PRINT_DEGENERATE: print ('DEGENERATE: with chi=90, phi and mu are colinear: choosing' ' mu = %s/2 by adding % 7.3f to mu and to phi. ' '(delta=eta=0 only)' % (NUNAME, mu_diff * TODEG)) print ' original:', original return pos def _theta_and_qaz_from_detector_angles(delta, nu): # Equation 19: cos_2theta = cos(delta) * cos(nu) theta = acos(cos_2theta) / 2. sgn = sign(sin(2. * theta)) qaz = atan2(sgn * sin(delta), sgn * cos(delta) * sin(nu)) return theta, qaz class YouUbCalcStrategy(PaperSpecificUbCalcStrategy): def calculate_q_phi(self, pos): [MU, DELTA, NU, ETA, CHI, PHI] = create_you_matrices(*pos.totuple()) # Equation 12: Compute the momentum transfer vector in the lab frame y = matrix('0; 1; 0') q_lab = (NU * DELTA - I) * y # Transform this into the phi frame. return PHI.I * CHI.I * ETA.I * MU.I * q_lab UNREACHABLE_MSG = ( 'The current combination of constraints with %s = %.4f\n' 'prohibits a solution for the specified reflection.') class YouHklCalculator(HklCalculatorBase): def __init__(self, ubcalc, geometry, hardware, constraints, raiseExceptionsIfAnglesDoNotMapBackToHkl=True): HklCalculatorBase.__init__(self, ubcalc, geometry, hardware, raiseExceptionsIfAnglesDoNotMapBackToHkl) self._hardware = hardware # for checking limits only self.constraints = constraints self.parameter_manager = constraints # TODO: remove need for this attr def __str__(self): return self.constraints.__str__() def _get_n_phi(self): return self._ubcalc.n_phi def _get_ubmatrix(self): return self._getUBMatrix() # for consistency def repr_mode(self): return repr(self.constraints.all) def _anglesToHkl(self, pos, wavelength): """Calculate miller indices from position in radians. """ return youAnglesToHkl(pos, wavelength, self._get_ubmatrix()) def _anglesToVirtualAngles(self, pos, _wavelength): """Calculate pseudo-angles in radians from position in radians. Return theta, qaz, alpha, naz, tau, psi and beta in a dictionary. """ # depends on surface normal n_lab. mu, delta, nu, eta, chi, phi = pos.totuple() theta, qaz = _theta_and_qaz_from_detector_angles(delta, nu) # (19) [MU, _, _, ETA, CHI, PHI] = create_you_matrices(mu, delta, nu, eta, chi, phi) Z = MU * ETA * CHI * PHI n_lab = Z * self._get_n_phi() alpha = asin(bound((-n_lab[1, 0]))) naz = atan2(n_lab[0, 0], n_lab[2, 0]) # (20) cos_tau = cos(alpha) * cos(theta) * cos(naz - qaz) + \ sin(alpha) * sin(theta) tau = acos(bound(cos_tau)) # (23) # Compute Tau using the dot product directly (THIS ALSO WORKS) # q_lab = ( (NU * DELTA - I ) * matrix([[0],[1],[0]]) # norm = norm(q_lab) # q_lab = matrix([[1],[0],[0]]) if norm == 0 else q_lab * (1/norm) # tau_from_dot_product = acos(bound(dot3(q_lab, n_lab))) sin_beta = 2 * sin(theta) * cos(tau) - sin(alpha) beta = asin(bound(sin_beta)) # (24) psi = next(self._calc_psi(alpha, theta, tau, qaz, naz)) return {'theta': theta, 'qaz': qaz, 'alpha': alpha, 'naz': naz, 'tau': tau, 'psi': psi, 'beta': beta} def _choose_single_solution(self, pos_virtual_angles_pairs_in_degrees): if len(pos_virtual_angles_pairs_in_degrees) == 1: return pos_virtual_angles_pairs_in_degrees[0] absolute_distances = [] for pos_, _ in pos_virtual_angles_pairs_in_degrees: absolute_distances.append(sum([abs(pos_.totuple()[i]) for i in (0, 3, 4, 5)])) shortest_solution_index = absolute_distances.index( min(absolute_distances)) pos, virtual_angles = pos_virtual_angles_pairs_in_degrees[shortest_solution_index] if logger.isEnabledFor(logging.INFO): msg = ('Multiple sample solutions found (choosing solution with ' 'shortest distance to all-zeros position):\n') i = 0 for (pos_, _), distance in zip(pos_virtual_angles_pairs_in_degrees, absolute_distances): msg += '*' if i == shortest_solution_index else '.' msg += ('mu=% 7.3f, delta=% 7.3f, nu=% 7.3f, eta=% 7.3f, chi=% 7.3f, phi=% 7.3f' % pos_.totuple()) msg += ' (distance=% 4.3f)\n' % (distance * TODEG) i += 1 msg += ':\n' logger.info(msg) return pos, virtual_angles def hklToAngles(self, h, k, l, wavelength, return_all_solutions=False): """ Return verified Position and all virtual angles in degrees from h, k & l and wavelength in Angstroms. The calculated Position is verified by checking that it maps back using anglesToHkl() to the requested hkl value. Those virtual angles fixed or generated while calculating the position are verified by by checking that they map back using anglesToVirtualAngles to the virtual angles for the given position. Throws a DiffcalcException if either check fails and raiseExceptionsIfAnglesDoNotMapBackToHkl is True, otherwise displays a warning. """ pos_virtual_angles_pairs = self._hklToAngles(h, k, l, wavelength, return_all_solutions) # in rad assert pos_virtual_angles_pairs pos_virtual_angles_pairs_in_degrees = [] for pos, virtual_angles in pos_virtual_angles_pairs: # to degrees: pos.changeToDegrees() for key, val in virtual_angles.items(): if val is not None: virtual_angles[key] = val * TODEG self._verify_pos_map_to_hkl(h, k, l, wavelength, pos) pos_virtual_angles_pairs_in_degrees.append((pos, virtual_angles)) if return_all_solutions: return pos_virtual_angles_pairs_in_degrees else: pos, virtual_angles = self._choose_single_solution(pos_virtual_angles_pairs_in_degrees) return pos, virtual_angles def hkl_to_all_angles(self, h, k, l, wavelength): return self.hklToAngles(h, k, l, wavelength, True) def _hklToAngles(self, h, k, l, wavelength, return_all_solutions=False): """(pos, virtualAngles) = hklToAngles(h, k, l, wavelength) --- with Position object pos and the virtual angles returned in degrees. Some modes may not calculate all virtual angles. """ if not self.constraints.is_fully_constrained(): raise DiffcalcException( "Diffcalc is not fully constrained.\n" "Type 'help con' for instructions") if not self.constraints.is_current_mode_implemented(): raise DiffcalcException( "Sorry, the selected constraint combination is valid but " "is not implemented. Type 'help con' for implemented combinations") # constraints are dictionaries ref_constraint = self.constraints.reference if ref_constraint: ref_constraint_name, ref_constraint_value = ref_constraint.items()[0] det_constraint = self.constraints.detector naz_constraint = self.constraints.naz samp_constraints = self.constraints.sample assert not (det_constraint and naz_constraint), ( "Two 'detector' constraints given") h_phi = self._get_ubmatrix() * matrix([[h], [k], [l]]) theta = self._calc_theta(h_phi, wavelength) tau = angle_between_vectors(h_phi, self._get_n_phi()) if is_small(sin(tau)) and ref_constraint: if ref_constraint_name == 'psi': raise DiffcalcException("Azimuthal angle 'psi' is undefined as reference and scattering vectors parallel.\n" "Please constrain one of the sample angles or choose different reference vector orientation.") elif ref_constraint_name == 'a_eq_b': raise DiffcalcException("Reference constraint 'a_eq_b' is redundant as reference and scattering vectors are parallel.\n" "Please constrain one of the sample angles or choose different reference vector orientation.") ### Reference constraint column ### if ref_constraint: # An angle for the reference vector (n) is given (Section 5.2) alpha, _ = self._calc_remaining_reference_angles( ref_constraint_name, ref_constraint_value, theta, tau) solution_tuples = [] if det_constraint or naz_constraint: if len(samp_constraints) == 1: for qaz, naz, delta, nu in self._calc_det_angles_given_det_or_naz_constraint( det_constraint, naz_constraint, theta, tau, alpha): for mu, eta, chi, phi in self._calc_sample_angles_from_one_sample_constraint( samp_constraints, h_phi, theta, alpha, qaz, naz): solution_tuples.append((mu, delta, nu, eta, chi, phi)) elif len(samp_constraints) == 2: if det_constraint: det_constraint_name, det_constraint_val = det_constraint.items()[0] for delta, nu, qaz in self._calc_remaining_detector_angles(det_constraint_name, det_constraint_val, theta): for mu, eta, chi, phi in self._calc_sample_angles_given_two_sample_and_detector( samp_constraints, qaz, theta, h_phi, self._get_n_phi()): solution_tuples.append((mu, delta, nu, eta, chi, phi)) else: raise DiffcalcException( 'No code yet to handle this combination of detector and sample constraints!') elif len(samp_constraints) == 2: if ref_constraint_name == 'psi': psi_vals = [ref_constraint_value,] else: psi_vals = self._calc_psi(alpha, theta, tau) for psi in psi_vals: angles = list(self._calc_sample_given_two_sample_and_reference( samp_constraints, h_phi, theta, psi)) solution_tuples.extend(angles) elif len(samp_constraints) == 3: for angles in self._calc_angles_given_three_sample_constraints( h, k, l, wavelength, return_all_solutions, samp_constraints, h_phi, theta): solution_tuples.append(angles) if not solution_tuples: raise DiffcalcException('No solutions were found. ' 'Please consider using an alternative set of constraints.') tidy_solutions = [_tidy_degenerate_solutions(YouPosition(*pos, unit='RAD'), self.constraints).totuple() for pos in solution_tuples] merged_solution_tuples = set(self._filter_angle_limits(tidy_solutions, not return_all_solutions)) if not merged_solution_tuples: raise DiffcalcException('No solutions were found matching existing hardware limits. ' 'Please consider using an alternative set of constraints.') #def _find_duplicate_angles(el): # idx, tpl = el # for tmp_tpl in filtered_solutions[idx:]: # if False not in [abs(x-y) < SMALL for x,y in zip(tmp_tpl, tpl)]: # return False # return True #merged_solution_tuples = filter(_find_duplicate_angles, enumerate(filtered_solutions, 1)) position_pseudo_angles_pairs = self._create_position_pseudo_angles_pairs(wavelength, merged_solution_tuples) if not position_pseudo_angles_pairs: raise DiffcalcException('No solutions were found. Please check hardware limits and ' 'consider using an alternative pseudo-angle constraints.') return position_pseudo_angles_pairs def _create_position_pseudo_angles_pairs(self, wavelength, merged_solution_tuples): position_pseudo_angles_pairs = [] for pos in merged_solution_tuples: # Create position position = YouPosition(*pos, unit='RAD') #position = _tidy_degenerate_solutions(position, self.constraints) #if position.phi <= -pi + SMALL: # position.phi += 2 * pi # pseudo angles calculated along the way were for the initial solution # and may be invalid for the chosen solution TODO: anglesToHkl need no # longer check the pseudo_angles as they will be generated with the # same function and it will prove nothing pseudo_angles = self._anglesToVirtualAngles(position, wavelength) is_sol = True for constraint in [self.constraints.reference, self.constraints.detector, self.constraints.naz]: try: constraint_name, constraint_value = constraint.items()[0] if constraint_name == 'a_eq_b': if not is_small(pseudo_angles['alpha'] - pseudo_angles['beta']): is_sol = False break else: if not is_small(constraint_value - pseudo_angles[constraint_name]): is_sol = False break except: continue if is_sol: position_pseudo_angles_pairs.append((position, pseudo_angles)) return position_pseudo_angles_pairs def _calc_theta(self, h_phi, wavelength): """Calculate theta using Equation1 """ q_length = norm(h_phi) if is_small(q_length): raise DiffcalcException('Reflection is unreachable as |Q| is too small') wavevector = 2 * pi / wavelength try: theta = asin(bound(q_length / (2 * wavevector))) except AssertionError: raise DiffcalcException( 'Reflection is unreachable as |Q| is too long') if is_small(cos(theta)): raise DiffcalcException( 'Reflection is unreachable as theta angle is too close to 90 deg') return theta def _calc_psi(self, alpha, theta, tau, qaz=None, naz=None): """Calculate psi from Eq. (18), (25) and (28) """ sin_tau = sin(tau) cos_theta = cos(theta) if is_small(sin_tau): # The reference vector is parallel to the scattering vector yield float('nan') elif is_small(cos_theta): # Reflection is unreachable as theta angle is too close to 90 deg yield float('nan') elif is_small(sin(theta)): # Reflection is unreachable as |Q| is too small yield float('nan') else: cos_psi = ((cos(tau) * sin(theta) - sin(alpha)) / cos_theta) # (28) if qaz is None or naz is None : try: acos_psi = acos(bound(cos_psi / sin_tau)) if is_small(acos_psi): yield 0. else: for psi in [acos_psi, -acos_psi]: yield psi except AssertionError: print ('WARNING: Diffcalc could not calculate an azimuth (psi)') yield float('nan') else: sin_psi = cos(alpha) * sin(qaz - naz) sgn = sign(sin_tau) eps = sin_psi**2 + cos_psi**2 sigma_ = eps/sin_tau**2 - 1 if not is_small(sigma_): print ('WARNING: Diffcalc could not calculate a unique azimuth ' '(psi) because of loss of accuracy in numerical calculation') yield float('nan') else: psi = atan2(sgn * sin_psi, sgn * cos_psi) yield psi def _calc_remaining_reference_angles(self, name, value, theta, tau): """Return alpha and beta given one of a_eq_b, alpha, beta or psi """ if name == 'psi': psi = value # Equation 26 for alpha sin_alpha = (cos(tau) * sin(theta) - cos(theta) * sin(tau) * cos(psi)) if abs(sin_alpha) > 1 + SMALL: raise DiffcalcException(UNREACHABLE_MSG % (name, value * TODEG)) alpha = asin(bound(sin_alpha)) # Equation 27 for beta sin_beta = cos(tau) * sin(theta) + cos(theta) * sin(tau) * cos(psi) if abs(sin_beta) > 1 + SMALL: raise DiffcalcException(UNREACHABLE_MSG % (name, value * TODEG)) beta = asin(bound(sin_beta)) elif name == 'a_eq_b': alpha = beta = asin(cos(tau) * sin(theta)) # (24) elif name == 'alpha': alpha = value # (24) sin_beta = 2 * sin(theta) * cos(tau) - sin(alpha) if abs(sin_beta) > 1 + SMALL: raise DiffcalcException(UNREACHABLE_MSG % (name, value * TODEG)) beta = asin(sin_beta) elif name == 'beta': beta = value sin_alpha = 2 * sin(theta) * cos(tau) - sin(beta) # (24) if abs(sin_alpha) > 1 + SMALL: raise DiffcalcException(UNREACHABLE_MSG % (name, value * TODEG)) alpha = asin(sin_alpha) return alpha, beta def _calc_det_angles_given_det_or_naz_constraint( self, det_constraint, naz_constraint, theta, tau, alpha): assert det_constraint or naz_constraint try: naz_qaz_angle = _calc_angle_between_naz_and_qaz(theta, alpha, tau) except AssertionError: return if det_constraint: # One of the detector angles is given (Section 5.1) det_constraint_name, det_constraint = det_constraint.items()[0] for delta, nu, qaz in self._calc_remaining_detector_angles( det_constraint_name, det_constraint, theta): if is_small(naz_qaz_angle): naz_angles = [qaz,] else: naz_angles = [qaz - naz_qaz_angle, qaz + naz_qaz_angle] for naz in naz_angles: yield qaz, naz, delta, nu elif naz_constraint: # The 'detector' angle naz is given: det_constraint_name, det_constraint = naz_constraint.items()[0] naz_name, naz = det_constraint_name, det_constraint assert naz_name == 'naz' if is_small(naz_qaz_angle): qaz_angles = [naz,] else: qaz_angles = [naz - naz_qaz_angle, naz + naz_qaz_angle] for qaz in qaz_angles: for delta, nu, _ in self._calc_remaining_detector_angles( 'qaz', qaz, theta): yield qaz, naz, delta, nu def _calc_remaining_detector_angles(self, constraint_name, constraint_value, theta): """Return delta, nu and qaz given one detector angle """ # (section 5.1) # Find qaz using various derivations of 17 and 18 sin_2theta = sin(2 * theta) cos_2theta = cos(2 * theta) if is_small(sin_2theta): raise DiffcalcException( 'No meaningful scattering vector (Q) can be found when ' 'theta is so small (%.4f).' % theta * TODEG) if constraint_name == 'delta': delta = constraint_value try: asin_qaz = asin(bound(sin(delta) / sin_2theta)) # (17 & 18) except AssertionError: return cos_delta = cos(delta) if is_small(cos_delta): #raise DiffcalcException( # 'The %s and %s circles are redundant when delta is constrained to %.0f degrees.' # 'Please change delta constraint or use 4-circle mode.' % (NUNAME, 'mu', delta * TODEG)) print (('DEGENERATE: with delta=90, %s is degenerate: choosing ' '%s = 0 (allowed because %s is unconstrained)') % (NUNAME, NUNAME, NUNAME)) acos_nu = 1. else: try: acos_nu = acos(bound(cos_2theta / cos_delta)) except AssertionError: return if is_small(cos(asin_qaz)): qaz_angles = [sign(asin_qaz) * pi / 2.,] else: qaz_angles = [asin_qaz, pi - asin_qaz] if is_small(acos_nu): nu_angles = [0.,] else: nu_angles = [acos_nu, -acos_nu] for qaz, nu in product(qaz_angles, nu_angles): sgn_ref = sign(sin_2theta) * sign(cos(qaz)) sgn_ratio = sign(sin(nu)) * sign(cos_delta) if sgn_ref == sgn_ratio: yield delta, nu, qaz elif constraint_name == NUNAME: nu = constraint_value cos_nu = cos(nu) if is_small(cos_nu): raise DiffcalcException( 'The %s circle constraint to %.0f degrees is redundant.' 'Please change this constraint or use 4-circle mode.' % (NUNAME, nu * TODEG)) cos_delta = cos_2theta / cos(nu) cos_qaz = cos_delta * sin(nu) / sin_2theta try: acos_delta = acos(bound(cos_delta)) acos_qaz = acos(bound(cos_qaz)) except AssertionError: return if is_small(acos_qaz): qaz_angles = [0.,] else: qaz_angles = [acos_qaz, -acos_qaz] if is_small(acos_delta): delta_angles = [0.,] else: delta_angles = [acos_delta, -acos_delta] for qaz, delta in product(qaz_angles, delta_angles): sgn_ref = sign(sin(delta)) sgn_ratio = sign(sin(qaz)) * sign(sin_2theta) if sgn_ref == sgn_ratio: yield delta, nu, qaz elif constraint_name == 'qaz': qaz = constraint_value asin_delta = asin(sin(qaz) * sin_2theta) if is_small(cos(asin_delta)): delta_angles = [sign(asin_delta) * pi / 2.,] else: delta_angles = [asin_delta, pi - asin_delta] for delta in delta_angles: cos_delta = cos(delta) if is_small(cos_delta): print (('DEGENERATE: with delta=90, %s is degenerate: choosing ' '%s = 0 (allowed because %s is unconstrained)') % (NUNAME, NUNAME, NUNAME)) #raise DiffcalcException( # 'The %s circle is redundant when delta is at %.0f degrees.' # 'Please change detector constraint or use 4-circle mode.' % (NUNAME, delta * TODEG)) nu = 0. else: sgn_delta = sign(cos_delta) nu = atan2(sgn_delta * sin_2theta * cos(qaz), sgn_delta * cos_2theta) yield delta, nu, qaz else: raise DiffcalcException( constraint_name + ' is not an explicit detector angle ' '(naz cannot be handled here)') def _calc_sample_angles_from_one_sample_constraint( self, samp_constraints, h_phi, theta, alpha, qaz, naz): sample_constraint_name, sample_value = samp_constraints.items()[0] q_lab = matrix([[cos(theta) * sin(qaz)], [-sin(theta)], [cos(theta) * cos(qaz)]]) # (18) n_lab = matrix([[cos(alpha) * sin(naz)], [-sin(alpha)], [cos(alpha) * cos(naz)]]) # (20) mu_eta_chi_phi_tuples = list(self._calc_remaining_sample_angles( sample_constraint_name, sample_value, q_lab, n_lab, h_phi, self._get_n_phi())) return mu_eta_chi_phi_tuples def _calc_sample_given_two_sample_and_reference( self, samp_constraints, h_phi, theta, psi): for angles in self._calc_sample_angles_given_two_sample_and_reference( samp_constraints, psi, theta, h_phi, self._get_n_phi()): qaz, psi, mu, eta, chi, phi = angles values_in_deg = tuple(v * TODEG for v in angles) logger.info('Initial angles: xi=%.3f, psi=%.3f, mu=%.3f, ' 'eta=%.3f, chi=%.3f, phi=%.3f' % values_in_deg) # Try to find a solution for each possible transformed xi logger.info("") msg = "---Trying psi=%.3f, qaz=%.3f" % (psi * TODEG, qaz * TODEG) logger.info(msg) for delta, nu, _ in self._calc_remaining_detector_angles('qaz', qaz, theta): logger.info("delta=%.3f, %s=%.3f", delta * TODEG, NUNAME, nu * TODEG) #for mu, eta, chi, phi in self._generate_sample_solutions( # mu, eta, chi, phi, samp_constraints.keys(), delta, # nu, wavelength, (h, k, l), ref_constraint_name, # ref_constraint_value): yield mu, delta, nu, eta, chi, phi def _calc_remaining_sample_angles(self, constraint_name, constraint_value, q_lab, n_lab, q_phi, n_phi): """Return phi, chi, eta and mu, given one of these""" # (section 5.3) N_lab = _calc_N(q_lab, n_lab) N_phi = _calc_N(q_phi, n_phi) Z = N_lab * N_phi.T if constraint_name == 'mu': # (35) mu = constraint_value V = calcMU(mu).I * N_lab * N_phi.T try: acos_chi = acos(bound(V[2, 2])) except AssertionError: return if is_small(sin(acos_chi)): # chi ~= 0 or 180 and therefor phi || eta The solutions for phi # and eta here will be valid but will be chosen unpredictably. # Choose eta=0: # # tan(phi+eta)=v12/v11 from docs/extensions_to_yous_paper.wxm chi = acos_chi eta = 0. phi = atan2(-V[1, 0], V[1, 1]) logger.debug( 'Eta and phi cannot be chosen uniquely with chi so close ' 'to 0 or 180. Returning phi=%.3f and eta=%.3f', phi * TODEG, eta * TODEG) yield mu, eta, chi, phi else: for chi in [acos_chi, -acos_chi]: sgn = sign(sin(chi)) phi = atan2(-sgn * V[2, 1], -sgn * V[2, 0]) eta = atan2(-sgn * V[1, 2], sgn * V[0, 2]) yield mu, eta, chi, phi elif constraint_name == 'phi': # (37) phi = constraint_value V = N_lab * N_phi.I * calcPHI(phi).T try: asin_eta = asin(bound(V[0, 1])) except AssertionError: return if is_small(cos(asin_eta)): raise DiffcalcException('Chi and mu cannot be chosen uniquely ' 'with eta so close to +/-90.') for eta in [asin_eta, pi - asin_eta]: sgn = sign(cos(eta)) mu = atan2(sgn * V[2, 1], sgn * V[1, 1]) chi = atan2(sgn * V[0, 2], sgn * V[0, 0]) yield mu, eta, chi, phi elif constraint_name in ('eta', 'chi'): if constraint_name == 'eta': # (39) eta = constraint_value cos_eta = cos(eta) if is_small(cos_eta): #TODO: Not likely to happen in real world!? raise DiffcalcException( 'Chi and mu cannot be chosen uniquely with eta ' 'constrained so close to +-90.') try: asin_chi = asin(bound(Z[0, 2] / cos_eta)) except AssertionError: return all_eta = [eta,] all_chi = [asin_chi, pi - asin_chi] else: # constraint_name == 'chi' # (40) chi = constraint_value sin_chi = sin(chi) if is_small(sin_chi): raise DiffcalcException( 'Eta and phi cannot be chosen uniquely with chi ' 'constrained so close to 0. (Please contact developer ' 'if this case is useful for you)') try: acos_eta = acos(bound(Z[0, 2] / sin_chi)) except AssertionError: return all_eta = [acos_eta, -acos_eta] all_chi = [chi,] for chi, eta in product(all_chi, all_eta): top_for_mu = Z[2, 2] * sin(eta) * sin(chi) + Z[1, 2] * cos(chi) bot_for_mu = -Z[2, 2] * cos(chi) + Z[1, 2] * sin(eta) * sin(chi) if is_small(top_for_mu) and is_small(bot_for_mu): # chi == +-90, eta == 0/180 and therefore phi || mu cos(chi) == # 0 and sin(eta) == 0 Experience shows that even though e.g. # the z[2, 2] and z[1, 2] values used to calculate mu may be # basically 0 (1e-34) their ratio in every case tested so far # still remains valid and using them will result in a phi # solution that is continuous with neighbouring positions. # # We cannot test phi minus mu here unfortunately as the final # phi and mu solutions have not yet been chosen (they may be # +-x or 180+-x). Otherwise we could choose a sensible solution # here if the one found was incorrect. # tan(phi+eta)=v12/v11 from extensions_to_yous_paper.wxm phi_minus_mu = -atan2(Z[2, 0], Z[1, 1]) logger.debug( 'Mu and phi cannot be chosen uniquely with chi so close ' 'to +/-90 and eta so close 0 or 180.\n After the final ' 'solution has been chose phi-mu should equal: %.3f', phi_minus_mu * TODEG) mu = atan2(-top_for_mu, -bot_for_mu) # (41) top_for_phi = Z[0, 1] * cos(eta) * cos(chi) - Z[0, 0] * sin(eta) bot_for_phi = Z[0, 1] * sin(eta) + Z[0, 0] * cos(eta) * cos(chi) phi = atan2(top_for_phi, bot_for_phi) # (42) # if is_small(bot_for_phi) and is_small(top_for_phi): # raise DiffcalcException( # 'phi=%.3f cannot be known with confidence as top and ' # 'bottom are both close to zero. chi=%.3f, eta=%.3f' # % (mu * TODEG, chi * TODEG, eta * TODEG)) yield mu, eta, chi, phi else: raise DiffcalcException('Given angle must be one of phi, chi, eta or mu') def _calc_angles_given_three_sample_constraints( self, h, k, l, wavelength, return_all_solutions, samp_constraints, h_phi, theta): if not 'mu' in samp_constraints: eta_ = self.constraints.sample['eta'] chi_ = self.constraints.sample['chi'] phi_ = self.constraints.sample['phi'] try: two_mu_qaz_pairs = _mu_and_qaz_from_eta_chi_phi(eta_, chi_, phi_, theta, h_phi) except AssertionError: return else: raise DiffcalcException( 'No code yet to handle this combination of 3 sample constraints!') # TODO: Code duplicated above for mu_, qaz in two_mu_qaz_pairs: logger.debug("--- Trying mu_:%.f qaz_%.f", mu_ * TODEG, qaz * TODEG) for delta, nu, _ in self._calc_remaining_detector_angles('qaz', qaz, theta): logger.info("delta=%.3f, %s=%.3f", delta * TODEG, NUNAME, nu * TODEG) yield mu_, delta, nu, eta_, chi_, phi_ def _calc_sample_angles_given_two_sample_and_reference( self, samp_constraints, psi, theta, q_phi, n_phi): """Available combinations: chi, phi, reference mu, eta, reference, chi, eta, reference chi, mu, reference """ def __get_phi_and_qaz(chi, eta, mu): a = sin(chi) * cos(eta) b = sin(chi) * sin(eta) * sin(mu) - cos(chi) * cos(mu) #atan2_xi = atan2(V[2, 2] * a + V[2, 0] * b, # V[2, 0] * a - V[2, 2] * b) # (54) qaz = atan2(V[2, 0] * a - V[2, 2] * b, -V[2, 2] * a - V[2, 0] * b) # (54) a = sin(chi) * sin(mu) - cos(mu) * cos(chi) * sin(eta) b = cos(mu) * cos(eta) phi = atan2(V[1, 1] * a - V[0, 1] * b, V[0, 1] * a + V[1, 1] * b) # (55) # if is_small(mu+pi/2) and is_small(eta) and False: # phi_general = phi # # solved in extensions_to_yous_paper.wxm # phi = atan2(V[1, 1], V[0, 1]) # logger.info("phi = %.3f or %.3f (std)", # phi*TODEG, phi_general*TODEG ) return qaz, phi N_phi = _calc_N(q_phi, n_phi) THETA = z_rotation(-theta) PSI = x_rotation(psi) if 'chi' in samp_constraints and 'phi' in samp_constraints: chi = samp_constraints['chi'] phi = samp_constraints['phi'] CHI = calcCHI(chi) PHI = calcPHI(phi) V = CHI * PHI * N_phi * PSI.T * THETA.T # (46) #atan2_xi = atan2(-V[2, 0], V[2, 2]) #atan2_eta = atan2(-V[0, 1], V[1, 1]) #atan2_mu = atan2(-V[2, 1], sqrt(V[2, 2] ** 2 + V[2, 0] ** 2)) try: asin_mu = asin(bound(-V[2, 1])) except AssertionError: return for mu in [asin_mu, pi - asin_mu]: sgn_cosmu = sign(cos(mu)) #xi = atan2(-sgn_cosmu * V[2, 0], sgn_cosmu * V[2, 2]) qaz = atan2(sgn_cosmu * V[2, 2], sgn_cosmu * V[2, 0], ) eta = atan2(-sgn_cosmu * V[0, 1], sgn_cosmu * V[1, 1]) yield qaz, psi, mu, eta, chi, phi elif 'mu' in samp_constraints and 'eta' in samp_constraints: mu = samp_constraints['mu'] eta = samp_constraints['eta'] V = N_phi * PSI.T * THETA.T # (49) try: bot = bound(-V[2, 1] / sqrt(sin(eta) ** 2 * cos(mu) ** 2 + sin(mu) ** 2)) except AssertionError: return if is_small(cos(mu) * sin(eta)): eps = atan2(sin(eta) * cos(mu), sin(mu)) chi_vals = [eps + acos(bot), eps - acos(bot)] else: eps = atan2(sin(mu), sin(eta) * cos(mu)) chi_vals = [asin(bot) - eps, pi - asin(bot) - eps] # (52) ## Choose final chi solution here to obtain compatable xi and mu ## TODO: This temporary solution works only for one case used on i07 ## Return a list of possible solutions? #if is_small(eta) and is_small(mu + pi / 2): # for chi in _generate_transformed_values(chi_orig): # if pi / 2 <= chi < pi: # break #else: # chi = chi_orig for chi in chi_vals: qaz, phi = __get_phi_and_qaz(chi, eta, mu) yield qaz, psi, mu, eta, chi, phi elif 'chi' in samp_constraints and 'eta' in samp_constraints: chi = samp_constraints['chi'] eta = samp_constraints['eta'] V = N_phi * PSI.T * THETA.T # (49) try: bot = bound(-V[2, 1] / sqrt(sin(eta) ** 2 * sin(chi) ** 2 + cos(chi) ** 2)) except AssertionError: return if is_small(cos(chi)): eps = atan2(cos(chi), sin(chi) * sin(eta)) mu_vals = [eps + acos(bot), eps - acos(bot)] else: eps = atan2(sin(chi) * sin(eta), cos(chi)) mu_vals = [asin(bot) - eps, pi - asin(bot) - eps] # (52) for mu in mu_vals: qaz, phi = __get_phi_and_qaz(chi, eta, mu) yield qaz, psi, mu, eta, chi, phi elif 'chi' in samp_constraints and 'mu' in samp_constraints: chi = samp_constraints['chi'] mu = samp_constraints['mu'] V = N_phi * PSI.T * THETA.T # (49) try: asin_eta = asin(bound((-V[2, 1] - cos(chi) * sin(mu)) / (sin(chi) * cos(mu)))) except AssertionError: return for eta in [asin_eta, pi - asin_eta]: qaz, phi = __get_phi_and_qaz(chi, eta, mu) yield qaz, psi, mu, eta, chi, phi else: raise DiffcalcException( 'No code yet to handle this combination of 2 sample ' 'constraints and one reference!:' + str(samp_constraints)) def _calc_sample_angles_given_two_sample_and_detector( self, samp_constraints, qaz, theta, q_phi, n_phi): """Available combinations: chi, phi, detector mu, eta, detector mu, phi, detector """ N_phi = _calc_N(q_phi, n_phi) if 'mu' in samp_constraints and 'eta' in samp_constraints: mu = samp_constraints['mu'] eta = samp_constraints['eta'] F = y_rotation(qaz - pi/2.) THETA = z_rotation(-theta) V = calcETA(eta).T * calcMU(mu).T * F * THETA # (56) try: bot = bound(-V[1, 0] / sqrt(N_phi[0, 0]**2 + N_phi[1, 0]**2)) eps = atan2(N_phi[1, 0], N_phi[0, 0]) phi_vals = [asin(bot) + eps, pi - asin(bot) + eps] # (59) except (AssertionError, ZeroDivisionError): # For the case of (00l) reflection, where N_phi[0,0] = N_phi[1,0] = 0 chi = atan2(V[0, 0] * N_phi[2, 0], V[2, 0] * N_phi[2, 0]) # (57) sgn_denom = sign(N_phi[1, 1] * N_phi[0, 2] - N_phi[1, 2] * N_phi[0, 1]) sin_phi = V[1, 1] * N_phi[1, 2] - V[1, 2] * N_phi[1, 1] cos_phi = V[1, 1] * N_phi[0, 2] - V[1, 2] * N_phi[0, 1] phi = atan2(sin_phi * sgn_denom, cos_phi * sgn_denom) yield mu, eta, chi, phi return for phi in phi_vals: a = N_phi[0, 0] * cos(phi) + N_phi[1, 0] * sin(phi) chi = atan2(N_phi[2, 0] * V[0, 0] - a * V[2, 0], N_phi[2, 0] * V[2, 0] + a * V[0, 0]) # (60) yield mu, eta, chi, phi elif 'chi' in samp_constraints and 'phi' in samp_constraints: chi = samp_constraints['chi'] phi = samp_constraints['phi'] CHI = calcCHI(chi) PHI = calcPHI(phi) V = CHI * PHI * N_phi # (62) try: bot = bound(V[2, 0] / sqrt(cos(qaz) ** 2 * cos(theta) ** 2 + sin(theta) ** 2)) except AssertionError: return eps = atan2(-cos(qaz) * cos(theta), sin(theta)) for mu in [asin(bot) + eps, pi - asin(bot) + eps]: a = cos(theta) * sin(qaz) b = -cos(theta) * sin(mu) * cos(qaz) + cos(mu) * sin(theta) eta = atan2(V[1, 0] * a + V[0, 0] * b, V[0, 0] * a - V[1, 0]* b) #a = -cos(mu) * cos(qaz) * sin(theta) + sin(mu) * cos(theta) #b = cos(mu) * sin(qaz) #psi = atan2(-V[2, 2] * a - V[2, 1] * b, V[2, 1] * a - V[2, 2] * b) yield mu, eta, chi, phi elif 'mu' in samp_constraints and 'phi' in samp_constraints: mu = samp_constraints['mu'] phi = samp_constraints['phi'] F = y_rotation(qaz - pi/2.) THETA = z_rotation(-theta) V = calcMU(mu).T * F * THETA E = calcPHI(phi) * N_phi try: bot = bound(-V[2, 0] / sqrt(E[0, 0]**2 + E[2, 0]**2)) except AssertionError: return eps = atan2(E[2, 0], E[0, 0]) for chi in [asin(bot) + eps, pi - asin(bot) + eps]: a = E[0, 0] * cos(chi) + E[2, 0] * sin(chi) eta = atan2(V[0, 0] * E[1, 0] - V[1, 0] * a, V[0, 0] * a + V[1, 0] * E[1, 0]) yield mu, eta, chi, phi else: raise DiffcalcException( 'No code yet to handle this combination of 2 sample ' 'constraints and one detector!:' + str(samp_constraints)) def _filter_angle_limits(self, possible_solutions, filter_out_of_limits=True): res = [] angle_names = self._hardware.get_axes_names() for possible_solution in possible_solutions: hw_sol = [] hw_possible_solution = self._geometry.internal_position_to_physical_angles(YouPosition(*possible_solution, unit='RAD')) for name, value in zip(angle_names, hw_possible_solution): hw_sol.append(self._hardware.cut_angle(name, value)) if filter_out_of_limits: is_in_limits = all([self._hardware.is_axis_value_within_limits(name, value) for name, value in zip(angle_names, hw_sol)]) else: is_in_limits = True if is_in_limits: sol = self._geometry.physical_angles_to_internal_position(tuple(hw_sol)) sol.changeToRadians() res.append(sol.totuple()) return res def _mu_and_qaz_from_eta_chi_phi(eta, chi, phi, theta, h_phi): h_phi_norm = normalised(h_phi) # (68,69) h1, h2, h3 = h_phi_norm[0, 0], h_phi_norm[1, 0], h_phi_norm[2, 0] a = sin(chi) * h2 * sin(phi) + sin(chi) * h1 * cos(phi) - cos(chi) * h3 b = (- cos(chi) * sin(eta) * h2 * sin(phi) - cos(eta) * h1 * sin(phi) + cos(eta) * h2 * cos(phi) - cos(chi) * sin(eta) * h1 * cos(phi) - sin(chi) * sin(eta) * h3) c = -sin(theta) sin_bit = bound(c / sqrt(a * a + b * b)) mu1 = asin(sin_bit) - atan2(b, a) mu2 = pi - asin(sin_bit) - atan2(b, a) mu1 = cut_at_minus_pi(mu1) mu2 = cut_at_minus_pi(mu2) # TODO: This special case should be *removed* when the general case has shown # to encompass it. It exists as fallback for a particular i16 experiment in # May 2013 --RobW. # if eta == chi == 0: # logger.debug("Testing against simplified equations for eta == chi == 0") # a = - h3 # b = - h1 * sin(phi) + h2 * cos(phi) # sin_bit = bound(c / sqrt(a * a + b * b)) # mu_simplified = pi - asin(sin_bit) - atan2(b, a) # mu_simplified = cut_at_minus_pi(mu_simplified) # if not ne(mu_simplified, mu): # raise AssertionError("mu_simplified != mu , %f!=%f" % (mu_simplified, mu)) [MU, _, _, ETA, CHI, PHI] = create_you_matrices(mu1, None, None, eta, chi, phi) h_lab = MU * ETA * CHI * PHI * h_phi # (11) qaz1 = atan2(h_lab[0, 0] , h_lab[2, 0]) [MU, _, _, ETA, CHI, PHI] = create_you_matrices(mu2, None, None, eta, chi, phi) h_lab = MU * ETA * CHI * PHI * h_phi # (11) qaz2 = atan2(h_lab[0, 0] , h_lab[2, 0]) return (mu1, qaz1) , (mu2, qaz2)