(* ::Package:: *)

LpyrNtrans::usage = "LpyrNtrans takes as input connectivity quadruples with a 
pyramidal N or P from arbitrary molecular structures. The pyramidal atom is 
element 4 in the conn list. It translates the base of the pyramid into the xy 
plane and then attaches a lone pair to the apex after which it backtranslates 
the triangle and the apex with LP back into the molecular coordinates. 
On exit the LP ccordinates are in the list tpyrn."


(*k = 1;
Clear[P];
P = 4;*)


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]]]] + rln);
   \[Beta] = ArcSin[rln/L];
   \[Gamma] = \[Pi]/2 - \[Alpha]/2; (*\[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;
   nk[[4,3]] = nk[[4,3]] + f[[k]]*(u2[[c[[4]]]] + rln); 
   U = p;
   U1 = nk;
   T = RotationMatrix[-T3, {1, 0, 0}];
   p = U.T;
   nk = U1.T;
   U = p;
   U1 = nk;
   T = RotationMatrix[-T2, {0, 0, 1}];
   p = U.T;
   nk = U1.T;
   U = p;
   U1 = nk;
   T = RotationMatrix[-T1, {1, 0, 0}];
   p = U.T;
   nk = U1.T;
   Clear[v, v1, v2];
   v  = Table[Join[p[[i]], {1}], {i, 1, P}];
   v1 = Table[Join[nk[[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;
   nk1 = v1.Q;
   p = Table[Take[p[[i]], 3], {i, 1, P}];
    nk = Table[Take[nk1[[i]], 3], {i, 1, P}];
   tpyrn[[k]] = nk[[4]]), {k, pn}];