// hcyl1_curv.java  compile    javac -cp . hcyl1_curv.java
//                  execute    java  -cp . hcyl1_curv
//                             writes      hcyl1_curv.dat
//                  execute    dat_to_stl hcyl1_curv.dat hcyl1_curv.stl
//
// curved version of cyl1_pipe  update phi,y,z per rc, then rotate phi about x
// make two, glue together, have pipe of a toroide.

public class hcyl1_curv  // make hcyl1_curv.dat then .stl for 3D printer
{
  double pi = 3.14159265358979323846;
  double Mx[][] = new double[3][3]; // rotation matrix
  double My[][] = new double[3][3];
  double Mz[][] = new double[3][3];

  public hcyl1_curv()
  {
    int num_points = 0;
    double data_points[][] = new double[20000][3];
    int num_polys = 0;
    int k = 0; // k index+1/4 = num_polys
    int data_polys[] = new int[20000]; // 3 pt1 pt2 pt3  +1 index
    double theta, x, y, z;
    double phi, xp, yp, zp, del; // rotated points
    double r = 10.0;
    double rc = 30.0; // curv cyl
    int pts = 17; // now 17  theta
    int tpts = 32; // now 32  phi
    double dtheta = pi/16.0; // now pi/16.0
    double dphi   = pi/16.0; // now pi/16.0;
    double vp[] = new double[3];
    double vr[] = new double[3];
    
    int pti = 0;
    int pt0=0;          // pt0  for triangles
    int pt1=0;          // next neighbor
    int pt2=0;          // for bottom
    int pt3=0;
    int ptin=tpts*pts;  // inside  pyin bottom
    
    System.out.println("hcyl1_curv.java writing hcyl1_curv.stl");
    System.out.println("r="+r+", pts="+pts+", dtheta="+dtheta);
    System.out.println("rc="+rc+", tpts="+tpts+", dphi="+dphi);

    r = 10.0;
    for(int j=0; j<tpts; j++) // curv
    {
      phi = (double)j*dphi;
      System.out.println("phi="+phi+" radians "+degrees(phi)+" degrees");
      for(int i=0; i<pts; i++) // generate vertices outside 0-4
      {
        theta = (double)i * dtheta;
        z = r * Math.sin(theta);
        x = r * Math.cos(theta);
        y = 0.0;
        makeMz(phi);
        vp[0] = x;
	vp[1] = y;
	vp[2] = z;
        matvec(Mz, vp, vr);
	xp = vr[0];
	yp = vr[1];
	zp = vr[2];
        yp = yp + rc*Math.sin(phi);
        xp = xp + rc*Math.cos(phi); // offset center on curv
        data_points[i+pti][0] = xp;
        data_points[i+pti][1] = yp;
        data_points[i+pti][2] = zp;
        num_points++;
        System.out.println(" ");
      } // end outside 0-4
      pti = pti + pts;
    } // end j, end generate outside points

    ptin = pti; // check
    System.out.println("start of inside ptin="+ptin);
    System.out.println("added points, num_points="+num_points);

    System.out.println(" ");
    System.out.println("generate triangles outside top");

    for(int j=0; j<tpts-1; j++) // curv
    {
      pt0 = j*pts;
      pt1 = (j+1)*pts;	
      for(int i=0; i<pts; i++) // generate triangles outside 1-tpts-1
      {
	if(i<pts-1)
	{
          data_polys[k++] = 3;
          data_polys[k++] = i+pt0+1;
          data_polys[k++] = i+pt1+2;
          data_polys[k++] = i+pt1+1;
          data_polys[k++] = 3;
          data_polys[k++] = i+pt0+1;
          data_polys[k++] = i+pt0+2;
          data_polys[k++] = i+pt1+2;
	} // no closing
      }
    } // end j
    // do closing
    pt0 = (tpts-1)*pts;
    pt1 = 0;	
    System.out.println("generate triangles outside cls pt0="+pt0+", pt1="+pt1);
    for(int i=0; i<pts-1; i++) // generate triangles outside tpts-1
    {
      data_polys[k++] = 3;
      data_polys[k++] = i+pt0+1;
      data_polys[k++] = i+pt1+2;
      data_polys[k++] = i+pt1+1;
      data_polys[k++] = 3;
      data_polys[k++] = i+pt0+1;
      data_polys[k++] = i+pt0+2;
      data_polys[k++] = i+pt1+2;
      num_polys = k/4;
    }
    
    System.out.println("added triangles, k="+k+", num_polys="+num_polys);

    ptin = pti;
    r = 8.0;
    for(int j=0; j<tpts; j++) // curv
    {
      phi = (double)j*dphi;
      for(int i=0; i<pts; i++) // generate vertices inside ptin 
      {
        theta = (double)i * dtheta;
        z = r * Math.sin(theta);
        x = r * Math.cos(theta);
	y = 0.0;
        System.out.println("phi="+phi+" radians "+degrees(phi)+" degrees");
        makeMz(phi);
        vp[0] = x;
	vp[1] = y;
	vp[2] = z;
        matvec(Mz, vp, vr);
	xp = vr[0];
	yp = vr[1];
	zp = vr[2];
        yp = yp + rc*Math.sin(phi);
        xp = xp + rc*Math.cos(phi); // offset center on curv
        data_points[i+pti][0] = xp;
        data_points[i+pti][1] = yp;
        data_points[i+pti][2] = zp;
        num_points++;
	System.out.println("inside phi="+phi+", xp="+xp+",yp="+yp+",zp="+zp);
        System.out.println(" ");
      }
      pti = pti + pts;
    } // end generate inside points

    System.out.println("added points, num_points="+num_points);

    System.out.println(" ");
    System.out.println("generate triangles inside");

    for(int j=0; j<tpts-1; j++) // curv
    {
      pt0 = ptin+j*pts;
      pt1 = ptin+(j+1)*pts;
      System.out.println("generate triangles inside pt0="+pt0+", pt1="+pt1);
      for(int i=0; i<pts-1; i++) // generate triangles inside 1-tpts
      {
	if(i<pts-1)
	{
          data_polys[k++] = 3;
          data_polys[k++] = i+pt0+1;
          data_polys[k++] = i+pt1+1;
          data_polys[k++] = i+pt1+2;
          data_polys[k++] = 3;
          data_polys[k++] = i+pt0+1;
          data_polys[k++] = i+pt1+2;
          data_polys[k++] = i+pt0+2;
	}
      } // end if  no closing
    } // end j

    // do closing
    pt0 = ptin +(tpts-1)*pts;
    pt1 = ptin;	
    System.out.println("generate triangles inside cls pt0="+pt0+", pt1="+pt1);
    for(int i=0; i<pts-1; i++) // generate triangles inside tpts-1
    {
      data_polys[k++] = 3;
      data_polys[k++] = i+pt0+1;
      data_polys[k++] = i+pt1+1;
      data_polys[k++] = i+pt1+2;
      data_polys[k++] = 3;
      data_polys[k++] = i+pt0+1;
      data_polys[k++] = i+pt1+2;
      data_polys[k++] = i+pt0+2;
      num_polys = k/4;
    }

    pt0 = 0;
    pt1 = (tpts)*pts;
    System.out.println("generate triangles out bottom pt0="+pt0+", pt1="+pt1);
    for(int i=0; i<tpts-1; i++) // generate triangles outside bottom
    {
      data_polys[k++] = 3;
      data_polys[k++] = pt0+1;
      data_polys[k++] = pt1+1+pts;
      data_polys[k++] = pt1+1;
      data_polys[k++] = 3;
      data_polys[k++] = pt0+1;
      data_polys[k++] = pt0+1+pts;
      data_polys[k++] = pt1+1+pts;
      pt0 = pt0+pts;
      pt1 = pt1+pts;
    }
    pt2 = 0;
    pt3 = (tpts)*pts;
    data_polys[k++] = 3;
    data_polys[k++] = pt0+1;
    data_polys[k++] = pt3+1;
    data_polys[k++] = pt1+1;
    data_polys[k++] = 3;
    data_polys[k++] = pt0+1;
    data_polys[k++] = pt2+1;
    data_polys[k++] = pt3+1;
    num_polys = k/4;
    
    System.out.println("added triangles, k="+k+", num_polys="+num_polys);

    System.out.println("generate triangles insi bottom pt0="+pt0+", pt1="+pt1);
    pt1 = pts-1;
    pt0 = (tpts)*pts+pts-1;
    for(int i=0; i<tpts-1; i++) // generate triangles inside bottom
    {
      data_polys[k++] = 3;
      data_polys[k++] = pt0+1;
      data_polys[k++] = pt1+1+pts;
      data_polys[k++] = pt1+1;
      data_polys[k++] = 3;
      data_polys[k++] = pt0+1;
      data_polys[k++] = pt0+1+pts;
      data_polys[k++] = pt1+1+pts;
      pt0 = pt0+pts;
      pt1 = pt1+pts;
    }
    pt3 = pts-1;
    pt2 = (tpts)*pts+pts-1;
    data_polys[k++] = 3;
    data_polys[k++] = pt0+1;
    data_polys[k++] = pt3+1;
    data_polys[k++] = pt1+1;
    data_polys[k++] = 3;
    data_polys[k++] = pt0+1;
    data_polys[k++] = pt2+1;
    data_polys[k++] = pt3+1;
    num_polys = k/4;
    
    System.out.println("added triangles, k="+k+", num_polys="+num_polys);



    System.out.println("generate triangles bottom pt0="+pt0+", pt1="+pt1);

    System.out.println("added triangles, k="+k+", num_polys="+num_polys);
    System.out.println(" ");
  
    new datwrite("hcyl1_curv.dat", data_points, num_points, data_polys, num_polys);
    
  } // end class hcyl1_curv

  void makeMx(double phi)
  {
    identm(Mx);
    Mx[1][1] = Math.cos(phi);
    Mx[2][2] = Math.cos(phi);
    Mx[2][1] = Math.sin(phi);
    Mx[1][2] = - Math.sin(phi);
    printm(Mx, "rotate about X Mx=");
  } // end makeMx

  void makeMy(double phi)
  {
    identm(My);
    My[0][0] = Math.cos(phi);
    My[2][2] = Math.cos(phi);
    My[0][2] = Math.sin(phi);
    My[2][0] = -Math.sin(phi);
    printm(My, "rotate about Y My=");
  }

  void makeMz(double phi)
  {
    identm(Mz);
    Mz[0][0] = Math.cos(phi);
    Mz[1][1] = Math.cos(phi);
    Mz[1][0] = Math.sin(phi);
    Mz[0][1] = -Math.sin(phi);
    printm(Mz, "rotate about Z Mz=");
  }

  void zerom(double A[][])
  {
    int i, j;
    for(i=0; i<3; i++)
      for(j=0; j<3; j++)
        A[i][j] = 0.0;
  } // end zerom

  void identm(double A[][])
  {
    int i, j;
    for(i=0; i<3; i++)
    {
      for(j=0; j<3; j++)
        A[i][j] = 0.0;
      A[i][i] = 1.0;
    }
  } // end identm

  void printm(double A[][], String txt)
  {
    int i, j;
    System.out.println(" ");
    System.out.println(txt);
    for(i=0; i<3; i++)
    {
      System.out.print("| ");
      for(j=0; j<3; j++)
	  System.out.print(A[i][j]+"   ");
      System.out.println(" |");
    }
  } // end printm

  void matmul(double A[][], double B[][], double C[][])
  {
    int i, j, k;
    for(i=0; i<3; i++)
    {
      for(j=0; j<3; j++)
      {
        C[i][j] = 0.0;
        for(k=0; k<3; k++)
          C[i][j] += A[i][k]*B[k][j];
      }
    }
  } // end matmul

  void matvec(double A[][], double V[], double W[])
  {
    int i, j;
    for(i=0; i<3; i++)
    {
      W[i] = 0.0;
      for(j=0; j<3; j++)
        W[i] += A[i][j]*V[j];
    } 
  } // end matvec

  double degrees(double a)
  {
    return 180.0*a/pi;
  }

  public static void main (String[] args)
  {
    new hcyl1_curv();
  }
} // end  hcyl1_curv.java