coords.h

//C Copyright (C) 2000-2006 Kevin Cowtan and University of York
//L
//L  This library is free software and is distributed under the terms
//L  and conditions of version 2.1 of the GNU Lesser General Public
//L  Licence (LGPL) with the following additional clause:
//L
//L     `You may also combine or link a "work that uses the Library" to
//L     produce a work containing portions of the Library, and distribute
//L     that work under terms of your choice, provided that you give
//L     prominent notice with each copy of the work that the specified
//L     version of the Library is used in it, and that you include or
//L     provide public access to the complete corresponding
//L     machine-readable source code for the Library including whatever
//L     changes were used in the work. (i.e. If you make changes to the
//L     Library you must distribute those, but you do not need to
//L     distribute source or object code to those portions of the work
//L     not covered by this licence.)'
//L
//L  Note that this clause grants an additional right and does not impose
//L  any additional restriction, and so does not affect compatibility
//L  with the GNU General Public Licence (GPL). If you wish to negotiate
//L  other terms, please contact the maintainer.
//L
//L  You can redistribute it and/or modify the library under the terms of
//L  the GNU Lesser General Public License as published by the Free Software
//L  Foundation; either version 2.1 of the License, or (at your option) any
//L  later version.
//L
//L  This library is distributed in the hope that it will be useful, but
//L  WITHOUT ANY WARRANTY; without even the implied warranty of
//L  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//L  Lesser General Public License for more details.
//L
//L  You should have received a copy of the CCP4 licence and/or GNU
//L  Lesser General Public License along with this library; if not, write
//L  to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
//L  The GNU Lesser General Public can also be obtained by writing to the
//L  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
//L  MA 02111-1307 USA


#ifndef CLIPPER_COORDS
#define CLIPPER_COORDS


#include "cell.h"
#include "spacegroup.h"
#include "clipper_stats.h"


namespace clipper
{
  // forward definitions
  class Grid; class Grid_sampling; class Grid_range;
  class Coord_grid; class Coord_map;
  class Coord_reci_frac; class Coord_reci_orth;
  class Coord_frac; class Coord_orth;
  class U_aniso_frac; class U_aniso_orth;



  class Resolution
  {
  public:
    inline Resolution() : resol(0.0) {}          
    explicit Resolution( const ftype& resol_ );  
    void init( const ftype& resol_ );            
    const ftype& limit() const;                  
    ftype invresolsq_limit() const;  
    bool is_null() const;            
  private:
    ftype resol;
  };



  class HKL_class
  {
  public:
    inline HKL_class() { epsilon_ = 0; allowed_ = 255; }
    HKL_class( const Spacegroup& spgr, const HKL& hkl );
    inline ftype epsilon() const { return ftype(epsilon_); }
    inline ftype epsilonc() const
      { if ( centric() ) return 2.0*ftype(epsilon_);
        else             return ftype(epsilon_); }
    inline ftype allowed() const { return ftype(allowed_) * (Util::pi()/12.0); }
    inline bool centric() const { return allowed_ != 255; } 
    inline bool sys_abs() const { return epsilon_ == 0; }   
  private:
    unsigned char epsilon_, allowed_;
  };



  class RTop_orth : public RTop<>
  {
  public:
    inline RTop_orth() {}
    inline explicit RTop_orth( const RTop<>& o ) : RTop<>( o ) {}
    inline explicit RTop_orth( const Mat33<>& r ) : RTop<>( r ) {}
    inline RTop_orth( const Mat33<>& r, const Vec3<>& t ) : RTop<>( r, t ) {}
    RTop_orth( const std::vector<Coord_orth>& src, const std::vector<Coord_orth>& tgt );
    RTop_orth( const std::vector<Coord_orth>& src, const std::vector<Coord_orth>& tgt, const std::vector<ftype>& wgt );
    template<class T> RTop_orth( const T& src, const T& tgt );
    RTop_frac rtop_frac( const Cell& cell ) const;
    RTop_orth inverse() const;
    Coord_orth axis_coordinate_near( const Coord_orth& centre ) const;
    Coord_orth screw_translation() const;
    static RTop_orth identity();
    static RTop_orth null();
  };


  class HKL : public Vec3<int>
  {
  public:
    inline HKL() {}                
    inline explicit HKL( const Vec3<int>& v ) :
      Vec3<int>( v ) {}            
    inline HKL( const int& h, const int& k, const int& l ) :
      Vec3<int>( h, k, l ) {}      
    inline const int& h() const { return (*this)[0]; }  
    inline const int& k() const { return (*this)[1]; }  
    inline const int& l() const { return (*this)[2]; }  
    inline int& h() { return (*this)[0]; }  
    inline int& k() { return (*this)[1]; }  
    inline int& l() { return (*this)[2]; }  

    inline ftype invresolsq( const Cell& cell ) const;
    inline Coord_reci_frac coord_reci_frac() const;
    inline Coord_reci_orth coord_reci_orth( const Cell& cell ) const;
    inline HKL transform( const Symop& op ) const;
    inline HKL transform( const Isymop& op ) const;
    inline ftype sym_phase_shift( const Symop& op ) const;
    String format() const;  
    friend inline HKL operator -(const HKL& h1)
      { return HKL( -h1.h(), -h1.k(), -h1.l() ); }
    friend inline HKL operator +(const HKL& h1, const HKL& h2)
      { return HKL( h1.h()+h2.h(), h1.k()+h2.k(), h1.l()+h2.l() ); }
    friend inline HKL operator -(const HKL& h1, const HKL& h2)
      { return HKL( h1.h()-h2.h(), h1.k()-h2.k(), h1.l()-h2.l() ); }
    friend inline HKL operator *(const int& s, const HKL& h1)
      { return HKL( s*h1.h(), s*h1.k(), s*h1.l() ); }
    friend inline HKL operator *(const Isymop& op, const HKL& h1)
      { return HKL( h1 * op.rot() ); }
  };


  class Coord_reci_orth : public Vec3<>
  {
  public:
    inline Coord_reci_orth() {}           
    inline explicit Coord_reci_orth( const Vec3<>& v ) :
      Vec3<>( v ) {}          
    inline Coord_reci_orth( const ftype& xs, const ftype& ys, const ftype& zs ) : Vec3<>( xs, ys, zs ) {} 
    inline const ftype& xs() const { return (*this)[0]; }  
    inline const ftype& ys() const { return (*this)[1]; }  
    inline const ftype& zs() const { return (*this)[2]; }  

    inline ftype invresolsq() const;
    inline Coord_reci_frac coord_reci_frac( const Cell& cell ) const;
    inline Coord_reci_orth transform( const RTop_orth& op ) const
      { return Coord_reci_orth( (*this) * op.rot() ); }
    String format() const;  
  };


  class Coord_reci_frac : public Vec3<>
  {
  public:
    inline Coord_reci_frac() {}           
    inline explicit Coord_reci_frac( const Vec3<>& v ) :
      Vec3<>( v ) {}          
    inline Coord_reci_frac( const ftype& us, const ftype& vs, const ftype& ws ) : Vec3<>( us, vs, ws ) {} 

    inline Coord_reci_frac( const HKL& hkl ) :
      Vec3<>( ftype(hkl[0]), ftype(hkl[1]), ftype(hkl[2]) ) {}
    inline HKL hkl() const
      { return HKL( Util::intr(us()), Util::intr(vs()), Util::intr(ws()) ); }
    inline ftype invresolsq( const Cell& cell ) const;
    inline const ftype& us() const { return (*this)[0]; }  
    inline const ftype& vs() const { return (*this)[1]; }  
    inline const ftype& ws() const { return (*this)[2]; }  

    inline Coord_reci_orth coord_reci_orth( const Cell& cell ) const;
    inline Coord_reci_frac transform( const RTop_frac& op ) const
      { return Coord_reci_frac( (*this) * op.rot() ); }
    String format() const;  
  };


  class Coord_grid : public Vec3<int>
  {
  public:
    inline Coord_grid() {}  

    inline explicit Coord_grid( const Vec3<int> v ) : Vec3<int>( v ) {}
    inline Coord_grid( const int& u, const int& v, const int& w ) :
      Vec3<int>(u,v,w) {}
    inline Coord_grid( const Grid& g, const int& index )
      { deindex( g, index ); }
    inline const int& u() const { return (*this)[0]; }  
    inline const int& v() const { return (*this)[1]; }  
    inline const int& w() const { return (*this)[2]; }  
    inline int& u() { return (*this)[0]; }  
    inline int& v() { return (*this)[1]; }  
    inline int& w() { return (*this)[2]; }  

    inline Coord_map coord_map() const;
    inline Coord_frac coord_frac( const Grid_sampling& g ) const;
    inline Coord_grid transform( const Isymop& op ) const
      { return op * (*this); }

    inline Coord_grid unit( const Grid_sampling& g ) const;


    inline const Coord_grid& next( const Grid& g );

    inline const Coord_grid& next( const Grid_range& g );
    inline bool last( const Grid& g ) const;
    inline bool last( const Grid_range& g ) const;
    inline int index( const Grid& g ) const;
    inline void deindex( const Grid& g, const int& index );
    // grid indexing operator
    //inline int index( const Grid_range& g ) const;
    // grid deindexing operator
    //inline void deindex( const Grid_range& g, const int& index );

    String format() const;  
    friend inline Coord_grid operator -(const Coord_grid& r1)
      { return ( Coord_grid( -r1.u(), -r1.v(), -r1.w() ) ); }
    friend inline Coord_grid operator +(const Coord_grid& r1, const Coord_grid& r2) { return ( Coord_grid( r1.u()+r2.u(), r1.v()+r2.v(), r1.w()+r2.w() ) ); }
    friend inline Coord_grid operator -(const Coord_grid& r1, const Coord_grid& r2) { return ( Coord_grid( r1.u()-r2.u(), r1.v()-r2.v(), r1.w()-r2.w() ) ); }
    friend inline Coord_grid operator *(const int& s, const Coord_grid& r1)
      { return ( Coord_grid( s*r1.u(), s*r1.v(), s*r1.w() ) ); }
    friend inline int operator == (const Coord_grid& r1,  const Coord_grid& r2)
      { return (r1.u()==r2.u() && r1.v()==r2.v() && r1.w()==r2.w()); }
    friend inline int operator != (const Coord_grid& r1,  const Coord_grid& r2)
      { return (r1.u()!=r2.u() || r1.v()!=r2.v() || r1.w()!=r2.w()); }
    friend inline Coord_grid operator *(const Isymop& op, const Coord_grid& r1)
      { return Coord_grid( op.rot() * r1 + op.trn() ); }
  };


  class Coord_orth : public Vec3<>
  {
  public:
    inline Coord_orth() {}                
    inline explicit Coord_orth( const Vec3<>& v ) :
      Vec3<>( v ) {}          
    inline Coord_orth( const ftype& x, const ftype& y, const ftype& z ) :
      Vec3<>( x, y, z ) {}    

    Coord_orth( const Coord_orth& x1, const Coord_orth& x2, const Coord_orth& x3, const ftype& length, const ftype& angle, const ftype& torsion );
    inline const ftype& x() const { return (*this)[0]; }  
    inline const ftype& y() const { return (*this)[1]; }  
    inline const ftype& z() const { return (*this)[2]; }  

    inline ftype lengthsq() const;
    inline Coord_frac coord_frac( const Cell& cell ) const;
    inline Coord_orth transform( const RTop_orth& op ) const
      { return op*(*this); }
    String format() const;  

    static ftype length( const Coord_orth& x1, const Coord_orth& x2);
    static ftype angle( const Coord_orth& x1, const Coord_orth& x2, 
                        const Coord_orth& x3);
    static ftype torsion( const Coord_orth& x1, const Coord_orth& x2,
                          const Coord_orth& x3, const Coord_orth& x4);
    friend inline Coord_orth operator -(const Coord_orth& x1)
      { return Coord_orth( -x1.x(), -x1.y(), -x1.z() ); }
    friend inline Coord_orth operator +(const Coord_orth& x1, const Coord_orth& x2) { return Coord_orth( x1.x()+x2.x(), x1.y()+x2.y(), x1.z()+x2.z() ); }
    friend inline Coord_orth operator -(const Coord_orth& x1, const Coord_orth& x2) { return Coord_orth( x1.x()-x2.x(), x1.y()-x2.y(), x1.z()-x2.z() ); }
    friend inline Coord_orth operator *(const ftype& s, const Coord_orth& x1)
      { return Coord_orth( s*x1.x(), s*x1.y(), s*x1.z() ); }
    friend inline Coord_orth operator *(const RTop_orth& op, const Coord_orth& x1) { return Coord_orth( op.rot() * x1 + op.trn() ); }
  };


  class Coord_frac : public Vec3<>
  {
  public:
    inline Coord_frac() {}                
    inline explicit Coord_frac( const Vec3<>& v ) :
      Vec3<>( v ) {}          
    inline Coord_frac( const ftype& u, const ftype& v, const ftype& w ) :
      Vec3<>( u, v, w ) {}    
    inline const ftype& u() const { return (*this)[0]; }  
    inline const ftype& v() const { return (*this)[1]; }  
    inline const ftype& w() const { return (*this)[2]; }  

    inline ftype lengthsq( const Cell& cell ) const;
    inline Coord_orth coord_orth( const Cell& cell ) const;
    inline Coord_map coord_map( const Grid& g ) const;
    inline Coord_grid coord_grid( const Grid& g ) const;
    inline Coord_frac transform( const RTop_frac& op ) const
      { return op*(*this); }
    inline Coord_frac lattice_copy_zero() const
      { return Coord_frac(u()-rint(u()),v()-rint(v()),w()-rint(w())); }
    inline Coord_frac lattice_copy_unit() const
      { return Coord_frac(u()-floor(u()),v()-floor(v()),w()-floor(w())); }
    inline Coord_frac lattice_copy_near(const Coord_frac& n) const
      { return (*this-n).lattice_copy_zero()+n; }
    Coord_frac symmetry_copy_near(const Spacegroup& spgr, const Cell& cell, const Coord_frac& n) const;
    String format() const;  
    friend inline Coord_frac operator -(const Coord_frac& u1)
      { return Coord_frac( -u1.u(), -u1.v(), -u1.w() ); }
    friend inline Coord_frac operator +(const Coord_frac& u1, const Coord_frac& u2) { return Coord_frac( u1.u()+u2.u(), u1.v()+u2.v(), u1.w()+u2.w() ); }
    friend inline Coord_frac operator -(const Coord_frac& u1, const Coord_frac& u2) { return Coord_frac( u1.u()-u2.u(), u1.v()-u2.v(), u1.w()-u2.w() ); }
    friend inline Coord_frac operator *(const ftype& s, const Coord_frac& u1)
      { return Coord_frac( s*u1.u(), s*u1.v(), s*u1.w() ); }
    friend inline Coord_frac operator *(const RTop_frac& op, const Coord_frac& x1) { return Coord_frac( op.rot() * x1 + op.trn() ); }
  };


  class Coord_map : public Vec3<>
  {
  public:
    inline Coord_map() {}  

    inline explicit Coord_map( const Vec3<>& v ) :
      Vec3<>( v ) {}
    inline explicit Coord_map( const Coord_grid& c ) :
      Vec3<>( ftype(c[0]), ftype(c[1]), ftype(c[2]) ) {}
    inline Coord_map( const ftype& u, const ftype& v, const ftype& w ) :
      Vec3<>( u, v, w ) {}
    inline Coord_frac coord_frac( const Grid& g ) const;
    inline Coord_grid coord_grid() const { return Coord_grid( Util::intr((*this)[0]), Util::intr((*this)[1]), Util::intr((*this)[2]) ); }
    inline Coord_grid floor() const { return Coord_grid( Util::intf((*this)[0]), Util::intf((*this)[1]), Util::intf((*this)[2]) ); }
    inline Coord_grid ceil() const { return Coord_grid( Util::intc((*this)[0]), Util::intc((*this)[1]), Util::intc((*this)[2]) ); }
    inline const ftype& u() const { return (*this)[0]; }  
    inline const ftype& v() const { return (*this)[1]; }  
    inline const ftype& w() const { return (*this)[2]; }  
    String format() const;  
    friend inline Coord_map operator -(const Coord_map& u1)
      { return Coord_map( -u1.u(), -u1.v(), -u1.w() ); }
    friend inline Coord_map operator +(const Coord_map& u1, const Coord_map& u2)
      { return Coord_map( u1.u()+u2.u(), u1.v()+u2.v(), u1.w()+u2.w() ); }
    friend inline Coord_map operator -(const Coord_map& u1, const Coord_map& u2)
      { return Coord_map( u1.u()-u2.u(), u1.v()-u2.v(), u1.w()-u2.w() ); }
    friend inline Coord_map operator *(const ftype& s, const Coord_map& u1)
      { return Coord_map( s*u1.u(), s*u1.v(), s*u1.w() ); }
  };



  class U_aniso_orth : public Mat33sym<>
  {
  public:
    inline U_aniso_orth() {};
    inline explicit U_aniso_orth( const Mat33sym<>& m ) : Mat33sym<>(m) {}
    inline explicit U_aniso_orth( const ftype& u ) :
      Mat33sym<>( u, u, u, 0.0, 0.0, 0.0 ) {}
    U_aniso_orth( const ftype& u11, const ftype& u22, const ftype& u33,
                  const ftype& u12, const ftype& u13, const ftype& u23 ) :
      Mat33sym<>( u11, u22, u33, u12, u13, u23 ) {}
    ftype u_iso() const;
    U_aniso_frac u_aniso_frac( const Cell& cell ) const;
    U_aniso_orth transform( const RTop_orth& op ) const;
    friend U_aniso_orth operator +(const U_aniso_orth& u1, const U_aniso_orth& u2) { return U_aniso_orth( u1.mat00()+u2.mat00(), u1.mat11()+u2.mat11(), u1.mat22()+u2.mat22(), u1.mat01()+u2.mat01(), u1.mat02()+u2.mat02(), u1.mat12()+u2.mat12() ); }
    friend U_aniso_orth operator -(const U_aniso_orth& u) { return U_aniso_orth( -u.mat00(), -u.mat11(), -u.mat22(), -u.mat01(), -u.mat02(), -u.mat12() ); }
    friend U_aniso_orth operator *(const ftype& s, const U_aniso_orth& u) { return U_aniso_orth( s*u.mat00(), s*u.mat11(), s*u.mat22(), s*u.mat01(), s*u.mat02(), s*u.mat12() ); }
  };



  class U_aniso_frac : public Mat33sym<>
  {
  public:
    inline U_aniso_frac() {};
    inline explicit U_aniso_frac( const Mat33sym<>& m ) : Mat33sym<>(m) {}
    U_aniso_frac( const ftype& u11, const ftype& u22, const ftype& u33,
                  const ftype& u12, const ftype& u13, const ftype& u23 ) :
      Mat33sym<>( u11, u22, u33, u12, u13, u23 ) {}
    U_aniso_orth u_aniso_orth( const Cell& cell ) const;
    U_aniso_frac transform( const RTop_frac& op ) const;
    friend U_aniso_frac operator +(const U_aniso_frac& u1, const U_aniso_frac& u2) { return U_aniso_frac( u1.mat00()+u2.mat00(), u1.mat11()+u2.mat11(), u1.mat22()+u2.mat22(), u1.mat01()+u2.mat01(), u1.mat02()+u2.mat02(), u1.mat12()+u2.mat12() ); }
    friend U_aniso_frac operator -(const U_aniso_frac& u) { return U_aniso_frac( -u.mat00(), -u.mat11(), -u.mat22(), -u.mat01(), -u.mat02(), -u.mat12() ); }
    friend U_aniso_frac operator *(const ftype& s, const U_aniso_frac& u) { return U_aniso_frac( s*u.mat00(), s*u.mat11(), s*u.mat22(), s*u.mat01(), s*u.mat02(), s*u.mat12() ); }
  };



  class Grid : public Vec3<int>
  {
  public:
    inline Grid() {}                     
    inline Grid( const int& nu, const int& nv, const int& nw ) :
      Vec3<int>( nu, nv, nw ) {}  
    inline const int& nu() const { return (*this)[0]; } 
    inline const int& nv() const { return (*this)[1]; } 
    inline const int& nw() const { return (*this)[2]; } 

    inline int size() const { return nu()*nv()*nw(); }
    inline bool in_grid( Coord_grid g ) const { return (g.u() >= 0 && g.u() < nu() && g.v() >= 0 && g.v() < nv() && g.w() >= 0 && g.w() < nw()); }

    inline int index( const Coord_grid& c ) const { return c.index(*this); }
    inline Coord_grid deindex( const int& index ) const { return Coord_grid( *this, index ); }
    String format() const;  
    void debug() const;
  };



  class Grid_sampling : public Grid
  {
  public:
    inline Grid_sampling() : Grid(Grid(0,0,0)) {}
    inline Grid_sampling( const int& nu, const int& nv, const int& nw ) :
      Grid( nu, nv, nw ) {}
    Grid_sampling( const Spacegroup& spacegroup, const Cell& cell, const Resolution& resol, const ftype rate = 1.5 );
    void init( const Spacegroup& spacegroup, const Cell& cell, const Resolution& resol, const ftype rate = 1.5 );

    Mat33<> matrix_grid_frac() const;
    Mat33<> matrix_frac_grid() const;

    bool is_null() const;

    // inherited functions listed for documentation purposes
    //-- const int& nu() const;
    //-- const int& nv() const;
    //-- const int& nw() const;
    //-- int size() const;
    //-- int index( const Coord_grid& c ) const;
    //-- Coord_grid deindex( const int& index ) const;
    //-- const String format() const;
  };



  class HKL_sampling
  {
  public:
    HKL_sampling(); 

    HKL_sampling( const Cell& cell, const Resolution& resolution );
    HKL hkl_limit() const;
    Resolution resolution( const Cell& cell ) const;
    inline bool in_resolution( const HKL& h ) const
      { return ( m00*itype64(h.h()*h.h()) + m11*itype64(h.k()*h.k()) +
                 m22*itype64(h.l()*h.l()) + m01*itype64(h.h()*h.k()) +
                 m02*itype64(h.h()*h.l()) + m12*itype64(h.k()*h.l()) )
          <=   ( sqrt_limit_value*sqrt_limit_value ); }
    bool is_null() const;
    String format() const;  
    friend inline int operator == (const HKL_sampling& h1, const HKL_sampling& h2)
      { return ( h1.m00==h2.m00 && h1.m11==h2.m11 && h1.m22==h2.m22 &&
                 h1.m01==h2.m01 && h1.m02==h2.m02 && h1.m12==h2.m12 ); }
  private:
    static itype64 sqrt_limit_value;
    itype64 m00, m11, m22, m01, m02, m12;
  };



  class Grid_range : public Grid
  {
  public:
    inline Grid_range() {}
    Grid_range( const Coord_grid& min, const Coord_grid& max );
    Grid_range( const Grid& grid, const Coord_frac& min, const Coord_frac& max );
    Grid_range( const Cell& cell, const Grid& grid, const ftype& radius );
    const Coord_grid& min() const { return min_; }
    const Coord_grid& max() const { return max_; }
    void add_border( const int b );
    bool in_grid( Coord_grid g ) const { return (g.u() >= min_.u() && g.u() <= max_.u() && g.v() >= min_.v() && g.v() <= max_.v() && g.w() >= min_.w() && g.w() <= max_.w()); }

    int index( const Coord_grid& c ) const { return (c - min_).index(*this); } 
    Coord_grid deindex( const int& index ) const { return Coord_grid( *this, index ) + min_; }
  private:
    Coord_grid min_, max_;
  };
  typedef Grid_range Grid_map;



  class Atom
  {
  public:
    Atom() {}
    template<class T> Atom( const T& atom ) : element_(atom.element()), coord_orth_(atom.coord_orth()), u_aniso_orth_(atom.u_aniso_orth()) , occupancy_(atom.occupancy()), u_iso_(atom.u_iso()){}
    const String& element() const { return element_; }
    const Coord_orth& coord_orth() const { return coord_orth_; }
    const ftype& occupancy() const { return occupancy_; }
    const ftype& u_iso() const { return u_iso_; }
    const U_aniso_orth& u_aniso_orth() const { return u_aniso_orth_; }
    void set_element( const String& s );             
    void set_coord_orth( const Coord_orth& s );      
    void set_occupancy( const ftype& s );            
    void set_u_iso( const ftype& s );                
    void set_u_aniso_orth( const U_aniso_orth& s );  

    void transform( const RTop_orth rt );
    bool is_null() const { return coord_orth_.is_null(); }
    static Atom null();
  private:
    String element_;
    Coord_orth coord_orth_;
    U_aniso_orth u_aniso_orth_;
    ftype occupancy_, u_iso_;
  };



  class Atom_list : public std::vector<Atom>
  {
  public:
    Atom_list() {}
    Atom_list( const std::vector<Atom>& list ) : std::vector<Atom>( list ) {}
    template<class T> Atom_list( const T& list ) { for ( int i = 0; i < list.size(); i++ ) push_back( Atom( list[i] ) ); }
  };


  // some template function definitions

  template<class T> RTop_orth::RTop_orth( const T& src, const T& tgt )
  {
    std::vector<Coord_orth> vsrc( src.size() );
    std::vector<Coord_orth> vtgt( tgt.size() );
    for ( int i = 0; i < src.size(); i++ ) vsrc[i] = src[i].coord_orth();
    for ( int i = 0; i < tgt.size(); i++ ) vtgt[i] = tgt[i].coord_orth();
    (*this) = RTop_orth( vsrc, vtgt );
  }

  // some inline function definitions
  HKL HKL::transform( const Symop& op ) const
    { return Coord_reci_frac(*this).transform(op).hkl(); }
  HKL HKL::transform( const Isymop& op ) const
    { return op*(*this); }
  ftype HKL::sym_phase_shift( const Symop& op ) const
    { return -Util::twopi()*( Coord_reci_frac(*this) * op.trn() ); }

  const Coord_grid& Coord_grid::next( const Grid& g )
    { w()++; if ( w() >= g.nw() ) { w() = 0; v()++; if ( v() >= g.nv() ) { v() = 0; u()++; } } return *this; }
  const Coord_grid& Coord_grid::next( const Grid_range& g )
    { w()++; if ( w() > g.max().w() ) { w() = g.min().w(); v()++; if ( v() > g.max().v() ) { v() = g.min().v(); u()++; } } return *this; }
  bool Coord_grid::last( const Grid& g ) const
    { return ( u() >= g.nu() ); }
  bool Coord_grid::last( const Grid_range& g ) const
    { return ( u() > g.max().u() ); }
  int Coord_grid::index( const Grid& g ) const
    { return ( u()*g.nv() + v() )*g.nw() + w(); }
  void Coord_grid::deindex( const Grid& g, const int& index )
    { u() = index/(g.nv()*g.nw()); v() = (index/g.nw()) % g.nv(); w() = (index) % g.nw(); }

  Coord_grid Coord_grid::unit( const Grid_sampling& g ) const
    { return Coord_grid( Util::mod(u(), g.nu()), Util::mod(v(), g.nv()), Util::mod(w(), g.nw()) ); }
  Coord_map Coord_grid::coord_map() const
    { return Coord_map( *this ); }
  Coord_frac Coord_grid::coord_frac( const Grid_sampling& g ) const
    { return Coord_frac( ftype(u())/ftype(g.nu()), ftype(v())/ftype(g.nv()), ftype(w())/ftype(g.nw()) ); }

  ftype HKL::invresolsq( const Cell& cell ) const
    { return cell.metric_reci().lengthsq( Coord_reci_frac( *this ) ); }
  Coord_reci_frac HKL::coord_reci_frac() const
    { return Coord_reci_frac( *this ); }
  Coord_reci_orth HKL::coord_reci_orth( const Cell& cell ) const
    { return coord_reci_frac().coord_reci_orth( cell ); }
  ftype Coord_reci_orth::invresolsq() const
    { return xs()*xs() + ys()*ys() + zs()*zs(); }
  Coord_reci_frac Coord_reci_orth::coord_reci_frac( const Cell& cell ) const
    { return Coord_reci_frac( (*this) * cell.matrix_orth() ); }
  ftype Coord_reci_frac::invresolsq( const Cell& cell ) const
    { return cell.metric_reci().lengthsq( *this ); }
  Coord_reci_orth Coord_reci_frac::coord_reci_orth( const Cell& cell ) const
    { return Coord_reci_orth( (*this) * cell.matrix_frac() ); }

  ftype Coord_orth::lengthsq() const
    { return x()*x()+y()*y()+z()*z(); }
  Coord_frac Coord_orth::coord_frac( const Cell& cell ) const
    { return Coord_frac( cell.matrix_frac() * (*this) ); }
  ftype Coord_frac::lengthsq( const Cell& cell ) const
    { return cell.metric_real().lengthsq( *this ); }
  Coord_orth Coord_frac::coord_orth( const Cell& cell ) const
    { return Coord_orth( cell.matrix_orth() * (*this) ); }
  Coord_map Coord_frac::coord_map( const Grid& g ) const
    { return Coord_map( u()*ftype(g.nu()), v()*ftype(g.nv()), w()*ftype(g.nw()) ); }
  Coord_grid Coord_frac::coord_grid( const Grid& g ) const
    { return Coord_grid( Util::intr(u()*ftype(g.nu())), Util::intr(v()*ftype(g.nv())), Util::intr(w()*ftype(g.nw())) ); }
  Coord_frac Coord_map::coord_frac( const Grid& g ) const
    { return Coord_frac( u()/ftype(g.nu()), v()/ftype(g.nv()), w()/ftype(g.nw()) ); }

} // namespace clipper

#endif