(* ::Package:: *)

PItrans::usage = "PItrans is similar to PIflattrans. It does not transform the whole
molecular skeleton into a plane and is thus suited for molecules with double bonds
in arbitrary orientations: C=C, C=O, and others e.g. C=N.
Input: Table t of the atomic coordinates, connection table connc of double bonded
atoms with the target atom at first position, cj its length, rpo, rpi the radii of 
the \[Pi] clouds above and below their plane. If different atoms, then rpo={a,b,c,...}, 
rpi={a,b,c,...}.  The tables tpo, tpi contain the coordinates of the \[Pi] clouds on exit."


P=3; c={0,0,0,0,0}; k=1;


Do[(c = connc[[k]]; (* Target 1 *) 
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 = U.T; 
mk = p; 
ml = p; 
mk[[1,3]] = rpo[[k]]; 
ml[[1,3]] = -rpi[[k]]; 
U = p; 
U1 = mk; 
U2 = ml; 
T = RotationMatrix[-T3, {1, 0, 0}]; 
p = U.T; 
mk = U1.T; 
ml = U2.T; 
U = p; 
U1 = mk; 
U2 = ml; 
T = RotationMatrix[-T2, {0, 0, 1}]; 
p = U.T; 
mk = U1.T; 
ml = U2.T; 
U = p; 
U1 = mk; 
U2 = ml; 
T = RotationMatrix[-T1, {1, 0, 0}]; 
p = U.T; 
mk = U1.T; 
ml = U2.T; 
Clear[v, v1, v2]; 
v = Table[Join[p[[i]], {1}], {i, 1, P}];
v1 = Table[Join[mk[[i]], {1}], {i, 1, P}]; 
v2 = Table[Join[ml[[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; 
mk = v1.Q; 
ml = v2.Q; 
p = Table[Take[p[[i]], 3], {i, 1, P}]; 
mk = Table[Take[mk[[i]], 3], {i, 1, P}]; 
ml = Table[Take[ml[[i]], 3], {i, 1, P}];
tpo[[k]] = mk[[1]]; 
tpu[[k]] = ml[[1]]), {k, cj}];