'*****************************************************
'*            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.: "Mathematiques 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                  *
'*                                                   *
'*                     Basic version by J-P Moreau.  *
'*                           (www.jpmoreau.fr)       *
'*****************************************************
defdbl a-h,o-z
defint i-n

'    input variables:
'    a,b,c,d,e,f: DOUBLE;   (coefficients of cartesian equation)

'    output variables:
'       ityp: 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,
'       xla,xlb: 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).

'    other internal variables:
'    delta,u,xl,xm,x2,y2: DOUBLE;
'    i,j,ityp: INTEGER;


'begin main program

  cls
  print
  print "  REDUCTION OF CONICALS"
  print
  input "    a  = ", a
  input "    2b = ", b
  input "    c  = ", c
  input "    2d = ", d
  input "    2e = ", e
  input "    f  = ", f
  print
  b=b/2: d=d/2: e=e/2
  delta=a*c-b*b
  if delta=0 then
    gosub 1000       'call Parabola
  else
    xc=(b*e-d*c)/delta
    yc=(b*d-a*e)/delta
    f = f + a*xc*xc+2*b*xc*yc+c*yc*yc+2*d*xc+2*e*yc
    if f=0 then
      if delta>0 then
        ityp=7        'one point
      else
        ityp=6        'two lines
        goto 100
      end if
    end if
    u=sqr((a-c)*(a-c)+4*b*b)
    xl=(a+c-u)/2: xm=(a+c+u)/2
    if a=c and b=0 then
      if f*a>=0 then
        ityp=8   'no conical
        goto 100
      else
        ityp=4   'circle
        xla=sqr(-f/a)
        ex=1#
        goto 100
      end if
    end if
    if a<c and b=0 then
      x2=1#: y2=0#: xv1=1#: yv1=0#
    else
      xv1=b: yv1=1#-a
      u=sqr(xv1*xv1+yv1*yv1)
      x2=xv1/u: y2=yv1/u
    end if
    if delta<0 then
      gosub 2000     'call Hyperbola
    else
      gosub 3000:    'call Ellipse
    end if
  end if 'else if delta=0

' print results
100  if ityp=1 then
        print "  Type: Ellipse"
        print
        print "  Center                :  x="; xc;"  y="; yc
        print "  Direction big axis    :  x="; xv1;"  y="; yv1
        print "  Direction small axis  :  x="; xv2;"  y="; yv2
        print "  Half length big axis  :  "; xla
        print "  Half length small axis:  "; xlb
        print "  First focus           :  x="; xf1;"  y="; yf1
        print "  Second focus          :  x="; xf2;"  y="; yf2
        print "  Excentricity          : "; ex
     end if
     if ityp=2 then
        if a+c=0 then
          print "  Type: Equilateral Hyperbola"
        else
          print "  Type: Hyperbola"
        end if
        print
        print "  Center                :  x="; xc;"  y="; yc
        print "  Direction first axis  :  x="; xv1;"  y="; yv1
        print "  Direction second axis :  x="; xv2;"  y="; yv2
        print "  First summit          :  x="; xs1;"  y="; ys1
        print "  Second summit         :  x="; xs2;"  y="; ys2
        print "  First focus           :  x="; xf1;"  y="; yf1
        print "  Second focus          :  x="; xf2;"  y="; yf2
        print "  Excentricity          : "; ex
     end if
     if ityp=3 then
        print "  Type: Parabola"
        print
        print "  Center                :  x="; xc; "  y="; yc
        print "  Symmetry direction    :  x="; xv1;"  y="; yv1
        print "  Focus                 :  x="; xf1;"  y="; yf1
        print "  Parameter             : "; ex
     end if
     if ityp=4 then
        print "  Type: Circle"
        print
        print "  Center                :  x="; xc;"  y="; yc
        print "  Radius                :  "; xla
        print "  Excentricity          : "; ex
     end if
     if ityp=5 then print "  Conical degenerated into a line."
     if ityp=6 then print "  Conical degenerated into two lines."
     if ityp=7 then print "  Conical degenerated into a single point."
     if ityp=8 then print "  Not a conical !"

  print

END 'of main program


'Subroutine Parabola
1000 ityp=3
    if c=0 and a=0 then 'the parabola does not exist or is degenerated into a line
      if d=0 and e=0 then
        typ=8 'no conical
      else
        typ=5 'line
        return
      end if
    end if
    if a=0 then
      x2=1: y2=0: xl=0: xm=1
    else
      if a<0 then
        f=-f: e=-e: d=-d
        c=-c: b=-b: a=-a
      end if
      xl=sqr(a): xm=sqr(c)
      if b<0 then xm=-xm
      u=sqr(a+c)
      x2=xm/u: y2=-xl/u
      f=f*u:c=(a+c)*u
      u=d*xm-e*xl: e=d*xl+e*xm: d=u
    end if
    if d=0 then
      if e*e<c*f then
        typ=8
      else
        typ=6  'two lines
        return
      end if
    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=xm: yv1=-xl
    end if
return


'Subroutine Hyperbola
2000 ityp=2
    xla=sqr(-abs(f)/xl): xlb=sqr(abs(f)/xm):
    if f<0 then
      u=xla: xla=xlb: xlb=u
      u=x2: x2=-y2: y2=u
      u=xv1: xv1=-yv1: yv1=u
    end if
    xv2=-yv1: yv2=xv1
    u=sqr(xla*xla+xlb*xlb)
    xs1=xc+x2*xla: ys1=yc+y2*xla
    xs2=xc-x2*xla: ys2=yc-y2*xla
    xf1=xc+x2*u: yf1=yc+y2*u
    xf2=xc-x2*u: yf2=yc-y2*u
    ex=u/xla
return


'Subroutine Ellipse
3000  ityp=1
    if f*xl>0 then
      typ=8  'no conical
      return
    end if
    xla=sqr(-f/xl): xlb=sqr(-f/xm)
    if xl<0 then
      u=xla: xla=xlb: xlb=u
      u=x2: x2=-y2: y2=u
      u=xv1: xv1=-yv1: yv1=u
    end if
    xv2=-yv1: yv2=xv1
    u=sqr(xla*xla-xlb*xlb)
    xf1=xc+x2*u: yf1=yc+y2*u
    xf2=xc-x2*u: yf2=yc-y2*u
    ex=u/xla
return

'end of file conical.bas