// ellipse3.java    compile    javac -cp . ellipse3.java
//                  execute    java  -cp . ellipse3
//
//                 . x4,y4  vertical axis  must be perpendicular 
//            .    |    .                  to horizontal axis
//          .      |       .
//   x1,y1 .----------------. x2,y2  horizontal axis  does not have to be
//          .      |       .                    parallel to the X axis
//            .    |    .      xc,yc is center
//                 . x3,y3      x,y  is point on ellipse3

public class ellipse3
{
  double x1 =  0.0; // rotated test
  double y1 =  0.0; // center not 0,0  is 2,2
  double x2 =  4.0;
  double y2 =  4.0;
  double x3 =  1.0;
  double y3 =  3.0;
  double x4 =  3.0;
  double y4 =  1.0;

  public ellipse3()
  {
    double a, b, xc, yc, x, y, chk, t, dt, area;
    double trot, r, ang, xr, yr;
    double dx12, dy12, dx34, dy34, len12, len34, ang12, ang34, dot13, cross13;
    int nt;
    int num_points = 0; // for datwrite
    int num_polys = 0;
    int poly = 0;
    int polys[] = new int[1000];
    double points[][] = new double[200][3];
    
    System.out.println("ellipse3.java running");
    System.out.println("x1,y1 ="+x1+","+y1+"  one end of major axis");
    System.out.println("x2,y2 ="+x2+","+y2+"  other end of major axis");
    System.out.println("x3,y3 ="+x3+","+y3+"  one end of minor axis");
    System.out.println("x4,y4 ="+x4+","+y4+"  other end of minor axis");
    xc = (x1+x2)/2.0;
    yc = (y3+y4)/2.0;
    System.out.println("xc,yc ="+xc+","+yc+"  center of ellipse3");
    a = len(x1,y1,x2,y2)/2.0;
    b = len(x3,y3,x4,y4)/2.0;
    System.out.println("a = "+a+" length of major axis");
    System.out.println("b = "+b+" length of minor axis");
    area = Math.PI*a*b;
    System.out.println("area = "+area);
    
    if(Math.abs(y1-y2)+Math.abs(x3-x4)<0.001) // horizontal vertical
    {
      System.out.println("check point x,y on ellipse3 (x*x)/(a*a) + (y*y)/(b*b) = 1.0");
      chk = (x1*x1)/(a*a) + (y1*y1)/(b*b);
      System.out.println("x1,y1 chk="+chk);
      chk = (x2*x2)/(a*a) + (y2*y2)/(b*b);
      System.out.println("x2,y2 chk="+chk); 
      chk = (x3*x3)/(a*a) + (y3*y3)/(b*b);
      System.out.println("x3,y3 chk="+chk);
      chk = (x4*x4)/(a*a) + (y4*y4)/(b*b);
      System.out.println("x4,y4 chk="+chk);
      System.out.println("given x, y =sqrt((1.0-((x*x)/(a*a))) * (b*b)) = ell(x,a,b)");
      System.out.println("generate some points, could be in a loop");
      x = (x1);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = (x1+x2/8.0);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = (x1+x2/4.0);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = (x1+x2/2.0);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = xc;
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = (x1/2.0+x2);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = (x1/4.0+x2);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = (x1/8.0+x2);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      x = (x2);
      y = ell(x,a,b);
      System.out.println("point on ellipse3 x,y = "+x+","+y);
      System.out.println(" ");
    }
    System.out.println("check angles ");
    System.out.println("45 deg = "+Math.PI/4.0+" rad");
    dx12 = x2-x1;
    dy12 = y2-y1;
    ang12 = Math.atan2(dy12,dx12);
    dx34 = x4-x3;
    dy34 = y4-y3;
    ang34 = Math.atan2(dy34, dx34);
    System.out.println("ang12 = "+ang12+"  ang34 = "+ang34);
    len12 = len(x1,y1,x2,y2);
    len34 = len(x3,y3,x4,y4);
    System.out.println("length 12 = "+len12+"  length34 = "+len34);
    dot13 = dot(dx12, dy12, dx34, dy34)/(len12*len34);
    System.out.println("normalized dot product = "+dot13);
    cross13 = cross(dx12, dy12, dx34, dy34)/(len12*len34);
    System.out.println("cross product = "+cross13);

    System.out.println("use x=xc+a*cos(t) y=yc+b*sin(t) ");
    for(int i=0; i<200; i++)
    {
      points[i][2] = 0.0; // no z
    }
    dt = Math.PI/16.0;
    trot = Math.PI/4.0;
    System.out.println("dt = "+dt+"  trot = "+trot);
    t = 0.0;
    nt = 32;
    points[0][0] = x1;
    points[0][1] = y1;
    points[1][0] = x2;
    points[1][1] = y2;
    points[2][0] = x3;
    points[2][1] = y3;
    points[3][0] = x4;
    points[3][1] = y4;
    num_points = 4;
    polys[poly++] = 2; // major axis
    polys[poly++] = 1;
    polys[poly++] = 2;
    polys[poly++] = 2; // minor axis
    polys[poly++] = 3;
    polys[poly++] = 4;

    for(int i=0; i<nt; i++)
    {
      x = px(a, xc, t);
      y = py(b, yc, t);
      r = len(x, y, xc, yc);
      ang = Math.atan2(y-yc,x-xc);
      System.out.println("x,y = "+x+","+y);
      System.out.println("r = "+r+"  ang = "+ang);
      xr = xc + r*Math.cos(trot+ang);
      yr = yc + r*Math.sin(trot+ang);
      System.out.println("xr,yr = "+xr+","+yr);
      points[num_points][0] = xr;
      points[num_points][1] = yr;
      num_points++;
      if(i>0)
      {
        polys[poly++] = 2;
        polys[poly++] = num_points-1;
        polys[poly++] = num_points;
      }
      t = t+dt;
    }
    polys[poly++] = 2; // close last to first
    polys[poly++] = num_points;
    polys[poly++] = 5;
    
    num_polys = poly/3; // all lines
    
    new datwrite("ellipse3.dat", points, num_points, polys, num_polys);
    System.out.println("ellipse3.java wrote allipse.dat");
    System.out.println("ellipse3.java finished"); 
  }

  double len(double xa, double ya, double xb, double yb) // length a to b
  {
    return Math.sqrt((xa-xb)*(xa-xb) + (ya-yb)*(ya-yb));
  } // end len

  double ell(double x, double a, double b)
  {
    return Math.sqrt((1.0-((x*x)/(a*a))) * (b*b));
  } // end ell

  double px(double a, double xc, double t)
  {
    return xc+a*Math.cos(t); 
  } // end px

  double py(double b, double yc, double t)
  {
    return yc+b*Math.sin(t); 
  } // end px

  double dot(double dx1, double dy1, double dx2, double dy2)
  {                               // not normalized
    return dx1*dx2 + dy1*dy2;
  } // end dot

  double cross(double dx1, double dy1, double dx2, double dy2)
  {                               // not normalized
    return dx1*dy2 - dy1*dx2;
  } // end cross

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