HP Forums

Full Version: Conic Section...
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
This code intention is to bring conic sections utility to Prime via Advanced Graphing and also as a standalone application. The application works(it identifies the conic section and calculates its parameters) .

Ecc= Eccentricity
Fcs=Foci
Vtx=Vertex
Dtx=Directrix expression
Sym=Symmetry expression
Asy=Asymptotes expressions
Ctr=Center
LRe=Latus Rectum Length
Radius= Radius
eqRθ = Simplified Canonical Expression; Rotating Angle






usage;
plot(X,Y) in Adv. Graphing if there exist a conic in the plot;
(2*X*Y)+X^2+(Y^2)-3*X=0
vget()
or
Conics(expression in x,y) as standalone application.

examples;
Conics((2*x*y)+x^2+(y^2)-3*x)
Conics(x^2+y^2-2);
Conics(x^2+y-2);


Code:

#pragma mode( separator(.,;) integer(h32) )
Conics();
EXPORT vget()
BEGIN
LOCAL I,func,fund;
LOCAL funcs:={}; 
FOR I FROM 1 TO 10 DO
 IF Advanced_Graphing.ISCHECK(I MOD 10) THEN
    func:=LOWER(string(expr(CAS.cat("'V'",(I MOD 10))))); 
    fund:=CAS.expr(LEFT(func,(INSTRING(func,"=")-1)) +"-1*"+ "("+ (RIGHT(func,(DIM(func)-INSTRING(func,"=")))) +")"); 
    IF TYPE(denom(normal(fund)))==8 THEN fund:=numer(normal(fund)); END;
    fund:=approx(CAS.expand(fund)); 
    fund:=(STRING(fund)); 
IF SIZE(CAS.coeff(fund,"'x'"))==1 AND SIZE(CAS.coeff(fund,"'y'"))==1 AND TYPE(lcoeff(CAS.lcoeff(fund,"'x'"),"'y'"))==0 THEN
       funcs[size(funcs)+1]:=CAS.expr(fund);
END;    
IF SIZE(CAS.coeff((fund),"'x'"))<4 AND SIZE(CAS.coeff((fund),"'y'"))<4 AND (TYPE(CAS.lcoeff(fund,"'x'"))==0 AND TYPE(CAS.lcoeff(fund,"'y'"))==0) THEN    
       funcs[size(funcs)+1]:=(CAS.expr(fund));
END; 
END;    
END;
IF SIZE(funcs)>=1 THEN Conics(CAS.exact(funcs[1])); END;
END;


#CAS
Conics(ff)
BEGIN
LOCAL Aq,Bq,Cq,Dq,Eqq,Fq;
LOCAL A,BB,Dx,Sx,LRa;
LOCAL g,c,AA,AAA,r1,ffx;
LOCAL LL,LLL,LLLL,LLLLL;
LOCAL aa,bb,MM,CC,DD,RRR;
IF AAngle<>1 THEN AAngle:=1; END;
LL:={}; LLL:={}; LLLL:={}; LLLLL:={}; aa:={}; bb:={}; M25:={}; M26:={}; r1:={};
IF TYPE(lcoeff(ff,x))==0 THEN A:=SIGN(lcoeff(ff,x)); ELSE A:=1; END;
IF has(ff,x) THEN aa:=A*(coeff(ff,x)); ELSE return ("NOT A CONIC1"); END;
IF has(ff,y) THEN bb:=A*(coeff(ff,y)); ELSE return ("NOT A CONIC2"); END;
IF SIZE(aa)==0 OR SIZE(bb)==0 THEN          return ("NOT A CONIC3"); END;
IF SIZE(aa)>3  OR SIZE(bb)>3  THEN          return ("NOT A CONIC4"); END;
IF (SIZE(aa)==2 AND SIZE(bb)==2) AND (aa[2])==0 THEN return ("NOT A CONIC5"); END;
IF (SIZE(aa)==2 AND SIZE(bb)==2) AND (TYPE(aa[1])==0 AND TYPE(bb[1])==0) THEN return ("line"); END;
IFTE( SIZE(aa)==2,aa={0,aa[1],aa[2]},Aq=aa[1]);
Aq:=aa[1];
IF has(aa[2],y) THEN BB:=coeff(aa[2],y); Bq:=BB[1]; Dq:=BB[2]; ELSE Bq:=0; Dq:=aa[2]; END;
IF has(aa[3],y) THEN BB:=coeff(aa[3],y); Fq:=BB[SIZE(BB)]; ELSE Fq:=aa[3]; END;
IFTE( SIZE(bb)==2,bb={0,bb[1],bb[2]},Cq=bb[1]);
Cq:=bb[1];
IF has(bb[2],x) THEN BB:=coeff(bb[2],x); Eqq:=BB[2]; ELSE Eqq:=bb[2]; END;
IF has(bb[3],x) THEN BB:=coeff(bb[3],x); Fq:=BB[SIZE(BB)]; ELSE Fq:=bb[3]; END;
MM:=sum(sum([[Aq, Bq/2, Dq/2],[Bq/2, Cq, Eqq/2],[Dq/2, Eqq/2, Fq]]));
IF TYPE(MM)<>0 THEN return ("NOT A CONIC6"); END;
MM:=DET([[Aq, Bq/2, Dq/2],[Bq/2, Cq, Eqq/2],[Dq/2, Eqq/2, Fq]]);
IF MM==0 AND ((Dq^2 +Eqq^2 > ((Aq+Cq)*Fq)) OR (Dq^2 +Eqq^2 < ((Aq+Cq)*Fq)) OR (Dq^2 +Eqq^2 == ((Aq+Cq)*Fq)))  THEN return ("line"); END; 

LL:={(Bq*Eqq-2*Cq*Dq)/(4*Aq*Cq-Bq^2),(-2*Aq*Eqq+Bq*Dq)/(4*Aq*Cq-Bq^2)};
IF Bq==0 THEN  BB:=(Aq*(Dq/(2*Aq))^2)+(Cq*(Eqq/(2*Cq))^2)-Fq; RRR:=({(1/(aa[1]/BB)),(1/(bb[1]/BB))}); LLL:=ABS(RRR); END;
IF Bq<>0 THEN 
IF Aq<>Cq THEN r:=(ATAN(Bq/(Aq-Cq)))/2; AA:=approx(COS(r)); AAA:=approx(SIN(r)); ELSE r:=pi/4; AA:=COS(r); AAA:=SIN(r); END; 
CC:=AA*'x'- AAA*'y'; DD:=AAA*'x'+AA*'y'; 
ffx:=simplify((Aq*(CC)^2)+(Bq*(CC)*(DD))+(Cq*(DD)^2)+Dq*(CC)+Eqq*(DD)+Fq); 
aa:=SIGN(lcoeff(ffx,x))*(coeff(ffx,x)); IFTE( SIZE(aa)==2,aa={0,aa[1],aa[2]},0);
IF has(aa[2],y) THEN BB:=coeff(aa[2],y); Bq:=BB[1]; ELSE Bq:=0; END;
IF Bq<>0 THEN
aa:=SIGN(lcoeff(ffx,x))*(coeff(ffx,x)); IFTE( SIZE(aa)==2,aa={0,aa[1],aa[2]},0);
Aq:=aa[1];
IF has(aa[2],y) THEN BB:=coeff(aa[2],y); Bq:=ROUND(BB[1],2); Dq:=BB[2]; ELSE Bq:=0; Dq:=aa[2]; END;
IF has(aa[3],y) THEN BB:=coeff(aa[3],y); Fq:=BB[SIZE(BB)]; ELSE Fq:=aa[3]; END;
bb:=SIGN(lcoeff(ffx,x))*(coeff(ffx,y)); IFTE( SIZE(bb)==2,bb={0,bb[1],bb[2]},0);
Cq:=bb[1];
IF has(bb[2],x) THEN BB:=coeff(bb[2],x); Eqq:=BB[2]; ELSE Eqq:=bb[2]; END;
IF has(bb[3],x) THEN BB:=coeff(bb[3],x); Fq:=BB[SIZE(BB)]; ELSE Fq:=bb[3]; END;
ffx:=simplify((Aq*(x)^2)+(Bq*(x)*(y))+(Cq*(y)^2)+Dq*(x)+Eqq*(y)+Fq); 
END;
IF Bq==0 THEN M26:=Conics(ffx); ELSE return (M25); END;
LLLL:=M25[2,2]; LLLLL:=M25[3,2]; LL:=M26[4,2];
IF BBq<>0 THEN
M26[2,2]:={{AA*LLLL[1,1]-AAA*LLLL[1,2],AAA*LLLL[1,1]+AA*LLLL[1,2]},{AA*LLLL[2,1]-AAA*LLLL[2,2],AAA*LLLL[2,1]+AA*LLLL[2,2]}};
IF BBq<0 THEN
M26[3,2]:={{AA*LLLLL[1,1]-AAA*LLLLL[1,2],AAA*LLLLL[1,1]+AA*LLLLL[1,2]},{AA*LLLLL[2,1]-AAA*LLLLL[2,2],AAA*LLLLL[2,1]+AA*LLLLL[2,2]},{AA*LLLLL[3,1]-AAA*LLLLL[3,2],AAA*LLLLL[3,1]+AA*LLLLL[3,2]},{AA*LLLLL[4,1]-AAA*LLLLL[4,2],AAA*LLLLL[4,1]+AA*LLLLL[4,2]} };
END;
IF BBq>0 THEN
M26[3,2]:={{AA*LLLLL[1,1]-AAA*LLLLL[1,2],AAA*LLLLL[1,1]+AA*LLLLL[1,2]},{AA*LLLLL[2,1]-AAA*LLLLL[2,2],AAA*LLLLL[2,1]+AA*LLLLL[2,2]} };
END;
M26[4,2]:={AA*LL[1]-AAA*LL[2],AAA*LL[1]+AA*LL[2]};
IF rowDim(M26)==5 THEN M26[5,2]:=normal(Sxh); END;
END;

IF BBq==0 THEN
M26[2,2]:={AA*LLLL[1]-AAA*LLLL[2],AAA*LLLL[1]+AA*LLLL[2]};
M26[3,2]:={AA*LLLLL[1]-AAA*LLLLL[2],AAA*LLLLL[1]+AA*LLLLL[2]};
M26[4,2]:=Dxh;
M26[5,2]:=Sxh;
END;
M25:=ADDROW(M25,["eqRθ",{simplify(exact(ffx)),r*1_(rad)}],rowDim(M25)+1);
return (M26);
END;


MM:=DET([[Aq, Bq/2, Dq/2],[Bq/2, Cq, Eqq/2],[Dq/2, Eqq/2, Fq]]);
BB:=sqrt(((Aq-Cq)^2)+(Bq^2)); g:=(sqrt((2*BB)/(Aq+Cq+BB))); 
IF MM>0 THEN g:=-1; g:=(sqrt((2*BB)/(g*(Aq+Cq)+BB))); END;
IF MM<0 THEN g:=1;  g:=(sqrt((2*BB)/(g*(Aq+Cq)+BB))); END;
IF MM==0 AND ((Dq^2 +Eqq^2 > ((Aq+Cq)*Fq)) OR (Dq^2 +Eqq^2 < ((Aq+Cq)*Fq)) OR (Dq^2 +Eqq^2 == ((Aq+Cq)*Fq)))  THEN return ("line"); END; 
BBq:=(Bq^2)-4*Aq*Cq;

CASE
IF BBq<0 THEN IF (SIZE(aa)==3 AND SIZE(bb)==3) AND (Aq==Cq) AND Bq==0 THEN 
MSGBOX("circle"); r:=(sqrt(((LL[1])^2+(LL[2])^2)-Fq)); 
return (M25:=[["Radius", r ],["Ctr",LL ]]);
ELSE MSGBOX("ellipse"); 
c:=(sqrt(ABS(LLL[2]-LLL[1]))); 
IF (RRR[1])<(RRR[2]) THEN LLLL:={{LL[1],LL[2]+c},{LL[1],LL[2]-c}}; LLLLL:={{LL[1],LL[2]+sqrt(LLL[2])},{LL[1],LL[2]-sqrt(LLL[2])},{LL[1]+sqrt(LLL[1]),LL[2]},{LL[1]-sqrt(LLL[1]),LL[2]}};
ELSE
LLLL:={{LL[1]+c,LL[2]},{LL[1]-c,LL[2]}}; LLLLL:={{LL[1]+sqrt(LLL[1]),LL[2]},{LL[1]-sqrt(LLL[1]),LL[2]},{LL[1],LL[2]+sqrt(LLL[2])},{LL[1],LL[2]-sqrt(LLL[2])}}; END;
END;
M25:=[["Ecc", g],["Fcs", LLLL ], ["Vtx",LLLLL ],["Ctr", LL ]];
END;

IF BBq>0 THEN 
IF (Aq == -Cq) OR (-Aq == Cq) THEN MSGBOX("rectangular hyperbola"); ELSE MSGBOX("hyperbola"); END;
c:=sqrt(ABS(LLL[2])+ABS(LLL[1]));
IF (RRR[1]) < (RRR[2]) THEN
LLLL:={{LL[1],LL[2]+c},{LL[1],LL[2]-c}}; LLLLL:={{LL[1],LL[2]+sqrt(LLL[2])},{LL[1],LL[2]-sqrt(LLL[2])}}; 
Sx:={expr(cat("(",LL[2]+sqrt(ABS(LLL[2]/LLL[1]))*(x-LL[1]),")","-y")), expr(cat("(",LL[2]-sqrt(ABS(LLL[2]/LLL[1]))*(x-LL[1]),")","-y"))};
Aq:=sqrt(ABS(LLL[2]/LLL[1])); Bq:=1; Cq:=(LL[2]-(sqrt(ABS(LLL[2]/LLL[1]))*LL[1]))/sqrt(Aq^2+Bq^2); Dq:=(LL[2]+(sqrt(ABS(LLL[2]/LLL[1]))*LL[1]))/sqrt(Aq^2+Bq^2);
IFTE(Aq<>0, r1={r-ATAN(Bq/Aq),r+ATAN(Bq/Aq)}, r1={r+(pi/2),r+(pi/2)});
Sxh:={x*COS(r1[1])+y*SIN(r1[1])+Cq,x*COS(r1[2])+y*SIN(r1[2])-Dq};
ELSE
LLLL:={{LL[1]+c,LL[2]},{LL[1]-c,LL[2]}}; LLLLL:={{LL[1]+sqrt(LLL[1]),LL[2]},{LL[1]-sqrt(LLL[1]),LL[2]}};
Sx:={expr(cat("(",LL[2]+sqrt(ABS(LLL[2]/LLL[1]))*(x-LL[1]),")","-y")), expr(cat("(",LL[2]-sqrt(ABS(LLL[2]/LLL[1]))*(x-LL[1]),")","-y"))};
Aq:=sqrt(ABS(LLL[2]/LLL[1])); Bq:=1; Cq:=(LL[2]-(sqrt(ABS(LLL[2]/LLL[1]))*LL[1]))/sqrt(Aq^2+Bq^2); Dq:=(LL[2]+(sqrt(ABS(LLL[2]/LLL[1]))*LL[1]))/sqrt(Aq^2+Bq^2);
IFTE(Aq<>0, r1={r-ATAN(Bq/Aq),r+ATAN(Bq/Aq)}, r1={r+(pi/2),r+(pi/2)});
Sxh:={x*COS(r1[1])+y*SIN(r1[1])+Cq,x*COS(r1[2])+y*SIN(r1[2])-Dq};
END;
LLLLL:=sub(LLLLL,1,2);
M25:=[["Ecc", g],["Fcs", LLLL ], ["Vtx",LLLLL ],["Ctr", LL ], ["Asy", Sx ]];
END;

IF BBq==0 THEN MSGBOX("parabola"); g:=1;
IF aa[1]<>0 THEN bb:={bb[2],bb[3]}; aa:=SIGN(aa[1])*aa;  bb:=(SIGN(bb[1])*bb); BB:=bb[1]; bb:=bb/bb[1];
LLLLL[1]:=(-aa[2]/(aa[1]*2)); LLLLL[2]:=subst(-bb[2],x,LLLLL[1]);  LLLL:={LLLLL[1],LLLLL[2]+(SIGN(lcoeff(-bb[2])))*((BB/4/aa[1]))};
Dx:=expr(cat("(",LLLLL[2]-(SIGN(lcoeff(-bb[2])))*(BB/4/aa[1]),")","-y")); Sx:=expr(cat(LLLLL[1],"-x")); LRa:=BB/aa[1];
Dxh:=(x*COS(r+(pi/2))+y*SIN(r+(pi/2))-(LLLLL[2]-(SIGN(lcoeff(-bb[2])))*(BB/4/aa[1]))); Sxh:=(x*COS(r)+y*SIN(r))-LLLLL[1];
ELSE
aa:={aa[2],aa[3]}; bb:=SIGN(bb[1])*bb; aa:=(SIGN(aa[1])*aa); BB:=aa[1]; aa:=aa/aa[1]; 
LLLLL[2]:=-bb[2]/(bb[1]*2); LLLLL[1]:=subst(-aa[2],y,LLLLL[2]);; LLLL:={LLLLL[1]+(SIGN(lcoeff(-aa[2],y)))*(BB/4/bb[1]),LLLLL[2]};
Dx:=expr(cat("(",LLLLL[1]-(SIGN(lcoeff(-aa[2],y)))*(BB/4/bb[1]),")","-x")); Sx:=expr(cat(LLLLL[2],"-y")); LRa:=BB/bb[1];
Dxh:=(x*COS(r)+y*SIN(r))-(LLLLL[1]-(SIGN(lcoeff(-aa[2],y)))*(BB/4/bb[1])); Sxh:=(x*COS(r+(pi/2))+y*SIN(r+(pi/2)))-LLLLL[2];
END;
M25:=[["Ecc", g],["Fcs", LLLL ], ["Vtx",LLLLL ],["Dtx", Dx ], ["Sym", Sx ], ["LRe", LRa]];
END;
END;
M25;
END;
#end
Isolate line 50:

IF has(aa[3],y) THEN BB:=coeff(aa[3],y); Fq:=BB[SIZE(BB)]; ELSE Fq:=aa[3]; END;

That is where it crashes on my tests.

I have had trouble with the coeff() command. I think this is where your trouble might be. It seems to be quirky, and there may have been some changes since this:

http://www.hpmuseum.org/forum/newreply.php?tid=10200

-Dale-
That's because aa has only two elements, has(aa[3],y) should stop the program and error, but aa[3] is inside has() and since there is no exceptions on the Prime, an error string object is passed to has and this causes the crash. I have commited a fix for the crash, but your code will still have to handle the size of aa (for example you might insert a 0 at the begnning if size(aa)==2).
Thanks guys....the above code has the fix;

why CAS and HOME application interaction parse [[]] as {{}} in CAS mode;?
and sometimes it does NOT execute;
http://www.hpmuseum.org/forum/thread-10107.html

example;
vget() in CAS
The updated/Complete above code calculates Parameters of a Conic Section; M25 and M26


Ecc= Eccentricity
Fcs=Foci
Vtx=Vertex
Dtx=Directrix expression
Sym=Symmetry expression
Asy=Asymptotes expressions
Ctr=Center
LRe=Latus Rectum Length
Radius= Radius
eqRθ = Simplified Canonical Expression; Rotating Angle
Reference URL's