bondangle:=proc(XA1,XA2,XA3)
#
description `This computes the bond angle in degrees between two bonds`,
`{A1,A2} and {A2,A3}, where the positions of these three atoms are`,
`given as column vector aguments.`;
#
evalf(180/Pi*arccos(DotProduct((XA1-XA2)/Norm(XA1-XA2,2),
(XA3-XA2)/Norm(XA3-XA2,2))));
end proc;

wedgeangle:=proc(XA1,XA2,XA3,XA4)
#
local u,V1,v1,V2,v2,z,phi;
#
description `This computes a wedge angle in degrees`,
`in the range 0 to 360 for the wedge between the two triangles`,
`{A1,A2,A3} and {A2,A3,A4} with respect to the orientation`,
`[A1,A2,A3,A4]. The column vectors for the positions of these`,
`four atoms are given as arguments.`;
#
u:=(XA3-XA2)/Norm(XA3-XA2,2);
V1:=XA1-XA2-u*(DotProduct(u,XA1-XA2));
v1:=V1/Norm(V1,2);
V2:=XA4-XA2-u*(DotProduct(u,XA4-XA2));
v2:=V2/Norm(V2,2);
z:=DotProduct(v1,v2)+I*DotProduct(CrossProduct(v1,v2),u);
phi:=evalf(180/Pi*argument(z));
if phi<0. then phi:=phi+360.0 end if;
return phi
end proc;

fiveMD2flap:=proc(b0,b1,b2,b3,b4,q,P,S,Gamma)
#
local bav,x0i,x1i,x2i,x3i,x4i,y0i,y1i,y2i,y3i,y4i,
eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,soln,
c0,c1,c2,c3,c4,c5,i,x0s,x1s,x2s,x3s,x4s,y0s,y1s,y2s,y3s,y4s,
z0s,z1s,z2s,z3s,z4s,XC0,XC1,XC2,XC3,XC4,tau1,tau4,phi4302,phi1203;
#
description `This converts Marzec-Day internal coordinates`,
`for a five membered ring into flap angle internal coordinates.`, 
`b0, .., b4 are the bond lengths in angstroms`,
`q,P are the Cremer-Pople pseudorotation amplitude and phase`,
`S,Gamma are the Marzec-Day elongation amplitude and phase`,
`Besides the bond lengths, the flap angle internal coordinates are:`,
`tau1,tau4: bond angles with vertices C1 and C4 respectively`,
`phi1203,phi4302: flap angles with axes C0-C2 and C0-C3 respectively.`,
`It is assumed that the user has previously typed:`,
`with(LinearAlgebra);.`;
#
bav:=evalf((b0+b1+b2+b3+b4)/5);
x0i:=0.;
y0i:=bav;
x1i:=evalf(bav*cos(Pi/2-2*Pi*1/5));
y1i:=evalf(bav*sin(Pi/2-2*Pi*1/5));
x2i:=evalf(bav*cos(Pi/2-2*Pi*2/5));
y2i:=evalf(bav*sin(Pi/2-2*Pi*2/5));
x3i:=evalf(bav*cos(Pi/2-2*Pi*3/5));
y3i:=evalf(bav*sin(Pi/2-2*Pi*3/5));
x4i:=evalf(bav*cos(Pi/2-2*Pi*4/5));
y4i:=evalf(bav*sin(Pi/2-2*Pi*4/5));
eqn1:=x0+x1+x2+x3+x4=0;
eqn2:=y0+y1+y2+y3+y4=0;
eqn3:=x0=0;
eqn4:=S^2*cos(2*Gamma)=x0^2-y0^2+x1^2-y1^2+x2^2-y2^2+x3^2-y3^2+x4^2-y4^2;
eqn5:=S^2*sin(2*Gamma)=2*(x0*y0+x1*y1+x2*y2+x3*y3+x4*y4);
c0:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*0/5+2*Pi/5))^2;
c1:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*1/5+2*Pi/5))^2;
c2:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*2/5+2*Pi/5))^2;
c3:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*3/5+2*Pi/5))^2;
c4:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*4/5+2*Pi/5))^2;
eqn6:=b0^2-q^2*c0=(x1-x0)^2+(y1-y0)^2;
eqn7:=b1^2-q^2*c1=(x2-x1)^2+(y2-y1)^2;
eqn8:=b2^2-q^2*c2=(x3-x2)^2+(y3-y2)^2;
eqn9:=b3^2-q^2*c3=(x4-x3)^2+(y4-y3)^2;
eqn10:=b4^2-q^2*c4=(x0-x4)^2+(y0-y4)^2;
soln:=fsolve({eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10},
{x0=x0i,x1=x1i,x2=x2i,x3=x3i,x4=x4i,y0=y0i,y1=y1i,y2=y2i,y3=y3i,y4=y4i});
for i from 1 to 10 do
  if op(1,soln[i])=x0 then x0s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x1 then x1s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x2 then x2s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x3 then x3s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x4 then x4s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y0 then y0s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y1 then y1s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y2 then y2s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y3 then y3s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y4 then y4s:=op(2,soln[i]) end if;
end do;
z0s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*0/5));
z1s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*1/5));
z2s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*2/5));
z3s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*3/5));
z4s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*4/5));
XC0:=Vector(3,{1=x0s,2=y0s,3=z0s});
XC1:=Vector(3,{1=x1s,2=y1s,3=z1s});
XC2:=Vector(3,{1=x2s,2=y2s,3=z2s});
XC3:=Vector(3,{1=x3s,2=y3s,3=z3s});
XC4:=Vector(3,{1=x4s,2=y4s,3=z4s});
tau1:=bondangle(XC0,XC1,XC2);
tau4:=bondangle(XC3,XC4,XC0);
phi1203:=wedgeangle(XC1,XC2,XC0,XC3);
phi4302:=wedgeangle(XC4,XC3,XC0,XC2);
return tau1,tau4,phi1203,phi4302
end proc;

fiveMD2int:=proc(b0,b1,b2,b3,b4,q,P,S,Gamma)
#
local bav,x0i,x1i,x2i,x3i,x4i,y0i,y1i,y2i,y3i,y4i,
eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10,soln,
c0,c1,c2,c3,c4,c5,i,x0s,x1s,x2s,x3s,x4s,y0s,y1s,y2s,y3s,y4s,
z0s,z1s,z2s,z3s,z4s,XC0,XC1,XC2,XC3,XC4,tau0,tau1,tau2,tau3,tau4,
nu0,nu1,nu2,nu3,X1,X2,X3,rho,theta,phi,alpha,psi1,psi2,
theta0,phi0,phi01,phi02,theta1,phi11,phi12,theta2,phi21,phi22,
theta3,phi31,phi32,theta4,phi4,phi41,phi42;
#
description `This converts Marzec-Day internal coordinates`,
`for a five membered ring into ordinary internal coordinates.`, 
`b0, .., b4 are the bond lengths in angstroms`,
`q,P are the Cremer-Pople pseudorotation amplitude and phase`,
`S,Gamma are the Marzec-Day elongation amplitude and phase`,
`Besides the bond lengths, the ordinary internal coordinates are:`,
`tau1,tau2,tau3: bond angles with vertices C1, C2, C3 respectively`,
`nu0,nu1,nu2,nu3: the standard torsion angles with those names.`,
`It is assumed that the user has previously typed:`,
`with(LinearAlgebra);.`;
#
bav:=evalf((b0+b1+b2+b3+b4)/5);
x0i:=0.;
y0i:=bav;
x1i:=evalf(bav*cos(Pi/2-2*Pi*1/5));
y1i:=evalf(bav*sin(Pi/2-2*Pi*1/5));
x2i:=evalf(bav*cos(Pi/2-2*Pi*2/5));
y2i:=evalf(bav*sin(Pi/2-2*Pi*2/5));
x3i:=evalf(bav*cos(Pi/2-2*Pi*3/5));
y3i:=evalf(bav*sin(Pi/2-2*Pi*3/5));
x4i:=evalf(bav*cos(Pi/2-2*Pi*4/5));
y4i:=evalf(bav*sin(Pi/2-2*Pi*4/5));
eqn1:=x0+x1+x2+x3+x4=0;
eqn2:=y0+y1+y2+y3+y4=0;
eqn3:=x0=0;
eqn4:=S^2*cos(2*Gamma)=x0^2-y0^2+x1^2-y1^2+x2^2-y2^2+x3^2-y3^2+x4^2-y4^2;
eqn5:=S^2*sin(2*Gamma)=2*(x0*y0+x1*y1+x2*y2+x3*y3+x4*y4);
c0:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*0/5+2*Pi/5))^2;
c1:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*1/5+2*Pi/5))^2;
c2:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*2/5+2*Pi/5))^2;
c3:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*3/5+2*Pi/5))^2;
c4:=(2/5)*(2*sin(2*Pi/5))^2*(cos(P+4*Pi*4/5+2*Pi/5))^2;
eqn6:=b0^2-q^2*c0=(x1-x0)^2+(y1-y0)^2;
eqn7:=b1^2-q^2*c1=(x2-x1)^2+(y2-y1)^2;
eqn8:=b2^2-q^2*c2=(x3-x2)^2+(y3-y2)^2;
eqn9:=b3^2-q^2*c3=(x4-x3)^2+(y4-y3)^2;
eqn10:=b4^2-q^2*c4=(x0-x4)^2+(y0-y4)^2;
soln:=fsolve({eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8,eqn9,eqn10},
{x0=x0i,x1=x1i,x2=x2i,x3=x3i,x4=x4i,y0=y0i,y1=y1i,y2=y2i,y3=y3i,y4=y4i});
for i from 1 to 10 do
  if op(1,soln[i])=x0 then x0s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x1 then x1s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x2 then x2s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x3 then x3s:=op(2,soln[i]) end if;
  if op(1,soln[i])=x4 then x4s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y0 then y0s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y1 then y1s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y2 then y2s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y3 then y3s:=op(2,soln[i]) end if;
  if op(1,soln[i])=y4 then y4s:=op(2,soln[i]) end if;
end do;
z0s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*0/5));
z1s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*1/5));
z2s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*2/5));
z3s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*3/5));
z4s:=evalf(sqrt(2/5)*q*cos(P-Pi/2+4*Pi*4/5));
XC0:=Vector(3,{1=x0s,2=y0s,3=z0s});
XC1:=Vector(3,{1=x1s,2=y1s,3=z1s});
XC2:=Vector(3,{1=x2s,2=y2s,3=z2s});
XC3:=Vector(3,{1=x3s,2=y3s,3=z3s});
XC4:=Vector(3,{1=x4s,2=y4s,3=z4s});
X1:=Vector(3,{1=-5,2=-4,3=0});
X2:=Vector(3,{1=-5,2=-4,3=1});
X3:=Vector(3,{1=-4,2=-4,3=0});
tau0:=bondangle(XC4,XC0,XC1);
tau1:=bondangle(XC0,XC1,XC2);
tau2:=bondangle(XC1,XC2,XC3);
tau3:=bondangle(XC2,XC3,XC4);
tau4:=bondangle(XC3,XC4,XC0);
nu0:=wedgeangle(XC4,XC0,XC1,XC2);
nu1:=wedgeangle(XC0,XC1,XC2,XC3);
nu2:=wedgeangle(XC1,XC2,XC3,XC4);
nu3:=wedgeangle(XC2,XC3,XC4,XC0);
#
# Here we compute some exocyclic bond angles and torsion angles.
# This is why we need tau0, tau4, and nu0, nu3.
#
theta0:=exobondangle(tau0);
phi0:=exotorsion(theta0);
phi01:=nu0+phi0;
phi02:=nu0-phi0;
#
theta1:=exobondangle(tau1);
phi11:=-exotorsion(theta1);
phi12:=-phi11;
#
theta2:=exobondangle(tau2);
phi21:=-exotorsion(theta2);
phi22:=-phi21;
#
theta3:=exobondangle(tau3);
phi31:=-exotorsion(theta3);
phi32:=-phi31;
#
theta4:=exobondangle(tau4);
phi4:=exotorsion(theta4);
phi41:=nu3-phi4;
phi42:=nu3+phi4;
#
# Here we compute some polar coordinates to orient the molecule
# in the Marzec-Day standard frame.
#
rho:=Norm(XC0-X1,2);
theta:=bondangle(XC0,X1,X2);
phi:=wedgeangle(X3,X1,X2,XC0);
alpha:=bondangle(X1,XC0,XC1);
psi1:=wedgeangle(X2,X1,XC0,XC1);
psi2:=wedgeangle(X1,XC0,XC1,XC2);
#
return `radial arm rho=`,rho,`spherical angle theta=`,theta,
`orientation angle alpha=`,alpha,`exocyclic theta0=`,theta0,
`endocyclic tau1=`,tau1,`exocyclic theta1=`,theta1,`endocyclic tau2=`,tau2,
`exocyclic theta2=`,theta2,`endocyclic tau3=`,tau3,`exocyclic theta3=`,theta3,
`exocyclic theta4=`,theta4,`spherical wedge phi=`,phi,`orientation psi1=`,psi1,
`orientation psi2=`,psi2,`(exocyclic phi01=`,phi01,`exocyclic phi02=`,phi02,
`proline exocyclic wedge angle nu0-180=`,evalf(nu0-180.0),
`)(exocyclic phi11=`,phi11,`exocyclic phi12=`,phi12,`) endocyclic nu1=`,nu1,
`(exocyclic phi21=`,phi21,`exocyclic phi22=`,phi22,`) endocyclic nu2=`,nu2,
`(exocyclic phi31=`,phi31,`exocyclic phi32=`,phi32,
`) (exocyclic phi41=`,phi41,`exocyclic phi42=`,phi42,`)`
end proc;

Zsystem:=proc(b0,b1,b2,b3,b4,tau1,tau4)
#
local l03sq,l02sq,l03,l02,a203,a430,a120;
#
description `This computes the additional labels of the trees in the`,
`Z-system for a five-membered ring. We are given the five bonds lengths`,
`b0,..,b4, and the two angles tau1,tau4 (in degrees). It is necessary to compute`,
`l03, l02: distances between atom pairs {A0,A3} and {A0,A2} resp.`,
`a430,a203,a120: angles between the indicated bond pairs, where the`,
`middle index is the number of the vertex atom of the angle.`,
`Fortunately this is trivial trigonometry.`;
#
l03sq:=evalf(b3^2+b4^2-2*b3*b4*cos(tau4*Pi/180));
l03:=sqrt(l03sq);
l02sq:=evalf(b0^2+b1^2-2*b0*b1*cos(tau1*Pi/180));
l02:=sqrt(l02sq);
a203:=evalf(180/Pi*arccos((l02sq+l03sq-b2^2)/(2*l02*l03)));
a430:=evalf(180/Pi*arccos((b3^2+l03sq-b4^2)/(2*b3*l03)));
a120:=evalf(180/Pi*arccos((b1^2+l02sq-b0^2)/(2*b1*l02)));
return l02,l03,a120,a203,a430
end proc;

exobondangle:=proc(tau)
#
local theta;
#
description `Given an endocyclic bond angle tau (in degrees) we compute`,
`the common value of the exocyclic bond angles assuming the vertex of the`,
`endocyclic bond angle has four substituents.`;
#
theta:=evalf(180/Pi*arccos(-2/(1+sqrt(1+16/(1+cos(tau/180*Pi))))))
end proc;

exotorsion:=proc(theta)
#
local phi;
#
description `Given the exoyclic bond angle theta (in degrees) we compute`,
`the exocyclic torsion angle as measured from the plane of the endocyclic`,
`bond angle.`;
#
phi:=evalf(180/Pi*(Pi-1/2*arccos(cos(theta/180*Pi)/(1+cos(theta/180*Pi)))))
end proc;