#! /usr/bin/python3 -t

MODULE_NAME = "mformula_carbon_suboxide"
MODULE_DESC = "general routines for building carbon suboxide polymer formulas"
MODULE_VERS = "1.0"

MODULE_COPYRIGHT = "Copyright © 2009-07-11 by the State University of Campinas (UNICAMP)"

import sys
import os
import re
import copy
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_poly_carbon_suboxide_monomer(svg,Lsym,Msym,Rsym,Lang,Rang) :
  """Builds the formula for one full monomer unit of poly(carbon suboxide).
  
  The monomer consists of two half-units. A copy {A} of the monomer will
  automatically connect to a copy {B} that has been displaced in X by
  3.0 single-bond lengths. The dangling half-bonds will be terminated by
  'atoms' of type {Lsym} on the left and {Rsym} on the right. The bonds
  between the two hemimers will be marked by 'atoms' of type {Msym}. The
  dangling bonds are omitted if {Lsym,Msym,Rsym} are {None}, and are
  left uncapped if {Lsym,Msym,Rsym} are '.'. """

  fm = mformula.obj();

  sb = svg.rel_bond_length(1.0); # Relative length of single bonds.
  dx = 1.5*sb;              # X displacement between hemimers.

  fL = build_formula_poly_carbon_suboxide_hemimer(svg,Lsym,Msym,Lang,0);
  fR = build_formula_poly_carbon_suboxide_hemimer(svg,Msym,Rsym,0,Rang);
  
  m0 = fm.add_subformula(fL, 0,   0,  [0,0]);
  m1 = fm.add_subformula(fR, 1, 180,  [dx,0]);
  return fm;
  #----------------------------------------------------------------------
 
def build_formula_poly_carbon_suboxide_hemimer(svg,Lsym,Rsym,Lang,Rang) :
  """Builds the formula for one half-unit of poly(carbon suboxide).
  
  The unit is a zigzag-bent C3O2 molecule, roughly vertical from top to
  bottom. From top to bottom, the atoms are numbered O4,C2,C0,C1, and O3
  (don't ask why), connected as O=C-C=C-O. The X axis bisects the double
  C=C bond between atoms 0 and 1. That bond is centered on the origin
  and tilted so that the horizontal displacement between them C0 and C1
  is half a single bond length.
  
  There are dangling single half-bonds exiting from atoms C2 and C1
  towards the left, and from C0 and O3 towards the right.
  
  If {Lang} and {Rang} are zero, a copy {A} of this unit will
  automatically connect to a copy {B} that is displaced in X by 1.5
  single bond lengths and Y-flipped. The dangling half-bonds will be
  terminated by 'atoms' of type {Lsym} on the left and {Rsym} on the
  right. The dangling bonds are omitted if {Lsym,Rsym} are {None}." \ "
  If {Lang} is nonzero, the top O=C- unit (atoms O4 and C2) will be
  rotated by {+Lang} radians, while the left-pointing dangling bond on
  C0 will be rotated by {-Lang}. The left-pointing dangling bond on C2
  will be rotated aditionally for proper matching with the terminator
  unit, as defined by {build_formula_poly_carbon_suboxide_terminator}."
  \ " Analogously, if {Rang} is nonzero, the bottom -O (atom O3) will be
  rotated by {+Rang}, while the right-pointing dangling bond on C0 will
  be rotated {-Lang}. The right-pointing dangling bond on O3 will be
  rotated aditionally for proper matching with a terminator unit rotated
  180 degrees. """
  
  show_box = 0;
  
  fu = mformula.obj();  

  dCC1 = svg.rel_bond_length(1.0); # Relative length of single bonds.
  dCC2 = svg.rel_bond_length(2.0); # Relative length of double C-C bonds.
  dCO2 = svg.CObondlength*svg.rel_bond_length(2.0); # Relative length of double C-O bonds.

  # All angles below are measured clockwise from the +X axis.

  ah = 120*pi/180;    # Unbent angle of C0-C2 and O3-C1 bonds.

  # The middle C0=C1 bond is tilted so that its dx is that of a single bond at 6 degrees: 
  a0 = math.acos(-dCC1*cos(ah)/dCC2);   # Angle of double C0=C1 bond.
  
  # The bottom left dangling bond on C1 is tilted from horizontal by {Lang}:
  a5 = pi - Lang;          # Angle of left dangling bond on atom C1.
  
  # Compute the angle {a2} of C0-C2 bond and angle {a6} of dangling bond C2-X6:
  if Lang != 0 :
    # Assume a pentagon with one shorter side and two distinct angles {Lalfa,Lbeta}:
    Lbeta = a5 - a0;        # Internal angle adjacent to short side.
    Lalfa = pi - 2*Lbeta/3.0; # Internal angle between two long sides.
    a2 = pi - Lbeta + a0;   # Angle of C0-C2 bond.
    a6 = pi + a2 - Lalfa;   # Angle of C2-X6 bond
  else :
    # Assume a X-symmetric hexagon with 60 degree angles at top:
    a2 = ah;  # Angle of C0-C2 bond.
    a6 = pi;  # Angle of C2-X6 bond
  
  # The top right dangling bond C0-X7 is tilted from horizontal by {Rang}:
  a7 = - Rang;       # Angle of dangling bond C0-X7.
  
  # Compute the angles {a3} of C1-O3 bond and {a8} of dangling bond O3-X8:
  # The bottom C1-O3 bond is tilted by {Rang} from the basic 60 degrees:
  if Rang != 0 :
    # Assume a pentagon with one shorter side and two distinct angles {Ralfa,Rbeta}:
    Rbeta = pi + a7 - a0;   # Internal angle adjacent to short side.
    Ralfa = pi - 2*Rbeta/3; # Internal angle between two long sides.
    a3 = a0 - Rbeta;        # Angle of C1-O3 bond.
    a8 = pi + a3 - Ralfa;   # Angle of O3-X8 bond
  else :
    # Assume a X-symmetric hexagon with 60 degree angles at bottom:
    a3 = ah - pi; # Angle of C1-O3 bond.
    a8 = 0;       # Angle of O3-X8 bond
    


  # Choose direction of C=O in top carbonyl to bisect the corner angle:
  a4 = (a2 + a6 - pi)/2.0;     # Angle of C2=O4 bond.

  # Atom positions:
  p0 = [ +dCC2*cos(a0)/2.0, +dCC2*sin(a0)/2.0 ];         # Middle carbon.
  p1 = [ -dCC2*cos(a0)/2.0, -dCC2*sin(a0)/2.0 ];         # Lower carbon.
  p2 = rn.add(p0, [ +dCC1*cos(a2), +dCC1*sin(a2) ]); # Upper carbon.
  p3 = rn.add(p1, [ +dCC1*cos(a3), +dCC1*sin(a3) ]); # Lower (ester) oxygen.
  p4 = rn.add(p2, [ +dCO2*cos(a4), +dCO2*sin(a4) ]); # Upper (carbonyl) oxygen.
  
  p5 = rn.add(p1, [ +dCC1*cos(a5)/2.0, +dCC1*sin(a5)/2.0 ]);  # Bot left dangling bond end.
  p6 = rn.add(p2, [ +dCC1*cos(a6)/2.0, +dCC1*sin(a6)/2.0 ]);  # Top left dangling bond end.
  p7 = rn.add(p0, [ +dCC1*cos(a7)/2.0, +dCC1*sin(a7)/2.0 ]);  # Top right dangling bond end.
  p8 = rn.add(p3, [ +dCC1*cos(a8)/2.0, +dCC1*sin(a8)/2.0 ]);  # Bot right dangling bond end.
  
  # Build a formula {fu} for the unit:
  k0 = fu.add_atom("C", p0, 0,0);
  k1 = fu.add_atom("C", p1, 0,0);
  k2 = fu.add_atom("C", p2, 0,0);
  k3 = fu.add_atom("O", p3, 0,0);
  k4 = fu.add_atom("O", p4, 0,0);
  
  fu.add_bond(k0,k1,2.0);
  fu.add_bond(k0,k2,1.0);
  fu.add_bond(k1,k3,1.0);
  fu.add_bond(k2,k4,2.0);
  
  if show_box :
    # Add a square bounding box and the origin:
    kbmm = fu.add_atom(".", [ -0.75*dCC1, -3.0*dCC1 ], 0,0);
    kbpm = fu.add_atom(".", [ +0.75*dCC1, -3.0*dCC1 ], 0,0);
    kbpp = fu.add_atom(".", [ +0.75*dCC1, +3.0*dCC1 ], 0,0);
    kbmp = fu.add_atom(".", [ -0.75*dCC1, +3.0*dCC1 ], 0,0);
    fu.add_bond(kbmm,kbmp,0.5);
    fu.add_bond(kbmp,kbpp,0.5);
    fu.add_bond(kbpp,kbpm,0.5);
    fu.add_bond(kbpm,kbmm,0.5);
    kbmp = fu.add_atom("H", [ 0, 0 ], 0,0);
  
  if Lsym != None :
    # Left dangling bonds:
    k5 = fu.add_atom(Lsym, p5, 0,0);
    k6 = fu.add_atom(Lsym, p6, 0,0);
    fu.add_bond(k1,k5,1.0);
    fu.add_bond(k2,k6,1.0);

  if Rsym != None :
    # Right dangling bonds:
    k7 = fu.add_atom(Rsym, p7, 0,0);
    k8 = fu.add_atom(Rsym, p8, 0,0);
    fu.add_bond(k0,k7,1.0);
    fu.add_bond(k3,k8,1.0);
  
  return fu;
  #----------------------------------------------------------------------

def build_formula_poly_carbon_suboxide_terminator(svg,Rsym,Rang,Rinv) :
  """Builds the formula for a terminator unit of poly(carbon suboxide).
  
  The unit is a boomerang-bent C3O2 molecule, with arms pointing roughly
  NW and SW. We consider first the case when {Rinv} is 0. The atoms are
  numbered from top to bottom O4, C2, C0, C1, and O3 (don't ask why),
  and connected as O4=C2=C0-C1=O3. There are dangling half-length single
  bonds on the right side, namely C0-X5 and C1-X6.
  
  If {Rang} is zero, the bond C0-C1 is vertical and centered on the X
  axis at 1/4 of a bond to the left, and the two dangling bonds are
  horizontal. Then a copy {A} of this unit will automatically connect to
  a copy {B} that is displaced in X by 1.5 single bond lengths and
  either X-flipped or rotated by 180 degrees, creating a cyclobutane
  ring with two ketene and two ketone groups.
  
  If {Rang} is nonzero, the bottom dangling bond C1-X6 is turned {Rang}
  radians from the horizontal. The whole terminator and the bond C0-X5
  are rotated so that the terminator will fit into a hemimer that was
  created with {build_formula_poly_carbon_suboxide_hemimer}, using the
  same angle as its {Lang} parameter, displaced in X by 1.5 single bond
  lengths.
  
  If {Rinv} is 1, the formula is mirrored across a line perpendicular to
  the {C0-C1} bond (so that the ketene is at the bottom). The connecting
  points {X5} and {X6} are exchanged but the unit still connects to the
  next one as before. """

  show_box = 0;
  
  # http://pubs.acs.org/doi/pdf/10.1021/ma60061a015
  
  fu = mformula.obj();

  dCC1 = svg.rel_bond_length(1.0); # Relative length of single bonds.
  dCC2 = svg.rel_bond_length(2.0); # Relative length of double bonds.
  dCO2 = svg.CObondlength*svg.rel_bond_length(2.0); # Relative length of double C-O bonds.

  # All angles below are measured clockwise from the +X axis.

  # Define the angles {a1--a6}:
  if Rang == 0 :
    # Assume a square ring:
    a1 = pi/2;      # Angle of C1-C0.
    a0 = pi/2;      # Angle of ring bond opposite to C1-C0.
    a2 = +0.75*pi;  # Angle of C0=C2.
    a3 = -0.75*pi;  # Angle of C1=O3.
    a4 = a2;        # Angle of C2=O4.
    a5 = 0;         # Angle of C0-X5.
    a6 = 0;         # Angle of C1-X6.
  else :
    # Assume a pentagon with one shorter side (in the next unit) and two distinct angles {Ralfa,Rbeta}:
    # Assume that the short side is tilted CW so that its dX is {dCC1/2}:
    a0 = math.acos(0.5*dCC1/dCC2); # Angle of short side of 5-cycle (in next hemimer).
    a6 = - Rang;               # Angle of C1-X6.
    Rbeta = pi + a6 - a0;      # Internal angle adjacent to short side.
    Ralfa = pi - 2*Rbeta/3.0;  # Internal angle between two long sides.
    a1 = Ralfa + a6;           # Angle of C1-C0 bond.
    a5 = a1 + Ralfa - pi;      # Angle of C0-X5.
    a2 = (a1 - a5 + pi)/2.0;   # Angle of C0=C2.
    a4 = a2;                   # Angle of C2=O4.
    a3 = (a1 + a6)/2.0 + pi;   # Angle of C1=O3.
    
  if Rinv :
    # Flip all the angles (except {a0}) in the direction C0-C1:
    a2 = 2*a1 - a2 - pi;
    a3 = 2*a1 - a3 - pi;
    a4 = 2*a1 - a4 - pi;
    a5 = 2*a1 - a5 - pi;
    a6 = 2*a1 - a6 - pi;
    a1 = a1 - pi;
    
  # Assumed shift between successive units:
  ushift = 1.5*dCC1;
    
  # Compute the position {q1} of the next units's lower ring atom:
  if Rang == 0 :
    # Assume a square upright ring with no double bonds, centered on the boundary between the units:
    assert a0 == pi/2;
    q1 = [ ushift/2.0 + dCC1/2.0, -dCC1/2.0 ];
  else :
    # Assume a pentagon with one double bond (in the next unit) and two distinct angles {Ralfa,Rbeta}:
    # Assume that the double bond is centered on the next unit's origin and tilted by {a0}:
    q1 = [ -dCC2*cos(a0)/2.0 + ushift, -dCC2*sin(a0)/2.0 ];

  # Compute positions {p0--p6}:
  if Rinv:
    # Carbon C0 is connected to {q1}, link angle {a5}:
    p0 = rn.add(q1, [ -dCC1*cos(a5), -dCC1*sin(a5) ]); # Carbon C0.
    p1 = rn.add(p0, [ -dCC1*cos(a1), -dCC1*sin(a1) ]); # Carbon C1.
  else :
    # Carbon C1 is connected to {q1}, link angle {a6}:
    p1 = rn.add(q1, [ -dCC1*cos(a6), -dCC1*sin(a6) ]); # Carbon C1.
    p0 = rn.add(p1, [ +dCC1*cos(a1), +dCC1*sin(a1) ]); # Carbon C0.
    
  p2 = rn.add(p0, [ +dCC1*cos(a2), +dCC1*sin(a2) ]); # Carbon C2.
  p3 = rn.add(p1, [ +dCO2*cos(a3), +dCO2*sin(a3) ]); # Carbonyl oxygen O3.
  p4 = rn.add(p2, [ +dCO2*cos(a4), +dCO2*sin(a4) ]); # Carbonyl oxygen O4.
  
  p5 = rn.add(p0, [ +dCC1*cos(a5)/2.0, +dCC1*sin(a5)/2.0 ]);  # End of top right dangling bond.
  p6 = rn.add(p1, [ +dCC1*cos(a6)/2.0, +dCC1*sin(a6)/2.0 ]);  # End of bot right dangling bond.
  
  # Build a formula {fu} for the unit:
  k0 = fu.add_atom("C", p0, 0,0);
  k1 = fu.add_atom("C", p1, 0,0);
  k2 = fu.add_atom("C", p2, 0,0);
  k3 = fu.add_atom("O", p3, 0,0);
  k4 = fu.add_atom("O", p4, 0,0);
  
  fu.add_bond(k0,k1,1.0);
  fu.add_bond(k0,k2,2.0);
  fu.add_bond(k1,k3,2.0);
  fu.add_bond(k2,k4,2.0);
  
  if show_box :
    # Add a square bounding box and the origin:
    kbmm = fu.add_atom(".", [ -0.75*dCC1, -3.0*dCC1 ], 0,0);
    kbpm = fu.add_atom(".", [ +0.75*dCC1, -3.0*dCC1 ], 0,0);
    kbpp = fu.add_atom(".", [ +0.75*dCC1, +3.0*dCC1 ], 0,0);
    kbmp = fu.add_atom(".", [ -0.75*dCC1, +3.0*dCC1 ], 0,0);
    fu.add_bond(kbmm,kbmp,0.5);
    fu.add_bond(kbmp,kbpp,0.5);
    fu.add_bond(kbpp,kbpm,0.5);
    fu.add_bond(kbpm,kbmm,0.5);
    kbmp = fu.add_atom("H", [ 0, 0 ], 0,0);

  if Rsym != None :
    # Right dangling bonds:
    k5 = fu.add_atom(Rsym, p5, 0,0);
    k6 = fu.add_atom(Rsym, p6, 0,0);
    fu.add_bond(k0,k5,1.0);
    fu.add_bond(k1,k6,1.0);
  
  return fu;
  #----------------------------------------------------------------------