// hcyl1_curv16.java  compile    javac -cp . hcyl1_curv16.java
//                  execute    java  -cp . hcyl1_curv16
//                             writes      hcyl1_curv16.dat
//                  execute    dat_to_stl hcyl1_curv16.dat hcyl1_curv16.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_curv16  // make hcyl1_curv16.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_curv16()
  {
    boolean debug = false;
    int num_points = 0;
    double data_points[][] = new double[5000][3];
    int num_polys = 0;
    int k = 0; // k index+1/4 = num_polys
    int data_polys[] = new int[5000]; // 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 = 9; // then 17  theta
    int tpts = 16; // then 32  phi
    double dtheta = pi/8.0; // then pi/16.0
    double dphi   = pi/8.0; // then 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_curv16.java writing hcyl1_curv16.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;
        System.out.println("theta="+theta+"  "+degrees(theta)+" degrees");
	//makeMx(phi);
	//makeMy(phi);
        makeMz(phi);
        vp[0] = x;
	vp[1] = y;
	vp[2] = z;
        System.out.println("un rotated vp about Z is (x="+
 		           vp[0]+", y="+vp[1]+", z="+vp[2]+")");
        matvec(Mz, vp, vr);
        System.out.println("rotating vp about Z is vr= (x"+
 		           vr[0]+", y="+vr[1]+", z="+vr[2]+")");
	xp = vr[0];
	yp = vr[1];
	zp = vr[2];
        if(i==0) // debug
	{
	  System.out.println("yp offest rc*Math.cos(phi)="+(rc*Math.cos(phi)));
	  System.out.println("xp offest rc*Math.sin(phi)="+(rc*Math.sin(phi)));
	}				    
        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++;
        // debug
	System.out.println("outside phi="+phi+", xp="+xp+",yp="+yp+",zp="+zp);
	System.out.println("pti="+pti+", z="+z);
	System.out.println("sin(phi)="+Math.sin(phi)+", cos(phi)="+Math.cos(phi));
        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");
    //
    //     40-39-38-37-36 (tpts-1)*pts=1    
    //      | /| /| /| /|
    //      |/ |/ |/ |/ |
    //      5- 4- 3- 2- 1  pt0+1
    //      | /| /| /| /|
    //      |/ |/ |/ |/ |
    //     10- 9- 8- 7- 6  pt1+1
    //      | /| /| /| /|
    //      |/ |/ |/ |/ |
    //     15-14-13-12-11  pt1+1
    //

    for(int j=0; j<tpts-1; j++) // curv
    {
      pt0 = j*pts;
      pt1 = (j+1)*pts;	
      System.out.println("generate triangles outside pt0="+pt0+", pt1="+pt1);
      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");
	// makeMx(phi);
        //makeMy(phi);
        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;
    }

    // do bottom
    //
    //   1-36-31-26-21-16-11- 6- 1  pt0+1
    //   | /| /| /| /| /| /| /| /|
    //   |/ |/ |/ |/ |/ |/ |/ |/ |
    //  41-75-71- 65-61-56-51-46-41  pt1+1
    //    
    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_curv16.dat", data_points, num_points, data_polys, num_polys);
    
  } // end class hcyl1_curv16

  //
  //   | 1   0     0   |  | cosb  0  sinb |  | cosc  -sinc  0 |  | X |
  //   |               |  |               |  |                |  |   |
  //   | 0 cosa  -sina |  |   0   1   0   |  | sinc   cosc  0 |  | Y |
  //   |               |  |               |  |                |  |   |
  //   | 0 sina  cosa  |  | -simb 0  cosb |  |   0     0    1 |  | Z |
  //
  // rotate about X axis,      Y axis,           Z axis
  //

  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_curv16();
  }
} // end  hcyl1_curv16.java