(* ::Package:: *)

OD1trans::usage = "OD1trans attaches lone pairs to O in C=O of arbitrary molecules, 
given three atoms with the target as second in the list coco. The triangle defines
the plane in which the LP's sit. Since the routine sometimes does not choose the
correct sign of the x or y coordinates, lists Okk and Okl must be given as input
with {{1,1,1},{-1,1,1}...} for either lobe accepting or changing a sign.
Input: t, oj, coco 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 = coco[[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}]