cctbx/crystal_orientation.h

#ifndef CCTBX_CRYSTL_ORIENT_H
#define CCTBX_CRYSTL_ORIENT_H

#include <cctbx/uctbx.h>
#include <cctbx/sgtbx/change_of_basis_op.h>

namespace cctbx {

  typedef scitbx::vec3<double> oc_vec3;
  typedef scitbx::mat3<double> oc_mat3;
  typedef scitbx::af::shared<oc_vec3> oc_vec3_list;

  enum basis_type { direct=0, reciprocal=1 };


  template <typename FloatType>
  oc_vec3
  vector_rotate_thru(oc_vec3 const& vec, oc_vec3 const& unit,
                     FloatType const& angle){
    //
    //Still need to add assertion that unit is indeed a unit vector
    //
    return unit*(unit*vec)
         +(vec-(unit*(unit*vec)))*std::cos(angle)
         -(vec.cross(unit))*std::sin(angle);
  }


  class crystal_orientation
  {
    public:

      crystal_orientation(){}


      explicit
      crystal_orientation(oc_mat3 const&, bool const&);

      cctbx::uctbx::unit_cell
      unit_cell() const;

      cctbx::uctbx::unit_cell
      unit_cell_inverse() const;

      oc_mat3
      direct_matrix() const;

      oc_mat3
      reciprocal_matrix() const;


      crystal_orientation
      change_basis(cctbx::sgtbx::change_of_basis_op const&) const;


      crystal_orientation
      change_basis(oc_mat3 const&) const;


      template <typename FloatType>
      crystal_orientation
      rotate_thru(oc_vec3 const& unit_axis, FloatType const& angle){

        //
        //still need to add an assertion that unit is a unit vector
        //
        oc_vec3_list abc; //astar,bstar,cstar
        for (int i=0; i<3; ++i){
           abc.push_back( Astar_.get_column(i) );
        }

        oc_vec3_list newabc;
        for (int i=0; i<3; ++i){
           newabc.push_back( vector_rotate_thru(abc[i],unit_axis,angle) );
        }

        oc_mat3 new_recip_matrix;
        for (int i=0; i<3; ++i){
           new_recip_matrix.set_column( i, newabc[i] );
        }

        return crystal_orientation(new_recip_matrix, cctbx::reciprocal);
      }


      inline double
      direct_mean_square_difference(crystal_orientation const& other) const
      {
        return cctbx::uctbx::mean_square_difference(
                 direct_matrix(),other.direct_matrix())/3;
      }


      inline double
      difference_Z_score(crystal_orientation const& other) const
      {
        oc_mat3 diff = (direct_matrix() - other.direct_matrix());
        double Z = 0.0;
        for (int idx = 0; idx<3; ++idx){
          oc_vec3 this_basis_vector(direct_matrix().get_row(idx));
          oc_vec3 diff_vector(diff.get_row(idx));
          Z += diff_vector.length()/(0.01*this_basis_vector.length());
        }
        return Z;
      }


      oc_mat3
      best_similarity_transformation(crystal_orientation const& other,
       double const& fractional_length_tolerance,
       int unimodular_generator_range=1) const;

    protected:
      oc_mat3 Astar_;

  }; //class crystal_orientation
} //namespace cctbx

#endif // CCTBX_CRYSTL_ORIENT_H

---