//---------------------------------------------------------------------------
//--- tube with sinus screw -------------------------------------------------
//---------------------------------------------------------------------------
// https://stackoverflow.com/a/54050883/2521214
//---------------------------------------------------------------------------
const int   ca= 20;     // points per slice
const int   cb=100;     // slices
const float r0=  0.3;   // minor screw radius
const float r1=  0.35;  // major screw radius
const float l1=  4.0;   // tube length
const float nz= 10.0;   // screws
//---------------------------------------------------------------------------
vec3 bend_pc; int bend_j0,bend_j1;  // just for debug draw
//---------------------------------------------------------------------------
// straight mesh
vec3 pnt0[ca][cb];  // vertex
vec3 nor0[ca][cb];  // normal
vec2 txr0[ca][cb];  // texcoord
vec3 mid0[cb];      // slice center
vec3 dir0[cb];      // slice central axis (normalized)
float len0[cb];     // slice arclength position
//---------------------------------------------------------------------------
// bended mesh
vec3 pnt1[ca][cb];  // vertex
vec3 nor1[ca][cb];  // normal
vec2 txr1[ca][cb];  // texcoord
vec3 mid1[cb];      // slice center
vec3 dir1[cb];      // slice central axis (normalized)
float len1[cb];     // slice arclength position
//---------------------------------------------------------------------------
void obj0_init()    // sin screw
    {
    int     i,j,i0,j0;
    float   a,b,l,da,db,dl,r,s,c,tx,ty;
    float   dtx=1.0/float(ca-1),dty=1.0/float(cb-1);
    vec3    u,v;
    // pnt,txr
    da=2.0*M_PI/float(ca-1);
    db=nz*2.0*M_PI/float(cb);
    dl=l1/float(cb);
    for (a=0.0,tx=0.0,i=0;i<ca;i++,a+=da,tx+=dtx)
        {
        s=sin(a);
        c=cos(a);
        for (l=-0.5*l1,b=0,ty=0.0,j=0;j<cb;j++,b+=db,l+=dl,ty+=dty)
            {
            r=r0+((r1-r0)*cos(a+b));
            pnt0[i][j].x=r*c;
            pnt0[i][j].y=r*s;
            pnt0[i][j].z=l;
            txr0[i][j].x=tx;
            txr0[i][j].y=ty;
            }
        }
    // mid,dir
    for (l=0.0,j=0;j<cb;j++,l+=dl)
        {
        mid0[j]=vec3(0.0,0.0, l-(0.5*l1));
        dir0[j]=vec3(0.0,0.0,dl);
        len0[j]=l;
        }
    // nor
    for (i0=ca-2,i=0;i<ca;i0=i,i++)
     for (j0=cb-1,j=0;j<cb;j0=j,j++)
        {
        u=pnt0[i][j]-pnt0[i0][j];
        v=pnt0[i][j]-pnt0[i][j0];
        nor0[i][j]=normalize(cross(u,v));
        }
    }
//---------------------------------------------------------------------------
void obj1_copy()    // obj1 = obj0
    {
    int i,j;
    for (i=0;i<ca;i++)
     for (j=0;j<cb;j++)
        {
        pnt1[i][j]=pnt0[i][j];
        txr1[i][j]=txr0[i][j];
        nor1[i][j]=nor0[i][j];
        }
    for (j=0;j<cb;j++)
        {
        mid1[j]=mid0[j];
        dir1[j]=dir0[j];
        len1[j]=len0[j];
        }
    }
//---------------------------------------------------------------------------
vec3 rotatex(vec3 p,vec3 p0,float a)
    {
    vec3 q; p-=p0;
    q.z=+(p.z*cos(a))+(p.y*sin(a));
    q.y=-(p.z*sin(a))+(p.y*cos(a));
    q.x=p.x;
    return q+p0;
    }
//---------------------------------------------------------------------------
vec3 rotatey(vec3 p,vec3 p0,float a)
    {
    vec3 q; p-=p0;
    q.x=+(p.x*cos(a))+(p.z*sin(a));
    q.z=-(p.x*sin(a))+(p.z*cos(a));
    q.y=p.y;
    return q+p0;
    }
//---------------------------------------------------------------------------
vec3 rotatez(vec3 p,vec3 p0,float a)
    {
    vec3 q; p-=p0;
    q.x=+(p.x*cos(a))+(p.y*sin(a));
    q.y=-(p.x*sin(a))+(p.y*cos(a));
    q.z=p.z;
    return q+p0;
    }
//---------------------------------------------------------------------------
void obj1_bendx(float l0,float l1,float ang)    // [units],[units],[rad] bend obj1 around x axis
    {
    int i,j,i0,j0,j1;
    float a,r,l;
    vec3 PC,p,u,v;
    vec3 P0,X0,Y0,Z0;
    // find start and end of bend
    for (j0= 0;(j0<cb)&&(len1[j0]<l0);j0++);
    for (j1=j0;(j1<cb)&&(len1[j1]<l1);j1++);
    if (j0>cb) return;          // no bend
    // coordinate system0
    P0=mid1[j0];
    Z0=normalize(dir1[j0]);
    X0=vec3(1.0,0.0,0.0);
    Y0=cross(Z0,X0);
    X0=cross(Y0,Z0);
    // bend center
    r=(l1-l0)/ang;
    PC=P0-(Y0*r);
    r=fabs(r);
    // just for debug draw
    bend_pc=PC;
    bend_j0=j0;
    bend_j1=j1;
    // bend <l0,l1)
    for (j=j0;j<cb;j++)
        {
        // arc length -> angle [rad] and length correction
        if (j<j1)
            {
            a=ang*(len1[j]-l0)/(l1-l0);
            p=Z0*(len1[j]-l0);
            }
        else{
            a=ang;
            p=Z0*(l1-l0);
            }
        // transform system0 -> system1
        mid1[j]=rotatex(mid1[j]-p,PC,a);
        dir1[j]=rotatex(dir1[j],vec3(0.0,0.0,0.0),a);
        for (i=0;i<ca;i++) pnt1[i][j]=rotatex(pnt1[i][j]-p,PC,a);
        }
    // nor
    for (i0=ca-2,i=0;i<ca;i0=i,i++)
     for (j0=cb-1,j=0;j<cb;j0=j,j++)
        {
        u=pnt1[i][j]-pnt1[i0][j];
        v=pnt1[i][j]-pnt1[i][j0];
        nor1[i][j]=normalize(cross(u,v));
        }
    }
//---------------------------------------------------------------------------
void obj0_draw()
    {
    int i,j;
    glColor3f(1.0,1.0,1.0);
    glEnable(GL_CULL_FACE);
    glFrontFace(GL_CW);
    for (i=0;i<ca-1;i++)
        {
        glBegin(GL_QUAD_STRIP);
        for (j=0;j<cb;j++)
            {
            glTexCoord2fv(txr0[i+1][j].dat);
            glNormal3fv  (nor0[i+1][j].dat);
            glVertex3fv  (pnt0[i+1][j].dat);
            glTexCoord2fv(txr0[i  ][j].dat);
            glNormal3fv  (nor0[i  ][j].dat);
            glVertex3fv  (pnt0[i  ][j].dat);
            }
        glEnd();
        }
    }
//---------------------------------------------------------------------------
void obj1_draw()
    {
    int i,j;
    glColor3f(1.0,1.0,1.0);
    glEnable(GL_CULL_FACE);
    glFrontFace(GL_CW);
    for (i=0;i<ca-1;i++)
        {
        glBegin(GL_QUAD_STRIP);
        for (j=0;j<cb;j++)
            {
            glTexCoord2fv(txr1[i+1][j].dat);
            glNormal3fv  (nor1[i+1][j].dat);
            glVertex3fv  (pnt1[i+1][j].dat);
            glTexCoord2fv(txr1[i  ][j].dat);
            glNormal3fv  (nor1[i  ][j].dat);
            glVertex3fv  (pnt1[i  ][j].dat);
            }
        glEnd();
        }
    }
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------