/****************************************************
*            REDUCTION  OF  CONICALS                *
*     of equation: ax^2+2bxy+cy^2+2dx+2ey+f=0       *
* ------------------------------------------------- *
* The program, knowing the coefficients a, b, c, d, *
* e and f of the cartesian equation, finds the      *
* caracteristics of the reduction elements:         *
* type of conical  (hyperbola,parabola,ellipse or   *
* circle) and center position,  focus position(s),  *
* axes and excentricity.                            *
* ------------------------------------------------- *
* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp *
* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"  *
* [BIBLI 05].                                       *
* ------------------------------------------------- *
* SAMPLE RUN:                                       *
*                                                   *
*  REDUCTION OF CONICALS                            *
*                                                   *
*    a  = 1                                         *
*    2b = 2                                         *
*    c  = 1                                         *
*    2d = -13                                       *
*    2e = -11                                       *
*    f  = 32                                        *
*                                                   *
*  Type: Parabola                                   *
*                                                   *
*  Center:  x=0.999999  y=5.000000                  *
*  Symmetry direction:  x=1  y=-1                   *
*  Focus:   x=1.124999  y=4.875000                  *
*  Parameter:           0.35355339                  *
*                                                   *
*                       C++ version by J-P Moreau.  *
*                           (www.jpmoreau.fr)       *
****************************************************/
#include <stdio.h>
#include <math.h>


// LABEL: result

     //input variables:
     double a,b,c,d,e,f;   // coefficients of cartesian equation

     /*output variables:
        typ: type of conical (integer)
             1: ellipse
             2: hyperbola
             3: parabola
             4: circle
             5: line
             6: two lines
             7: one point
             8: no conical at all.
        ex: excentricity (or parameter for a parabola),
        xc,yc: coordinates of center,
        la,lb: axis half-length for an ellipse (la, radius for a circle),
        xf1,yf1,xf2,yf2: focus coordinates (only one focus for a parabola),
        xs1,ys1,xs2,ys2: summit coordinates for an hyperbola,
        xv1,yv1,xv2,yv2: vector directions of symmetry axes
                         (only one direction for a parabola). 
     */
     double ex,xc,yc,la,lb,xf1,yf1,xf2,yf2,xs1,ys1,xs2,ys2,xv1,yv1,xv2,yv2;

     // other internal variables
     double delta,u,l,m,x2,y2;
     int i,j,typ;


  void Parabola()  {
    double x,y;
    typ=3;
    if (c==0 && a==0) { //the parabola does not exist or is degenerated into a line
      if (d==0 && e==0) typ=8;else typ=5;  //8=no conical, 5=line
      return;
    }
    if (a==0) {
      x2=1; y2=0; l=0; m=1;
    }
    else {
      if (a<0) {
        f=-f; e=-e; d=-d;
        c=-c; b=-b; a=-a;
      }
      l=sqrt(a); m=sqrt(c);
      if (b<0) m=-m;
      u=sqrt(a+c);
      x2=m/u; y2=-1/u;
      f=f*u;c=(a+c)*u;
      u=d*m-e*l; e=d*l+e*m; d=u;
    }
    if (d==0) {
      if (e*e<c*f) typ=8; else typ=6; //8=no conical, 6=two lines
      return;
    }
    else {
      x=(e*e-c*f)/2/c/d; y=-e/c;
      xc=x*x2-y*y2; yc=y*x2+x*y2;
      ex=-d/c;
      xf1=xc+ex*x2/2; yf1=yc+ex*y2/2;
      xv1=m; yv1=-l;
    }
  }  //Parabola()
                 
  
  void Hyperbola()  {
    typ=2;
    la=sqrt(-fabs(f)/l); 
	lb=sqrt(fabs(f)/m);
    if (f<0) {
      u=la; la=lb; lb=u;
      u=x2; x2=-y2; y2=u;
      u=xv1; xv1=-yv1; yv1=u;
    }
    xv2=-yv1; yv2=xv1;
    u=sqrt(la*la+lb*lb);
    xs1=xc+x2*la; ys1=yc+y2*la;
    xs2=xc-x2*la; ys2=yc-y2*la;
    xf1=xc+x2*u; yf1=yc+y2*u;
    xf2=xc-x2*u; yf2=yc-y2*u;
    ex=u/la;
  }

  void Ellipse() {
    typ=1;
    if (f*l>0) {typ=8; return;}  //8=no conical
    la=sqrt(-f/l); lb=sqrt(-f/m);
    if (l<0) {
      u=la; la=lb; lb=u;
      u=x2; x2=-y2; y2=u;
      u=xv1; xv1=-yv1; yv1=u;
    }
    xv2=-yv1; yv2=xv1;
    u=sqrt(la*la-lb*lb);
    xf1=xc+x2*u; yf1=yc+y2*u;
    xf2=xc-x2*u; yf2=yc-y2*u;
    ex=u/la;
  }


void main()  {

  printf("\n  REDUCTION OF CONICALS\n\n");
  printf("    a  = "); scanf("%lf",&a);
  printf("    2b = "); scanf("%lf",&b);
  printf("    c  = "); scanf("%lf",&c);
  printf("    2d = "); scanf("%lf",&d);
  printf("    2e = "); scanf("%lf",&e);
  printf("    f  = "); scanf("%lf",&f);
  printf("\n");
  b=b/2; d=d/2; e=e/2;
  delta=a*c-b*b;
  if (delta==0) Parabola();
  else {
    xc=(b*e-d*c)/delta;
    yc=(b*d-a*e)/delta;
    f += a*xc*xc+2*b*xc*yc+c*yc*yc+2*d*xc+2*e*yc;
    if (f==0) {
      if (delta>0) typ=7; else typ=6;  //7=one point, 6=two lines
      goto result;
    }
    u=sqrt((a-c)*(a-c)+4*b*b);
    l=(a+c-u)/2; m=(a+c+u)/2;
    if (a==c && b==0)
      if (f*a>=0) {
        typ=8; // no conical
        goto result;
      }
      else {
        typ=4; // circle
        la=sqrt(-f/a);
        ex=1;
        goto result;
      }
    if (a<c && b==0) {
      x2=1; y2=0; xv1=1; yv1=0;
    }
    else {
      xv1=b; yv1=1-a;
      u=sqrt(xv1*xv1+yv1*yv1);
      x2=xv1/u; y2=yv1/u;
    }
    if (delta<0) Hyperbola(); else Ellipse();
} // else if delta=0

  // print results:
result: switch(typ) {
    case 1: {
        printf("  Type: Ellipse\n\n");
        printf("  Center                :  x=%9.6f  y=%9.6f\n", xc, yc);
        printf("  Direction big axis    :  x=%5.2f  y=%5.2f\n", xv1, yv1);
        printf("  Direction small axis  :  x=%5.2f  y=%5.2f\n", xv2, yv2);
        printf("  Half length big axis  : %9.6f\n", la);
        printf("  Half length small axis: %9.6f\n", lb);
        printf("  First focus           :  x=%9.6f  y=%9.6f\n", xf1, yf1);
        printf("  Second focus          :  x=%9.6f  y=%9.6f\n", xf2, yf2);
        printf("  Excentricity          :  %10.8f\n", ex);
		break;
			}
    case 2: {
        if (a+c==0) printf("  Type: Equilateral Hyperbola\n\n");
               else printf("  Type: Hyperbola\n\n");
        printf("  Center                :  x=%9.6f  y=%9.6f\n", xc, yc);
        printf("  Direction first axis  :  x=%5.2f  y=%5.2f\n", xv1, yv1);
        printf("  Direction second axis :  x=%5.2f  y=%5.2f\n", xv2, yv2);
        printf("  First summit          :  x=%9.6f  y=%9.6f\n", xs1, ys1);
        printf("  Second summit         :  x=%9.6f  y=%9.6f\n", xs2, ys2);
        printf("  First focus           :  x=%9.6f  y=%9.6f\n", xf1, yf1);
        printf("  Second focus          :  x=%9.6f  y=%9.6f\n", xf2, yf2);
        printf("  Excentricity          :  %10.8f\n", ex);
		break;
			}
    case 3: {
        printf("  Type: Parabola\n\n");
        printf("  Center                :  x=%9.6f  y=%9.6f\n", xc, yc);
        printf("  Symmetry direction    :  x=%5.2f  y=%5.2f\n", xv1, yv1);
        printf("  Focus                 :  x=%9.6f  y=%9.6f\n", xf1, yf1);
        printf("  Parameter             :  %10.8f\n", ex);
		break;
			}
    case 4: {
        printf("  Type: Circle\n\n");
        printf("  Center                :  x=%9.6f  y=%9.6f\n", xc, yc);
        printf("  Radius                :  %9.6f\n", la);
        printf("  Excentricity          :  %10.8f\n", ex);
		break;
			}
    case 5: printf("  Conical degenerated into a line.\n");
	        break;
    case 6: printf("  Conical degenerated into two lines.\n");
	        break;
    case 7: printf("  Conical degenerated into a single point.\n");
	        break;
    case 8: printf("  Not a conical !\n");

		} //end switch(typ)

  printf("\n\n");

}

//end of file conical.cpp