--- a
+++ b/Tools/ModelUtilities/WrappingSurfaces/WrappingCylinder5PointFit.any
@@ -0,0 +1,60 @@
+ // required input:
+// Z_rotation
+//  AnyDrawNode drw = {};
+//AnyDrawRefFrame drw ={ScaleXYZ={1,1,1}*0.2;};
+
+//
+//AnyDrawSphere p1 = {Position = (.p[0]-.sRel)*.ARel;ScaleXYZ = 0.003*{1,1,1};};
+//AnyDrawSphere p2 = {Position = (.p[1]-.sRel)*.ARel;ScaleXYZ = 0.003*{1,1,1};};
+//AnyDrawSphere p3 = {Position = (.p[2]-.sRel)*.ARel;ScaleXYZ = 0.003*{1,1,1};};
+//AnyDrawSphere p4 = {Position = (.p[3]-.sRel)*.ARel;ScaleXYZ = 0.003*{1,1,1};};
+//AnyDrawSphere p5 = {Position = (.p[4]-.sRel)*.ARel;ScaleXYZ = 0.003*{1,1,1};};
+//
+//AnyDrawRefFrame drw = {ScaleXYZ = 0.1*{1,1,1}; V};
+
+
+// transform landmarks by applying transformation and scaling
+AnyFloat Q1 = p[0];
+AnyFloat Q2 = p[1];
+AnyFloat Q3 = p[2];
+AnyFloat Q4 = p[3];
+AnyFloat Q5 = p[4];
+
+AnyFloat Length = lengthscalefactor * vnorm(Q2 - Q1);
+
+// get direction from the first two points
+AnyFloat dir = (Q2 - Q1) / vnorm(Q2 - Q1);
+
+// projection of the points 3-5 to the plane through Q1 with normal dir
+AnyFloat Q3p = Q1 + (Q3 - Q1) - dir * (Q3 - Q1)' * dir;
+AnyFloat Q4p = Q1 + (Q4 - Q1) - dir * (Q4 - Q1)' * dir;
+AnyFloat Q5p = Q1 + (Q5 - Q1) - dir * (Q5 - Q1)' * dir;
+
+AnyFloat alpha = ((Q4p - Q5p) * (Q4p - Q5p)' * (Q3p - Q4p) * (Q3p - Q5p)') 
+/ (2 * cross(Q3p - Q4p, Q4p - Q5p) * cross(Q3p - Q4p, Q4p - Q5p)');
+AnyFloat beta = ((Q3p - Q5p) * (Q3p - Q5p)' * (Q4p - Q3p) * (Q4p - Q5p)') 
+/ (2 * cross(Q3p - Q4p, Q4p - Q5p) * cross(Q3p - Q4p, Q4p - Q5p)');
+AnyFloat gamma = ((Q3p - Q4p) * (Q3p - Q4p)' * (Q5p - Q3p) * (Q5p - Q4p)') 
+/ (2 * cross(Q3p - Q4p, Q4p - Q5p) * cross(Q3p - Q4p, Q4p - Q5p)');
+AnyVec3 U = (alpha * Q3p + beta * Q4p + gamma * Q5p);
+
+AnyFloat Radius = vnorm(U-Q3p);
+AnyVec3 dir1 = (Q3p - U) / vnorm(Q3p - U);
+AnyVec3 dir2tmp = (Q4p - U) - (Q4p - U) * dir1' * dir1;
+AnyVec3 dir2 = dir2tmp / vnorm(dir2tmp);
+
+AnyFloat  signtmp = cross(dir1,dir2)*dir';
+AnyFloat sign = iffun(gtfun(signtmp,0),1,-1);
+
+ARel = RotMat(U, U + dir1, U  +sign*dir2)*RotMat(Z_rotation,z);
+
+sRel = U- (Length-vnorm(Q2-Q1))/2 *dir;
+//AnyDrawRefFrame drw ={RGB={0,0,1};ScaleXYZ={1,1,1}*0.1;};
+
+
+
+AnySurfCylinder cyl = {
+  Radius=.Radius*.radiusscalefactor;
+  Length = abs(.Length);
+//  AnyDrawParamSurf drw={RGB={1,0,0};Opacity = 0.05;Visible = On;};
+};