#
#  crystal.rb
#
#  Created by Toshi Nagata.
#  Copyright 2012 Toshi Nagata. All rights reserved.
#
# This program 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 version 2 of the License.
# 
# This program 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.

#  Definition for use in ORTEP
class AtomRef
  def to_adc
    sym = self.symop
    if sym == nil
      idx = self.index + 1
      symcode = 55501
    else
      idx = sym[4] + 1
      symcode = (sym[1] + 5) * 10000 + (sym[2] + 5) * 1000 + (sym[3] + 5) * 100 + sym[0] + 1
    end
    return idx, symcode
  end
end

#  Used in bond_angle_with_sigma
module Math
  def acos_safe(arg)
    if arg <= -1.0
      return PI
    elsif arg >= 1.0
      return 0.0
    else
      return acos(arg)
    end
  end
end

class Molecule

def export_ortep(fp)

  #  Create atom list
  hidden = atom_group { |ap| !is_atom_visible(ap.index) }
  hydrogen = self.show_hydrogens
  expanded = self.show_expanded
  atomlist = atom_group { |ap|
    (ap.element != "H" || hydrogen) &&
    (ap.symop == nil || expanded) &&
    (!hidden.include?(ap.index))
  }

  #  Title
  fp.printf "%-78.78s\n", self.name + ": generated by Molby at " + Time.now.to_s

  #  Cell parameters
  cp = self.cell
  if cp == nil
    cp = [1, 1, 1, 90, 90, 90]
  end
  fp.printf "%9.3f%9.3f%9.3f%9.3f%9.3f%9.3f\n", cp[0], cp[1], cp[2], cp[3], cp[4], cp[5]
  
  #  Symmetry operations
  syms = self.symmetries
  if syms == nil || syms.length == 0
   fp.print "1             0  1  0  0              0  0  1  0              0  0  0  1\n"
  else
    syms.each_with_index { |s, i|
      a = s.to_a
      fp.printf "%s%14g%3g%3g%3g%15g%3g%3g%3g%15g%3g%3g%3g\n", (i == syms.length - 1 ? "1" : " "), a[9], a[0], a[1], a[2], a[10], a[3], a[4], a[5], a[11], a[6], a[7], a[8]
    }
  end

  #  Atoms (all symmetry unique atoms regardless they are visible or not)
  n = 0
  each_atom { |ap|
    break if ap.symop != nil
    fp.printf " %4.4s%22s%9.4f%9.4f%9.4f%9d\n", ap.name, "", ap.fract_x, ap.fract_y, ap.fract_z, 0
    an = ap.aniso
    if an != nil
      fp.printf " %8.5f%9.6f%9.6f%9.6f%9.6f%9.6f%9d\n", an[0], an[1], an[2], an[3], an[4], an[5], 0
    else
      t = ap.temp_factor
      t = 1.2 if t <= 0
      fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", t, 0.0, 0.0, 0.0, 0.0, 0.0, 6
    end
    n += 1
  }
  natoms_tep = n

  #  Special points to specify cartesian axes
  axis, angle = self.get_view_rotation
  tr = Transform.rotation(axis, -angle)
  org = self.get_view_center
  x = org + tr.column(0)
  y = org + tr.column(1)
  tr = self.cell_transform
  if tr
    tr = tr.inverse
  else
    tr = Transform.identity
  end
  org = tr * org
  x = tr * x
  y = tr * y
  fp.printf " CNTR                      %9.4f%9.4f%9.4f        0\n", org.x, org.y, org.z
  fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
  fp.printf " X                         %9.4f%9.4f%9.4f        0\n", x.x, x.y, x.z
  fp.printf " %8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6
  fp.printf " Y                         %9.4f%9.4f%9.4f        0\n", y.x, y.y, y.z
  fp.printf "1%8.3f%9g%9g%9g%9g%9g%9d\n", 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 6

  #  Initialize
  fp.print  "      201\n"

  #  Pen size
  fp.print  "      205        8\n"

  #  Paper size, margin, viewing distance
  fp.print  "      301      6.6      6.6        0      0.0\n"

  #  Sort atoms by symop and index
  acodes = Hash.new
  ig = IntGroup.new   #  Used for collecting contiguous atoms
  each_atom(atomlist) { |ap|
    idx, symcode = ap.to_adc
    acode = symcode * self.natoms + idx - 1
    acodes[acode] = ap.index
    ig.add(acode)
  }
  index2an = []       #  Index in Molecule to Index in ORTEP
  ig.each_with_index { |acode, i|
    index2an[acodes[acode]] = i + 1
  }
  i = 0
  adcs = []
  while (r = ig.range_at(i)) != nil
    s = r.first
    s = (s / self.natoms) + (s % self.natoms + 1) * 100000  #  Rebuild true ADC (atom index is in the upper digit)
    e = r.last
    e = (e / self.natoms) + (e % self.natoms + 1) * 100000
    if s < e
      adcs.push(s)
      adcs.push(-e)
    else
      adcs.push(s)
    end
    i += 1
  end
  k = 0

  #  Atom list
  adcs.each_with_index { |a, i|
    if k == 0
      fp.print "      401"
    end
    fp.printf "%9d", a
    k += 1
    if i == adcs.length - 1 || k == 6 || (k == 5 && i < adcs.length - 2 && adcs[i + 2] < 0)
      fp.print "\n"
      k = 0
    end
  }

  #  Axes
  fp.printf "      501%4d55501%4d55501%4d55501%4d55501%4d55501                 1\n", natoms_tep + 1, natoms_tep + 1, natoms_tep + 2, natoms_tep + 1, natoms_tep + 3
#  fp.print  "      502        1      0.0        2      0.0        3      0.0\n"
  
  #  Autoscale
  fp.print  "      604                               1.538\n"
  
  #  Explicit bonds
  bond_inst = Array.new(6) { [] }   #  Bonds for types 1 to 5
  bonds.each { |b|
    next if !atomlist.include?(b[0]) || !atomlist.include?(b[1])
    #  TODO: determination of bond types should be refined
    an1 = atoms[b[0]].atomic_number
    an2 = atoms[b[1]].atomic_number
    if an1 == 1 || an2 == 1
      btype = 1
    elsif an1 <= 8 && an2 <= 8
      btype = 3
    else
      btype = 5
    end
    bond_inst[btype].push(b[0], b[1])
  }

  #  Output bond specifications
  #  Avoid including too many ADCs in a single 811/821 instruction
  #  (Upper limit is 140. Here we divide at every 36 ADCs)
  output_bonds = proc { |icode|
    bond_inst.each_with_index { |inst, ii|
      next if inst.length == 0
      inst.each_with_index { |b, i|
        if i % 6 == 0
          fp.printf "  %d   %3s", (i >= inst.length - 6 || i % 36 == 30 ? 2 : 1), (i % 36 == 0 ? icode.to_s : "")
        end
        idx, scode = atoms[b].to_adc
        fp.printf "%9d", idx * 100000 + scode
        if i % 6 == 5 || i == inst.length - 1
          fp.print "\n"
          if i == inst.length - 1 || i % 36 == 35
            fp.printf "%21s%3d%12s%6.3f\n", "", ii, "", 0.05
          end
        end
      }
    }
  }

  fp.print "  0  1001     0.02\n"  #  Activate hidden line removal
  output_bonds.call(821)
  
  #  Atom types
  atom_inst = Array.new(5) { IntGroup.new }   #  Atoms for 714 (sphere), 712 (with axes), 711 (shade and axes)
  atomlist.each { |i|
    #  TODO: determination of atom types should be refined
    an1 = atoms[i].atomic_number
    if an1 == 1
      atype = 4
    elsif an1 <= 6
      atype = 2
    else
      atype = 1
    end
    idx, scode = atoms[i].to_adc
    atom_inst[atype].add(idx)
  }
  [4,2,1].each { |ii|
    inst = atom_inst[ii]
    i = 0
    while (r = inst.range_at(i)) != nil
      fp.printf "  1   %3d\n", 710 + ii
      fp.printf "%27s%9d%9d\n", "", r.first, r.last
      i += 1
    end
  }

  output_bonds.call(811)

  #  Close plot
  fp.print "      202\n"
  fp.print "       -1\n"
  
end

def savetep(filename)
  if natoms == 0
	raise MolbyError, "cannot save ORTEP input; the molecule is empty"
  end
  fp = open(filename, "wb")
  export_ortep(fp)
  fp.close
  return true
end

end

#  Best-fit planes
#  Ref. W. C. Hamilton, Acta Cryst. 1961, 14, 185-189
#       T. Ito, Acta Cryst 1981, A37, 621-624

#  An object to described the best-fit plane
#  
#  Contains the plane coefficients and the constant (a, b, c, d for ax + by + cz + d = 0),
#  the error matrix, and the metric tensor.
#
class Molby::Plane
  attr_accessor :molecule, :group, :coeff, :const, :error_matrix, :metric_tensor
  def initialize(mol, group, coeff, const, err, met)
    @molecule = mol
    @group = group
    @coeff = coeff
    @const = const
    @error_matrix = err
    @metric_tensor = met
    self
  end
  def sigma
    [sqrt(@error_matrix[0, 0]), sqrt(@error_matrix[1, 1]), sqrt(@error_matrix[2, 2]), sqrt(@error_matrix[3, 3])]
  end
  def inspect
    s = sprintf("Molby::Plane[\n coeff, const = [[%f, %f, %f], %f],\n", @coeff.x, @coeff.y, @coeff.z, @const)
    s += sprintf(" sigma = [[%10.4e, %10.4e, %10.4e], %10.4e],\n", *self.sigma)
    (0..3).each { |i|
      s += (i == 0 ? " error_matrix = [" : "     ")
      (0..i).each { |j|
        s += sprintf("%12.6e%s", @error_matrix[j, i], (j == i ? (i == 3 ? "],\n" : ",\n") : ","))
      }
    }
    s += sprintf(" molecule = %s\n", @molecule.inspect)
    s += sprintf(" group = %s\n", @group.inspect)
    (0..3).each { |i|
      s += (i == 0 ? " metric_tensor = [" : "     ")
      (0..3).each { |j|
        s += sprintf("%12.6e%s", @metric_tensor[j, i], (j == 3 ? (i == 3 ? "]]\n" : ",\n") : ","))
      }
    }
    s
  end
  def distance(ap)
    if ap.is_a?(AtomRef)
      fr = ap.fract_r
      sig = ap.sigma
    else
      fr = Vector3D[*ap]
      sig = Vector3D[0, 0, 0]
    end
    d = fr.dot(@coeff) + @const
    sig1 = (@coeff.x * sig.x) ** 2 + (@coeff.y * sig.y) ** 2 + (@coeff.z * sig.z) ** 2
    sig2 = LAMatrix.multiply("t", fr, @error_matrix, fr)[0, 0]
    if ap.is_a?(AtomRef) && ap.molecule == @molecule && @group.include?(ap.index)
      #  The atom defines the plane
      sig0 = sig1 - sig2
      sig0 = 0.0 if sig0 < 0.0
    else
      sig0 = sig1 + sig2
    end  
    return d, sqrt(sig0)
  end
  def dihedral(plane)
    e1 = @error_matrix.submatrix(0, 0, 3, 3)
    e2 = plane.error_matrix.submatrix(0, 0, 3, 3)
    m = @metric_tensor.submatrix(0, 0, 3, 3)
    c = plane.coeff
    cos_t = plane.coeff.dot(m * @coeff)
    if cos_t > 1.0
      cos_t = 1.0
    elsif cos_t < -1.0
      cos_t = -1.0
    end
    t = acos(cos_t)
    sig_t = (m * e1).trace + (m * e2).trace
    if sig_t < t * t
      w = 1.0 / sin(t)
      sig_t = w * w * (c.dot(LAMatrix.multiply(m, e1, m, c)) + @coeff.dot(LAMatrix.multiply(m, e2, m, @coeff)))
    end
    t *= 180.0 / PI
    sig_t = sqrt(sig_t) * 180.0 / PI
    return t, sig_t
  end
end

class Molecule

#  Calculate best-fit plane for the given atoms
#  Return value: a Molby::Plane object

def plane(group)

  #  Number of atoms
  dim = group.length

  #  Positional parameters and standard deviations
  x = []
  sig = []
  sig_min = 1e10
  each_atom(group) { |ap|
    x.push(ap.fract_r)
    sig.push(ap.sigma)
    if (s = ap.sigma_x) > 0.0 && s < sig_min
      sig_min = s
    end
    if (s = ap.sigma_y) > 0.0 && s < sig_min
      sig_min = s
    end
    if (s = ap.sigma_z) > 0.0 && s < sig_min
      sig_min = s
    end
  }
  if sig_min == 1e10
    sig_min = 1e-12
  end
  sig.each { |s|
    s.x = sig_min if s.x < sig_min
    s.y = sig_min if s.y < sig_min
    s.z = sig_min if s.z < sig_min
  }

  #  The metric tensor of the reciprocal lattice
  #  g[j, i] = (ai*).dot(aj*), where ai* and aj* are the reciprocal axis vectors
  t = self.cell_transform
  if t.nil?
    t = Transform.identity
  end
  g2inv = LAMatrix[t]
  g2 = g2inv.inverse
  g2[3, 3] = 0.0
  g2inv[3, 3] = 0.0
  g = LAMatrix.multiply("t", g2, g2)

  #  The variance-covariance matrices of the atomic parameters
  #  mm[k][n] is a 3x3 matrix describing the correlation between the atoms k and n,
  #  and its components are defined as: sigma_k[i] * sigma_n[j] * corr[k, i, n, j],
  #  where corr(k, i, n, j) is the correlation coefficients between the atomic parameters
  #  k[i] and n[j].
  mm = Array.new(dim) { Array.new(dim) }
  zero = LAMatrix.zero(3, 3)
  dim.times { |k|
    dim.times { |n|
      mkn = LAMatrix.new(3, 3)
      if k == n
        3.times { |i|
          3.times { |j|
            if i == j
              mkn[j, i] = sig[k][i] * sig[n][j]
            else
              #  Inter-coordinate correlation should be implemented here
            end
          }
        }
      else
        #  Inter-atomic correlation should be implemented here
      end
      mm[k][n] = (mkn == zero ? zero : mkn)
    }
  }

  #  The variance-covariance matrix of the atom-plance distances
  #  m[j, i] = v.transpose * mm[i][j] * v, where v is the plane coefficient vector
  #  The inverse of m is the weight matrix
  m = LAMatrix.new(dim, dim)
  
  #  The matrix representation of the atomic coordinates
  #  y[j, i] = x[i][j] (for j = 0..2), -1 (for j = 3)
  #  y * LAMatrix[a, b, c, d] gives the atom-plane distances for each atom
  y = LAMatrix.new(4, dim)
  dim.times { |i|
    y[0, i] = x[i].x
    y[1, i] = x[i].y
    y[2, i] = x[i].z
    y[3, i] = 1.0
  }

  #  The coefficients to be determined
  n0 = LAMatrix[1, 1, 1, 0]
  v = LAMatrix[1, 1, 1]     #  The coefficient part

  iter = 0
  while iter < 20

    iter += 1

    #  Set zero to the "constant" part, and normalize the "coefficient" part
    n0[0, 3] = 0.0
    n0 = g2 * n0
    n0.multiply!(1.0 / n0.fnorm)
    n0 = g2inv * n0
    3.times { |i| v[0, i] = n0[0, i] }

    #  Build the variance-covariance matrix    
    dim.times { |i|
      dim.times { |j|
        m[j, i] = LAMatrix.multiply("t", v, mm[i][j], v)[0, 0]
      }
    }
    c = LAMatrix.multiply("t", y, "i", m, y)

    #  Invert c: only the inverse is used in the following, so c is inversed destructively
    cinv = c.inverse!
 
    if iter == 1

      #  Determine the tentative solution, which is given by the eigenvector of cinv * g
      #  for the largest eigenvalue
      evals, evecs = (cinv * g).eigenvalues
      4.times { |i| n0[0, i] = evecs[3, i] }

    else

      #  Convert the coefficient vector to the reciprocal space
      h = g * n0
      
      #  Determine multiplier
      #  In this implementation, the sign of delta-n is opposite from that used in
      #  the reference
      lam = 1.0 / (LAMatrix.multiply("t", h, cinv, h)[0, 0])
      
      #  Solve the linearized equation
      #  (Is the equation 21 in the reference really correct? Shouldn't it read
      #   B = 1 - lambda * C.inverse * H* ? )
      b = LAMatrix.multiply(lam, cinv, g)
      b.sub!(LAMatrix.identity(4))

      dn = b * n0
      n0 += dn

      break if dn[0, 0] ** 2 + dn[0, 1] ** 2 + dn[0, 2] ** 2 < 1e-9

    end
  end

  #  Error matrix = b * cinv * b.transpose
  em = LAMatrix.multiply(b, cinv, "t", b)
  coeff = Vector3D[n0[0, 0], n0[0, 1], n0[0, 2]]
  const = n0[0, 3]

  return Molby::Plane.new(self, group, coeff, const, em, g)

end

def cmd_plane
  plane_settings = @plane_settings || Hash.new
  mol = self
  h = Dialog.run("Best-Fit Planes", "Close", nil) {
    refresh_proc = proc { |it|
      n = it[:tag][/\d/].to_i
      g = plane_settings["group#{n}"]
      if g
        str = g.inspect.sub!("IntGroup[", "").sub!("]", "")
        set_value("group#{n}", str)
        if n == 1 || n == 2
          p = mol.plane(g) rescue p = nil
          plane_settings["plane#{n}"] = p
          if p
            coeff = p.coeff
            const = p.const
            sig = p.sigma
            aps = (n == 1 ? "" : "'")
            str = sprintf("a%s = %f(%f)\nb%s = %f(%f)\nc%s = %f(%f)\nd%s = %f(%f)",
                          aps, coeff.x, sig[0],
                          aps, coeff.y, sig[1],
                          aps, coeff.z, sig[2],
                          aps, const, sig[3])
            set_value("result#{n}", str)
          else
            set_value("result#{n}", "")
          end
          p1 = plane_settings["plane1"]
          p2 = plane_settings["plane2"]
          if p1 && p2
            t, sig = p1.dihedral(p2)
            str = sprintf("%f(%f)", t, sig)
            set_value("dihedral", str)
          else
            set_value("dihedral", "")
          end
        else
          p = plane_settings["plane1"]
          if p
            str = ""
            mol.each_atom(g) { |ap|
              d, sig = p.distance(ap)
              str += sprintf("%d %f(%f)\n", ap.index, d, sig)
            }
            str.chomp!
          else
            str = ""
          end
          set_value("result#{n}", str)
        end
      else
        set_value("group#{n}", "")
        set_value("result#{n}", "")
      end
    }
    set_proc = proc { |it|
      n = it[:tag][/\d/].to_i
      sel = mol.selection
      if sel.count > 0
        str = sel.inspect.sub!("IntGroup[", "").sub!("]", "")
        set_value("group#{n}", str)
        plane_settings["group#{n}"] = sel
      else
        plane_settings["group#{n}"] = nil
      end
      refresh_proc.call(it)
    }
    text_proc = proc { |it|
      n = it[:tag][/\d/].to_i
      str = it[:value].gsub(/[^-.,0-9]/, "")  #  Remove unsane characters
      g = eval("IntGroup[#{str}]") rescue g = nil
      plane_settings["group#{n}"] = g
      refresh_proc.call(it)
    }
    layout(3,
      item(:text, :title=>"Plane 1 (ax + by + cz + d = 0)"),
      -1, -1,
      item(:text, :title=>"Atoms"),
      item(:textfield, :width=>240, :height=>32, :tag=>"group1", :action=>text_proc),
      item(:button, :title=>"Set Current Selection", :tag=>"button1", :action=>set_proc),
      item(:text, :title=>"Results"),
      item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result1"),     
      item(:button, :title=>"Recalculate", :tag=>"refresh1", :action=>refresh_proc),
      item(:line),
      -1, -1,
      item(:text, :title=>"Plane 2 (a'x + b'y + c'z + d' = 0)"),
      -1, -1,
      item(:text, :title=>"Atoms"),
      item(:textfield, :width=>240, :height=>32, :tag=>"group2", :action=>text_proc),
      item(:button, :title=>"Set Current Selection", :tag=>"button2", :action=>set_proc),
      item(:text, :title=>"Results"),
      item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result2"),
      item(:button, :title=>"Recalculate", :tag=>"refresh2", :action=>refresh_proc),
      item(:text, :title=>"Dihedral angle with Plane 1"), -1, -1,
      -1,
      item(:textfield, :width=>240, :height=>16, :tag=>"dihedral"), -1,
      item(:line),
      -1, -1,
      item(:text, :title=>"Distance from Plane 1"), -1, -1,
      item(:text, :title=>"Atoms"),
      item(:textfield, :width=>240, :height=>32, :tag=>"group3", :action=>text_proc),
      item(:button, :title=>"Set Current Selection", :tag=>"button3", :action=>set_proc),
      item(:text, :title=>"Results"),
      item(:textview, :width=>240, :height=>68, :editable=>false, :tag=>"result3"),
      item(:button, :title=>"Recalculate", :tag=>"refresh3", :action=>refresh_proc)
    )
    refresh_proc.call(item_with_tag("refresh1"))
    refresh_proc.call(item_with_tag("refresh2"))
    refresh_proc.call(item_with_tag("refresh3"))
  }
  @plane_settings = plane_settings
end

#  Calculate bond length and angles with standard deviations
#  args can be a single IntGroup (calculate bonds and angles including those atoms)
#  or arrays of 2 or 3 integers (explicitly specifying bonds and angles)
def bond_angle_with_sigma(*args)

  if args.length >= 2 || (args.length == 1 && !args[0].is_a?(IntGroup))
    #  Bonds and angles are explicitly specified
    bonds = []
    angles = []
    args.each { |arg|
      if arg.length == 2 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer)
        bonds.push(arg)
      elsif arg.length == 3 && arg[0].is_a?(Integer) && arg[1].is_a?(Integer) && arg[2].is_a?(Integer)
        angles.push(arg)
      else
        raise MolbyError, "Invalid argument #{arg.inspect}"
      end
    }
  else
    if args.length == 0
	  g = nil
	else
	  g = args[0]
	end
	bonds = self.bonds.select { |b|
	  (g == nil || (g.include?(b[0]) && g.include?(b[1]))) &&
	    (atoms[b[0]].symop == nil || atoms[b[1]].symop == nil)
	}
	angles = self.angles.select { |ang|
	  (g == nil || (g.include?(ang[0]) && g.include?(ang[1]) && g.include?(ang[2]))) &&
	    (atoms[ang[0]].symop == nil || atoms[ang[1]].symop == nil || atoms[ang[2]].symop == nil)
	}
  end

  #  A list of interatomic distance, its differential coefficients, and other
  #  useful quantities.
  #  The hash tag is a list [i, j] (i and j are the atom indices) or [i, j, k]
  #  (i, j, k are the atom indices, and r(ijk) is the distance between the atom i
  #  and the center point between the atoms j and k).
  #  The value is a list of following quantities:
  #    index 0: rij or r(ijk)
  #    index 1-9: d(rij)/d(xij), d(rij)/d(yij), d(rij)/d(zij),
  #      d(rij)/da, d(rij)/db, d(rij)/dc, d(rij)/d(alpha), d(rij)/d(beta), d(rij)/d(gamma)
  #    index 10: the list of the "base atom"
  #    index 11: the list of the transform matrices
  dlist = Hash.new

  #  A list of fractional coordinates and sigmas
  fract = []
  sigma = []
  each_atom { |ap|
    fract.push(ap.fract_r)
    sigma.push(ap.sigma)
  }

  #  A list of base atoms (for symmetry-related atoms) and transform matrices
  bases = []
  trans = []
  symcode = []
  trans_i = Transform.identity
  each_atom { |ap|
    sym = ap.symop
    bases.push(sym ? sym[4] : ap.index)
	if sym
	  tr = transform_for_symop(sym).transpose
      tr[3, 0] = tr[3, 1] = tr[3, 2] = 0.0
	else
	  tr = trans_i
	end
    trans.push(tr)
	symcode.push(sym ? sprintf("%d_%d%d%d", sym[0] + 1, sym[1] + 5, sym[2] + 5, sym[3] + 5) : ".")
  }

  #  Unit cell parameter
  cell = self.cell
  # $a, $b, $c, $alpha, $beta, $gamma, $sig_a, $sig_b, $sig_c, $sig_alpha, $sig_beta, $sig_gamma = self.cell
  cos_a = cos(cell[3] * PI / 180.0)
  cos_b = cos(cell[4] * PI / 180.0)
  cos_c = cos(cell[5] * PI / 180.0)
  abc = cell[0] * cell[1] * cos_c
  bca = cell[1] * cell[2] * cos_a
  cab = cell[2] * cell[0] * cos_b
  aa = cell[0] * cell[0]
  bb = cell[1] * cell[1]
  cc = cell[2] * cell[2]

  get_dlist = proc { |_i, _j, _k|
    if _k != nil
      if _j > _k
        _j, _k = _k, _j
      end
      _p = dlist[[_i, _j, _k]]
      return _p if _p != nil
      _p = (dlist[[_i, _j, _k]] = [])
      _vij = fract[_i] - (fract[_j] + fract[_k]) * 0.5
      _p[10] = [bases[_i], bases[_j], bases[_k]]
      _p[11] = [trans[_i], trans[_j], trans[_k]]
    else
      if _i > _j
        _i, _j = _j, _i
      end
      _p = dlist[[_i, _j]]
      return _p if _p != nil
      _p = (dlist[[_i, _j]] = [])
      _vij = fract[_i] - fract[_j]
      _p[10] = [bases[_i], bases[_j]]
      _p[11] = [trans[_i], trans[_j]]
    end
    _xij = _vij.x
    _yij = _vij.y
    _zij = _vij.z
    _dij = sqrt(aa * _xij * _xij + bb * _yij * _yij + cc * _zij * _zij + 2 * bca * _yij * _zij + 2 * cab * _zij * _xij + 2 * abc * _xij * _yij)
    _p[0] = _dij
    _p[1] = (aa * _xij + abc * _yij + cab * _zij) / _dij
    _p[2] = (bb * _yij + bca * _zij + abc * _xij) / _dij
    _p[3] = (cc * _zij + cab * _xij + bca * _yij) / _dij
    _p[4] = (cell[0] * _xij * _xij + cell[2] * cos_b * _zij * _xij + cell[1] * cos_c * _xij * _yij) / _dij
    _p[5] = (cell[1] * _yij * _yij + cell[0] * cos_c * _xij * _yij + cell[2] * cos_a * _yij * _zij) / _dij
    _p[6] = (cell[2] * _zij * _zij + cell[1] * cos_a * _yij * _zij + cell[0] * cos_b * _zij * _xij) / _dij
    _p[7] = (-cell[1] * cell[2] * sin(cell[3] * PI / 180.0) * _yij * _zij) * (PI / 180.0) / _dij
    _p[8] = (-cell[2] * cell[0] * sin(cell[4] * PI / 180.0) * _zij * _xij) * (PI / 180.0) / _dij
    _p[9] = (-cell[0] * cell[1] * sin(cell[5] * PI / 180.0) * _xij * _yij) * (PI / 180.0) / _dij
    return _p
  }

  diff_by_rn = proc { |_dl, _n, _ijk|
    #  dl = dlist(i, j)
    #  return value: Vector3D[ d(rij)/d(xn), d(rij)/d(yn), d(rij)/d(zn) ]
    #  If ijk is true, then dl is dlist(i, j, k)
    _dv = Vector3D[_dl[1], _dl[2], _dl[3]]
    _c = Vector3D[0, 0, 0]
    if _dl[10][0] == _n
      _c += _dl[11][0] * _dv
    end
    if _ijk
      if _dl[10][1] == _n
        _c -= _dl[11][1] * _dv * 0.5
      end
      if _dl[10][2] == _n
        _c -= _dl[11][2] * _dv * 0.5
      end
    else
      if _dl[10][1] == _n
        _c -= _dl[11][1] * _dv
      end
    end
    return _c
  }

  notate_with_sigma = proc { |_val, _sig|
    if _sig == 0.0
      return sprintf "%.4f", _val
    end
    _lg = log(_sig.abs / 1.95) / log(10.0)
    _n = -(_lg.floor)
    _val2 = (_val * (10 ** _n) + 0.5).floor * (0.1 ** _n)
    return sprintf "%.#{_n}f(%d)", _val2, (_sig * (10 ** _n) + 0.5).floor
  }

  results = []

  c = Vector3D[0, 0, 0]
  bonds.each { |b|
    i = b[0]
    j = b[1]
    if i > j
      i, j = j, i
    end
    p = get_dlist.call(i, j, nil)
    sig = 0.0
    p[10].uniq.each { |k|
      s = sigma[k]
      next unless s
      c = diff_by_rn.call(p, k, false)
      #  Test
      if nil
        apk = atoms[k]
        r = apk.fract_r
        d0 = calc_bond(i, j)
        apk.fract_r = r + Vector3D[0.00001, 0, 0]
        amend_by_symmetry(p[10])
        d1 = calc_bond(i, j)
        apk.fract_r = r + Vector3D[0, 0.00001, 0]
        amend_by_symmetry(p[10])
        d2 = calc_bond(i, j)
        apk.fract_r = r + Vector3D[0, 0, 0.00001]
        amend_by_symmetry(p[10])
        d3 = calc_bond(i, j)
        apk.fract_r = r
        amend_by_symmetry(p[10])
        printf " dr/dv[%d] cal %10.5g %10.5g %10.5g\n", k, c[0], c[1], c[2]
        printf " dr/dv[%d] est %10.5g %10.5g %10.5g\n", k, (d1 - d0) / 0.00001, (d2 - d0) / 0.00001, (d3 - d0) / 0.00001
      end
      sig += c[0] * c[0] * s[0] * s[0] + c[1] * c[1] * s[1] * s[1] + c[2] * c[2] * s[2] * s[2]
    }
    6.times { |n|
      sig += (p[4 + n] * cell[6 + n]) ** 2
    }
    #  Test
    if nil
      ncell = cell.dup
      d0 = calc_bond(i, j)
      dd = [0, 0, 0, 0, 0, 0]
      6.times { |n|
        ncell[n] += 0.0001
        set_cell(ncell, true)
        amend_by_symmetry(p[10])
        d1 = calc_bond(i, j)
        dd[n] = (d1 - d0) / 0.0001
        ncell[n] = cell[n]
      }
      set_cell(ncell, true)
      amend_by_symmetry(p[10])
      printf " dr/d{a,b,c} cal %10.5g %10.5g %10.5g\n", p[4], p[5], p[6]
      printf " dr/d{a,b,c} est %10.5g %10.5g %10.5g\n", dd[0], dd[1], dd[2]
      printf " dr/d{alpha,beta,gamma} cal %10.5g %10.5g %10.5g\n", p[7], p[8], p[9]
      printf " dr/d{alpha,beta,gamma} est %10.5g %10.5g %10.5g\n", dd[3], dd[4], dd[5]
    end
    sig = sqrt(sig)
    results.push([i, j, nil, notate_with_sigma.call(p[0], sig), symcode[i], symcode[j]])
  }

  angles.each { |ang|
    i = ang[0]
    j = ang[1]
    k = ang[2]
    p0 = get_dlist.call(i, j, nil)
    p1 = get_dlist.call(j, k, nil)
    p2 = get_dlist.call(i, k, nil)
    p3 = get_dlist.call(j, i, k)
    t123 = acos_safe((p0[0] ** 2 + p1[0] ** 2 - p2[0] ** 2) / (2 * p0[0] * p1[0]))
    t124 = acos_safe((p0[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p0[0] * p3[0]))
    t324 = acos_safe((p1[0] ** 2 + p3[0] ** 2 - p2[0] ** 2 * 0.25) / (2 * p1[0] * p3[0]))
    if nil
      printf "t123 = %.2f t124+t324 = %.2f t124 = %.2f t324 = %.2f\n", t123*180/PI, (t124+t324)*180/PI, t124*180/PI, t324*180/PI
    end
    dtdr12 = -(p0[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p0[0] ** 2) * sin(t124))
    dtdr23 = -(p1[0] ** 2 - p3[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * p3[0] * (p1[0] ** 2) * sin(t324))
    dtdr13 = p2[0] / (sin(t124) * 4 * p0[0] * p3[0]) + p2[0] / (sin(t324) * 4 * p1[0] * p3[0])
    dtdr24 = -(p3[0] ** 2 - p0[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * (p3[0] ** 2) * p0[0] * sin(t124)) - (p3[0] ** 2 - p1[0] ** 2 + p2[0] ** 2 * 0.25) / (2 * (p3[0] ** 2) * p1[0] * sin(t324))
    pp = (p0[10] + p1[10] + p2[10] + p3[10]).uniq
    sig = 0.0
    pp.each { |n|
      s = sigma[n]
      next unless s
      c = Vector3D[0, 0, 0]
      c += diff_by_rn.call(p0, n, false) * (dtdr12 * 180.0 / PI)
      c += diff_by_rn.call(p1, n, false) * (dtdr23 * 180.0 / PI)
      c += diff_by_rn.call(p2, n, false) * (dtdr13 * 180.0 / PI)
      c += diff_by_rn.call(p3, n, true) * (dtdr24 * 180.0 / PI)
      #  Test
      if nil
        apn = atoms[n]
        r = apn.fract_r
        t0 = calc_angle(i, j, k)
        apn.fract_r = r + Vector3D[0.00001, 0, 0]
        amend_by_symmetry(pp)
        t1 = calc_angle(i, j, k)
        apn.fract_r = r + Vector3D[0, 0.00001, 0]
        amend_by_symmetry(pp)
        t2 = calc_angle(i, j, k)
        apn.fract_r = r + Vector3D[0, 0, 0.00001]
        amend_by_symmetry(pp)
        t3 = calc_angle(i, j, k)
        apn.fract_r = r
        amend_by_symmetry(pp)
        printf " dt/dv[%d] cal %10.5g %10.5g %10.5g\n", n, c[0], c[1], c[2]
        printf " dt/dv[%d] est %10.5g %10.5g %10.5g\n", n, (t1 - t0) / 0.00001, (t2 - t0) / 0.00001, (t3 - t0) / 0.00001
      end
      sig += c[0] * c[0] * s[0] * s[0] + c[1] * c[1] * s[1] * s[1] + c[2] * c[2] * s[2] * s[2]
    }
    dd = dtdr12 * p0[4] + dtdr23 * p1[4] + dtdr13 * p2[4] + dtdr24 * p3[4]
    sig += dd * dd * cell[6] * cell[6]
    dd = dtdr12 * p0[5] + dtdr23 * p1[5] + dtdr13 * p2[5] + dtdr24 * p3[5]
    sig += dd * dd * cell[7] * cell[7]
    dd = dtdr12 * p0[6] + dtdr23 * p1[6] + dtdr13 * p2[6] + dtdr24 * p3[6]
    sig += dd * dd * cell[8] * cell[8]
    dd = dtdr12 * p0[7] + dtdr23 * p1[7] + dtdr13 * p2[7] + dtdr24 * p3[7]
    sig += dd * dd * cell[9] * cell[9]
    dd = dtdr12 * p0[8] + dtdr23 * p1[8] + dtdr13 * p2[8] + dtdr24 * p3[8]
    sig += dd * dd * cell[10] * cell[10]
    dd = dtdr12 * p0[9] + dtdr23 * p1[9] + dtdr13 * p2[9] + dtdr24 * p3[9]
    sig += dd * dd * cell[11] * cell[11]
    sig = sqrt(sig)
    results.push([i, j, k, notate_with_sigma.call(t123*180/PI, sig), symcode[i], symcode[j], symcode[k]])
  }
  results
end

def cmd_bond_angle_with_sigma
  if self.cell == nil
    error_message_box "Unit cell is not defined"
  elsif self.cell.length < 12
    error_message_box "Sigmas are not assigned for the unit cell parameters"
  else
    if selection.count == 0
	  a = bond_angle_with_sigma
	else
	  a = bond_angle_with_sigma(selection)
	end
	s = ""
	a.each { |e|
	  if e[2] == nil
	    ss = sprintf("%s %s  %s  %s %s\n", atoms[e[0]].name, atoms[e[1]].name, e[3], e[4], e[5])
	  else
	    ss = sprintf("%s %s %s  %s  %s %s %s\n", atoms[e[0]].name, atoms[e[1]].name, atoms[e[2]].name, e[3], e[4], e[5], e[6])
	  end
	  s += ss
	}
	print s
  end
end

if lookup_menu("Best-fit Planes...") < 0
  register_menu("", "")
  register_menu("Best-fit Planes...", :cmd_plane)
  register_menu("Bonds and angles with sigma", :cmd_bond_angle_with_sigma)
end

end
