#! /usr/bin/python3 -t 
# _*_ coding: iso-8859-1 _*_

MODULE_NAME = "mformula_hydrocarbon"
MODULE_DESC = "general routines for building hydrocarbon formulas"
MODULE_VERS = "1.0"

MODULE_COPYRIGHT = "Copyright © 2013-02-14 by the State University of Campinas (UNICAMP)"

import sys
import os
import re
import math; from math import sqrt,sin,cos,pi
# sys.path[1:0] = [ sys.path[0] + '/../lib', os.path.expandvars('${STOLFIHOME}/lib'), '.' ]
# sys.stderr.write(re.sub('[,]', ',\n', "%s: path = %r\n" % (PROG_NAME, sys.path)));
import rn
import mformula

def build_formula_radial_C_CHm_n(style,n,m,chc,chps,vps,den,subs) :
  """Builds the formula for an equilibrated radial C(CH[m])[n] molecule, radical or ion.
  
  The molecule drawing has radial {n}-fold symmetry; if {n} is even
  member 0 points to +X, if {n} is odd it points to +Y. The parameter
  {m} is the number of filled valences per peripheral carbon. The
  parameters {chc,den} are integers, {chps,vps} are lists of integers;
  the charge in the central atom is {chc/den}, each peripheral group {i}
  is connected by a bond with valence {vps[i]/den} and has charge
  {chps[i]/den}. 
  
  The free bonds are computed from {n,m,ch}. The
  parameter {subs} is a string with {m} letters, describing the
  substituents CCW on each radial carbon; currently only 'H' for
  hydrogen is accepted."""

  fm = mformula.obj(style);
  assert n >= 2;
  assert m == len(subs);
  
  # Compute radial bond valences {vp[i]}, total bond valence {vt}, periph charges {qp[i]}:
  vp = [None] *n;
  qp = [None] *n;
  up = [None] *n;
  vt = 0;
  for i in range(n) :
    vp[i] = (vps[i] + 0.0)/den;
    vt = vt + vps[i];
    qp[i] = (chps[i] + 0.0)/den;
    up[i] = 4 + (chps[i] - vps[i] + 0.0)/den - m;
    assert up[i] >= 0;
  vt = vt/den;
    
  # Compute central charge {qc} as float and dangling valence at center {uc}:
  qc = (chc + 0.0)/den;
  uc = 4 + qc -  vt; 
  assert uc >= 0;
  
  # Compute angle of group 0 {ang0} and angle increment {angs} (degrees):
  ang0 = (n % 2)*90.0;
  angs = 360.0/n;
  
  p0 = [0,0];
  k0 = fm.add_atom("C", p0, 0,0);
  if uc > 0 :
    # Add dangling valence on central atom:
    angu = ang0 + 0.5*angs;
    fm.add_free_bond(fm, k0,angu*pi/180,uc); 

  kk = [None] * n;     # {kk[i]} is the index of the C atom of group {i}.
  
  # Add the radial groups:
  sys.stderr.write("! build_formula_radial_C_CHm_n  n = %r  m = %r  qc = %r  uc = %r\n" % (n,m,qc,uc));
  for i in range(n) :
    sys.stderr.write("!   i = %d  qp = %r  vp = %r  up = %r\n" % (i,vp[i],qp[i],up[i]));

    # Compute distance {R} of peripheral carbon:
    dp = fm.bond_length(vp[i]);
    R = dp;

    # Add formula {fg} of the spoke:
    fg = build_formula_radial_C_CHm_n_elem(style,R, qp[i],up[i],m,subs);
    kk[i] = fm.add_subformula(fg, 0, ang0 + i*angs, [0,0]);
    
    # Add radial bonds with valency {vb}:
    fm.add_bond(k0+0,kk[i]+0,vp[i]);

  return fm;
  #----------------------------------------------------------------------
  
def build_formula_radial_C_CHm_n_elem(style,R, q,u,m,subs) :
  """Builds the formula for one radial unit of a radial hydrocarbon
  molecule, ion, or radical.
  
  The group is orizontal with carbon at {[R,0]}; it will have a charge
  {q}, unfilled valency {u}, and {m} substituents with valency 1."""
  
  fm = mformula.obj(style);
  # Compute length {db} of bond to substituent and angle increment {angs} (degrees):
  db = fm.bond_length(1);
  angs = 180/(m - 0.5);
  iC = fm.add_atom("C", [R, 0], q,0);
  for j in range(m) :
    angj = -90 + (j + 0.25)*angs; # angle of substituent in degrees
    ca = cos(angj*pi/180); sa = sin(angj*pi/180);
    iH = fm.add_atom(subs[j], [R + db*ca, 0 + db*sa], 0,0);
    fm.add_bond(iC,iH,1.0);
  if u > 0 :
    # Add dangling valence on peripheral carbon:
    angu =  -90 + 0.75*angs;
    fm.add_free_bond(iC,angu*pi/180,u);
    
  return fm;
  #----------------------------------------------------------------------

def compute_valence_for_radial_C_CHm_n(n,m,qc,qp) :
  """Computes the radial bond valency for {build_formula_radial_C_CHm_n}.
  
  Assumes {n} peripheral carbons with {m} substituents each, charge {qc}
  on central carbon, charge {qp} on each peripheral carbon. """
  
  return min((4.0 + qc)/n, 4 + qp - m);
  #----------------------------------------------------------------------