// sphere_div.java  divide a tetrahedron   
//                           octahedron, write a .dat file   
//               compile with datread.java and datwrite.java             

// sphere_div sphere_div3.dat 3  # last number is split, max=4 

//  tetrahedron   4 vertices     4 faces     6 edges  
//  split 1      10 vertices    16 faces    24 edges  
//  split 2      34 vertices    64 faces    96 edges  
//  split 3     130 vertices   256 faces   384 edges  
//  split 4     514 vertices  1024 faces  1536 edges  

//  octahedron    6 vertices     8 faces    12 edges  
//  split 1      18 vertices    32 faces    48 edges  
//  split 2      66 vertices   128 faces   192 edges  
//  split 3     258 vertices   512 faces   768 edges  
//  split 4    1026 vertices  2048 faces  3072 edges  

// num_points = vertices,  num_polys = faces 

public class sphere_div
{
  //  dpts [][]

  int num_points = 0; // for datread, datwrite and datclean 
  int num_polys = 0;
  int num_edges = 0;
  // dpts data_points[] = new dpts[2050]; // each has .x .y .z 
  double data_points[][] = new double[2050][3]; // each has 0 1 2 
  int  data_verts[] = new int[8200]; // 4*num_polys is count=3 index 
  int  link_verts[][] = new int[2050][3]; // links to three edges 
  int  edge_list[][] = new int[3100][3]; // old_vert, mid_vert, poly 
  double size;
  int status = -1;
  int splits = 1;
  double Pi = 3.141592653589793238462643383279502884197;
  double vertices[][] = new double[6][3]; // initial tetrahedron 

  //  void tetrahedron_compute(); defined below
  //  void octahedron_compute();
  //  void split(int splits);
  //  double d3len(double a[], double b[]);
  //  void vertex_compute();
  //  void polygon(int a, int b, int c);
  //  int insert_edge(int a, int b);
  //  void split_poly(int i);

  public sphere_div(String out_name, String s_splits)
  {
    splits = Integer.parseInt(s_splits);
    System.out.println("sphere_div.java splits="+splits+", and write to "+
                       out_name);
  
    split(splits);  // do the work, set up for datwrite 
  
    new datwrite(out_name, data_points, num_points, data_verts, num_polys);
    System.out.println("sphere_div.java has written "+ out_name);
  } // end sphere_div constructor

  double d3len(double a[], double b[])
  {
    double x1, y1, z1, x2, y2, z2;
    x1 = a[0];  y1 = a[1];  z1 = a[2];
    x2 = b[0];  y2 = b[1];  z2 = b[2];
    return Math.sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
  }

  void tetrahedron_compute()
  {
    double phia, the120, the, r;
    double phiaa  = -19.471220333; // generator 
  
    // based on unit sphere x*x+y*y+z*z=1           
    // x=r*cos(the)*cos(phi)           the 0 to 2Pi 
    // y=r*sin(the)*cos(phi)    phi -Pi/2 to Pi/2   
    // z=r*sin(phi)          radians=Pi*degrees/180 
    r = 1.0;
    phia = Pi*phiaa/180.0; // 1 set of three points 
    vertices[0][0] = 0.0;
    vertices[0][1] = 0.0;
    vertices[0][2] = r;
    System.out.println("i="+(0)+", x="+vertices[0][0]+", y="+vertices[0][1]+
                       ", z="+vertices[0][2]);
    the120 = Pi*120.0/180;
    the = 0.0;
    for(int i=1; i<4; i++)
    {
      vertices[i][0]=r*Math.cos(the)*Math.cos(phia);
      vertices[i][1]=r*Math.sin(the)*Math.cos(phia);
      vertices[i][2]=r*Math.sin(phia);
      System.out.println("i="+i+", x="+vertices[i][0]+", y="+vertices[i][1]+
                         ", z="+vertices[i][2]);
      the = the+the120;
    }

    for(int i=0; i<4; i++)
    {
      data_points[num_points][0] = vertices[i][0];
      data_points[num_points][1] = vertices[i][1];
      data_points[num_points][2] = vertices[i][2];
      num_points++;
    }

    // map vertices to 4 faces 
    polygon(0,1,2);
    polygon(0,2,3);
    polygon(0,3,1);
    polygon(1,2,3);
  }

  void octahedron_compute()
  {
    double phia, the90, the, r;
    r = 1.0;
    phia = 0.0; // 1 set of four points 
    the90 = Pi*90.0/180;
    vertices[0][0]=0.0;
    vertices[0][1]=0.0;
    vertices[0][2]=r;
    System.out.println("i=0, x="+vertices[0][0]+", y="+vertices[0][1]+
                       ", z="+vertices[0][2]);
    vertices[5][0]=0.0;
    vertices[5][1]=0.0;
    vertices[5][2]=-r;

    the = 0.0;
    for(int i=1; i<5; i++)
    {
      vertices[i][0]=r*Math.cos(the)*Math.cos(phia);
      vertices[i][1]=r*Math.sin(the)*Math.cos(phia);
      vertices[i][2]=r*Math.sin(phia);
      System.out.println("i="+i+", x="+vertices[i][0]+", y="+vertices[i][1]+
                         ", z="+vertices[i][2]);
      the = the+the90;
    }
    System.out.println("i=5, x="+vertices[5][1]+", y="+vertices[5][1]+
                       ", z="+vertices[5][2]);
    System.out.println("num_points="+num_points+", num_polys="+num_polys+
                       ", num_edges="+num_edges);
    num_polys = 0;
    System.out.println("num_points="+num_points+", num_polys="+num_polys +
                       ", num_edges="+num_edges);

    for(int i=0; i<6; i++)
    {
      data_points[num_points][0] = vertices[i][0];
      data_points[num_points][1] = vertices[i][1];
      data_points[num_points][2] = vertices[i][2];
      num_points++;
    }

    // map vertices to 8 faces 
    polygon(0,1,2);
    polygon(0,2,3);
    polygon(0,3,4);
    polygon(0,4,1);
    polygon(5,2,1);
    polygon(5,3,2);
    polygon(5,4,3);
    polygon(5,1,4);
  }

  void polygon(int a, int b, int c)
  {
    int i;

    // insert next polygon 
    i = 4*num_polys;
    data_verts[i]   = 3;
    data_verts[i+1] = a+1; // .dat file starts with 1, not zero 
    data_verts[i+2] = b+1;
    data_verts[i+3] = c+1;
    link_verts[num_polys][0] = insert_edge(a,b);
    link_verts[num_polys][1] = insert_edge(b,c);
    link_verts[num_polys][2] = insert_edge(c,a);
    num_polys++;
  }

  int insert_edge(int a, int b)
  {
    int first, second;

    if(a<b) { first=a; second=b; }
    else    { first=b; second=a; }
    for(int i=0; i<num_edges; i++)
    {
      if(edge_list[i][0]==first && edge_list[i][1]==second) return i;
    }
    edge_list[num_edges][0] = first; // raw, not +1 
    edge_list[num_edges][1] = second;
    edge_list[num_edges][2] = 0; // for mid point vertex 
    num_edges++;
    return num_edges-1;
  }

  void split_poly(int i) // make four faces out of one 
  {
    int a, b, c, dl, d, el, e, fl, f; // original a,b,c  mids  d,e,f 
    // new d,e,f   a,d,f   b,e,d  c,f,e   
    a  = data_verts[4*i+1]-1;
    b  = data_verts[4*i+2]-1;
    c  = data_verts[4*i+3]-1;
    dl = link_verts[i][0];
    d  = edge_list[dl][2];
    el = link_verts[i][1];
    e  = edge_list[el][2];
    fl = link_verts[i][2];
    f  = edge_list[fl][2];
    System.out.println("a="+a+", b="+b+", c="+c+", d="+d+", e="+e+", f="+f+
                       ", dl="+dl+", el="+el+", fl="+fl);
    data_verts[4*i+1] = d+1;
    data_verts[4*i+2] = e+1;
    data_verts[4*i+3] = f+1;
    data_verts[4*num_polys]   = 3;
    data_verts[4*num_polys+1] = a+1;
    data_verts[4*num_polys+2] = d+1;
    data_verts[4*num_polys+3] = f+1;
    num_polys++;
    data_verts[4*num_polys]   = 3;
    data_verts[4*num_polys+1] = b+1;
    data_verts[4*num_polys+2] = e+1;
    data_verts[4*num_polys+3] = d+1;
    num_polys++;
    data_verts[4*num_polys]   = 3;
    data_verts[4*num_polys+1] = c+1;
    data_verts[4*num_polys+2] = f+1;
    data_verts[4*num_polys+3] = e+1;
    num_polys++;
  }

  void split(int splits)
  {
    int j, k;
    double mid[] = new double[3];
    double origin[] = {0.0, 0.0, 0.0};
    double len, len_one;
  
    // tetrahedron_compute(); 
    octahedron_compute();
    System.out.println("num_points="+num_points+", num_polys="+num_polys+
                       ", num_edges="+num_edges);
    for(int s=0; s<splits; s++)
    {
      System.out.println("splits="+splits+", s="+s);
      // initial .dat structure built, now split 
      // build link_verts and edge list 
      num_edges = 0;
      for(int i=0; i<num_polys; i++)
      {
        link_verts[i][0] = insert_edge(data_verts[4*i+1]-1, data_verts[4*i+2]-1);
        link_verts[i][1] = insert_edge(data_verts[4*i+2]-1, data_verts[4*i+3]-1);
        link_verts[i][2] = insert_edge(data_verts[4*i+3]-1,data_verts[4*i+1]-1);
      }
      for(int i=0; i<num_edges; i++)
      {
        System.out.println("edge="+i+","+edge_list[i][0]+","+ edge_list[i][1]);
      }

      for(int i=0; i<num_edges; i++) // split edges 
      {
        j = edge_list[i][0];
        k = edge_list[i][1];
        mid[0]=(data_points[j][0]+data_points[k][0])/2.0; // x 
        mid[1]=(data_points[j][1]+data_points[k][1])/2.0; // y 
        mid[2]=(data_points[j][2]+data_points[k][2])/2.0; // z 
        len = d3len(mid, origin);
        mid[0] = mid[0]/len;
        mid[1] = mid[1]/len;
        mid[2] = mid[2]/len;
        System.out.println("midx="+mid[0] +", y="+mid[1]+", z="+mid[2]+
                           ", len="+len);
                
        edge_list[i][2] = num_points;
        data_points[num_points][0] = mid[0];
        data_points[num_points][1] = mid[1];
        data_points[num_points][2] = mid[2];
        num_points++;
      }

      System.out.println("num_points="+num_points+", num_polys="+num_polys+
                         ", num_edges="+num_edges);
      for(int i=0; i<num_edges; i++)
      {
        System.out.println("edge="+i+",  v1="+edge_list[i][0]+
                           ", v2="+edge_list[i][1]+", mid="+edge_list[i][2]);
      }
      for(int i=0; i<num_polys; i++)
      {
        System.out.println("link="+i+", d="+link_verts[i][0]+
                           ",  e="+link_verts[i][1]+", f="+link_verts[i][2]);
      }


      j = num_polys; // only loop through old 
      for(int i=0; i<j; i++) // split polys 
      {
        split_poly(i);
      }
      System.out.println("num_points="+num_points+", num_polys="+num_polys+
                         ", num_edges="+num_edges);
             
      for(int i=0; i<num_polys; i++)
      {
  	System.out.println("polygon("+(data_verts[4*i+1]-1)+","+
                          (data_verts[4*i+2]-1)+","+(data_verts[4*i+3]-1)+")");
      }
      System.out.println("num_points="+num_points+", num_polys="+num_polys+
                         ", num_edges="+num_edges);
    }
    num_edges = num_edges*4;
    System.out.println("num_points="+num_points+", num_polys="+num_polys+
                       ", num_edges="+num_edges);
  } // end split

  public static void main (String[] args) // file name, number of splits
  {
      new sphere_div(args[0],args[1]);
  }
} // end sphere_div.java