// inputs
length = 1;
pitch = 1;
alpha_i = Pi/3; // intermesh angle [0,Pi/2]
alpha_f = alpha_i; // flight angle
alpha_t = Pi/2 - alpha_f; // tip angle
alpha_r = alpha_t; // root angle

// based on unit screw diameter
Cd = Cos(alpha_i/2);  // centerline distance
Dr = 2*Cd-1; // root diameter
Hc = 1-Cd; // channel depth
Ih = Sin(alpha_i/2)/2; // intermesh height

lc = 0.1;
offset_x = Cd/2; // from x=0
offset_y = 0; // from y=0

// first screw

c1p0 = newp; // center point
Point(c1p0) = {-offset_x+0, offset_y+0, 0, lc};

angle = alpha_t/2;
xt1 = 1/2*Cos(angle);
yt1 = 1/2*Sin(angle);

angle = angle + alpha_f;
xr1 = Dr/2*Cos(angle);
yr1 = Dr/2*Sin(angle);

angle = angle + alpha_r;
xr2 = Dr/2*Cos(angle);
yr2 = Dr/2*Sin(angle);

angle = angle + alpha_f;
xt2 = 1/2*Cos(angle);
yt2 = 1/2*Sin(angle);

angle = angle + alpha_t;
xt3 = 1/2*Cos(angle);
yt3 = 1/2*Sin(angle);

angle = angle + alpha_f;
xr3 = Dr/2*Cos(angle);
yr3 = Dr/2*Sin(angle);

angle = angle + alpha_r;
xr4 = Dr/2*Cos(angle);
yr4 = Dr/2*Sin(angle);

angle = angle + alpha_f;
xt4 = 1/2*Cos(angle);
yt4 = 1/2*Sin(angle);

c1p1 = newp; Point(c1p1) = {-offset_x+xt1, offset_y+yt1, 0, lc};
c1p2 = newp; Point(c1p2) = {-offset_x+xr1, offset_y+yr1, 0, lc};
c1p3 = newp; Point(c1p3) = {-offset_x+xr2, offset_y+yr2, 0, lc};
c1p4 = newp; Point(c1p4) = {-offset_x+xt2, offset_y+yt2, 0, lc};
c1p5 = newp; Point(c1p5) = {-offset_x+xt3, offset_y+yt3, 0, lc};
c1p6 = newp; Point(c1p6) = {-offset_x+xr3, offset_y+yr3, 0, lc};
c1p7 = newp; Point(c1p7) = {-offset_x+xr4, offset_y+yr4, 0, lc};
c1p8 = newp; Point(c1p8) = {-offset_x+xt4, offset_y+yt4, 0, lc};

c1a = newc; Circle(c1a) = {c1p2,c1p0,c1p3};
c1b = newc; Circle(c1b) = {c1p4,c1p0,c1p5};
c1c = newc; Circle(c1c) = {c1p6,c1p0,c1p7};
c1d = newc; Circle(c1d) = {c1p8,c1p0,c1p1};

// flank areas
alpha = -1/2*((xt1-xr1)*(xt1-xr1)+(yt1-yr1)*(yt1-yr1))/(yt1-yr1);
beta = -(xt1-xr1)/(yt1-yr1);
x0 = xr1 + alpha*beta/(1+beta*beta)*(1-Sqrt(1-(1+1/(beta*beta))*(1-Cd*Cd/(alpha*alpha))));
y0 = yr1 - alpha - beta*(xr1-x0);

p = newp; Point(p) = {-offset_x+x0, offset_y+y0, 0, lc};
c1e = newc; Circle(c1e) = {c1p1,p,c1p2};

p = newp; Point(p) = {-offset_x-x0, offset_y+y0, 0, lc};
c1f = newc; Circle(c1f) = {c1p3,p,c1p4};

p = newp; Point(p) = {-offset_x-x0, offset_y-y0, 0, lc};
c1g = newc; Circle(c1g) = {c1p5,p,c1p6};

p = newp; Point(p) = {-offset_x+x0, offset_y-y0, 0, lc};
c1h = newc; Circle(c1h) = {c1p7,p,c1p8};

ll = newl; Line Loop(ll) = {c1f, c1b, c1g, c1c, c1h, c1d, c1e, c1a};
s = newl; Plane Surface(s) = {ll};

Extrude {{0,0,length}, {0,0,1}, {-offset_x, -offset_y, 0}, 2*Pi*pitch*length} {
  Surface{s};
}