(* ::Package:: *)

ODtrans::usage = "ODtrans attaches lone pairs to O in C=O of arbitrary molecules.
Input: t, oj, conn with target 2, Okk, Okl, rlo;
Output: tplo,tplu."


k=1;
Plota= Table[{0}, {i, 1, oj}]; 
Plotc= Plota;
Plotd =Plota; 
Plote= Plota;
Plotb0 =Plota; 


Do[(c = conn[[k]]; 
a1 = dc[[c[[2]],c[[1]]]]; a2 = dc[[c[[2]],c[[3]]]]; d1 = dc[[c[[1]],c[[3]]]]; 
\[Alpha] = ArcCos[(a1^2 + a2^2 - d1^2)/(2*a1*a2)]; 
L = u2[[c[[2]]]] + rlo; 
\[Beta] = ArcSin[rlo/L]; 
\[Gamma] = 3*\[Pi]/2 - \[Alpha]/2; 
v = Table[Join[t[[c[[i]]]], {1}], {i, 1, P}]; 
x2 = t[[c[[2]],1]]; 
y2 = t[[c[[2]],2]]; 
z2 = t[[c[[2]],3]]; 
Q = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {-x2, -y2, -z2, 1}}; 
p = v.Q; 
p = Table[Take[p[[i]], 3], {i, 1, P}]; 
T1 = ArcTan[p[[1,3]]/p[[1,2]]]; 
U = p; 
T = RotationMatrix[T1, {1, 0, 0}]; 
p = U.T; 
T2 = ArcTan[p[[1,2]]/p[[1,1]]]; 
U = p; 
T = RotationMatrix[T2, {0, 0, 1}]; 
p = U.T; 
T3 = ArcTan[p[[3,3]]/p[[3,2]]]; 
U = p; 
T = RotationMatrix[T3, {1, 0, 0}]; 
p = Chop[U.T]; 
nk = p; 
nl = p; 
nk[[2,1]] = L*Sin[\[Gamma] + 1.0472]; nk[[2,2]] = (-L)*Cos[\[Gamma] + 1.0472]; nk[[2,3]] = 0; 
nl[[2,1]] = L*Sin[\[Gamma] - 1.0472]; nl[[2,2]] = (-L)*Cos[\[Gamma] - 1.0472]; nl[[2,3]] = 0;
nk[[2]]=Okk[[k]]*nk[[2]]; nl[[2]]=Okl[[k]]*nl[[2]]; 
Plota[[k]] = Table[Graphics3D[{Opacity[0.5], Sphere[p[[i]], u2[[c[[i]]]]]}], {i, 1, P}]; 
Plotc [[k]] = Table[Graphics3D[{Tube[{p[[i]], p[[j]]}, 0.02]}], {i, 1, P - 1}, {j, i + 1, P}]; 
Plotd[[k]]  = Graphics3D[{Opacity[0.5], Lighter[Red], Sphere[nk[[2]], rlo]}]; 
Plote[[k]]  = Graphics3D[{Opacity[0.5], Lighter[Red], Sphere[nl[[2]], rlo]}]; 
Plotb0[[k]]  = Table[Graphics3D[Text[c[[i]], p[[i]], {1, 0}]], {i, 1, P}];
U = p; 
U1 = nk; 
U2 = nl; 
T = RotationMatrix[-T3, {1, 0, 0}]; 
p = U.T; 
nk = U1.T; 
nl = U2.T; 
U = p; 
U1 = nk; 
U2 = nl; 
T = RotationMatrix[-T2, {0, 0, 1}]; 
p = U.T; 
nk = U1.T; 
nl = U2.T; 
U = p; 
U1 = nk; 
U2 = nl; 
T = RotationMatrix[-T1, {1, 0, 0}]; 
p = U.T; 
nk = U1.T; 
nl = U2.T; 
Clear[v, v1, v2]; 
v = Table[Join[p[[i]], {1}], {i, 1, P}]; 
v1 = Table[Join[nk[[i]], {1}], {i, 1, P}]; 
v2 = Table[Join[nl[[i]], {1}], {i, 1, P}]; 
Q = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {x2, y2, z2, 1}}; 
p = v.Q; 
nk = v1.Q; 
nl = v2.Q; 
p = Table[Take[p[[i]], 3], {i, 1, P}]; 
nk = Table[Take[nk[[i]], 3], {i, 1, P}]; 
nl = Table[Take[nl[[i]], 3], {i, 1, P}]; 
tplo[[k]]= nk[[2]]; 
tplu[[k]]= nl[[2]]),{k,oj}]; 


Table[Show[Plotd[[i]],Plote[[i]], Plota[[i]], Plotb0[[i]], Plotc[[i]], {Axes -> Automatic, AxesLabel -> {xl, yl, zl}, 
        AspectRatio -> Automatic, Boxed -> False, PlotRange -> Automatic, SphericalRegion -> True}],{i,1,oj}]