(12-01-2022 02:11 AM)roadrunner Wrote: Three mistakes i found are:
line 220 should have a semicolon after continue
line 215 temp variable isn't defined
line 6, variables c and r are not defined
-road
Thanks! I still had problems because all function parameters are passed by value rather than reference but by making the arrays global, I bypassed that problem.
Code:
EP := 1e-15; // approximate zero
coef := makemat(0,6);
r := makemat(0,6);
export conic()
begin
local t;
// A B C D EE F
set_conic( 14, -4, 11, -44, -58, 71 );
put_conic( coef );
t := conic_type( coef );
if t <> 0 then
standard_conic( coef, t );
put_conic( r );
end;
end;
// set values of conic coefficient array
set_conic( a, b, c, d, ee, f )
begin
coef(1) := a;
coef(2) := b;
coef(3) := c;
coef(4) := d;
coef(5) := ee;
coef(6) := f;
end;
// determines the type of conic section and returns a type code
//
// 0 - degenerate
// 1 - circle,
// 2 - ellipse
// 3 - hyperbola
// 4 - parabola
//
conic_type( x )
begin
local str := "conic is ";
local t := 0;
local A := x(1);
local B := x(2);
local C := x(3);
local D := x(4);
local E := x(5);
local F := x(6);
// | F D/2 E/2 |
// M = | D/2 A B/2 |
// | E/2 B/2 C |
local M11 := A*C - B*B/4;
local M23 := D*E/4 - B*F/2;
local det := A*C*F + (B*D*E - B*B*F - C*D*D - A*E*E)/4;
if det == 0 then // conic is degenerate
t := 0;
if M11 > EP then
str := str + "point";
else
if M11 < -EP then
if A + C == 0 then
str := str + "two perpendicular lines";
else
str := str + "two intersecting lines";
end;
else
if M23 == 0 then
str := str + "a single line";
else
str := str + "two parallel lines";
end;
end;
end;
else
if M11 > EP then // conic is ellipse, parabola, or hyperbola
if A*det > 0 then // no solution
t := 0;
str := str + "nothing";
else
if (A == C) and (B == 0) then
str := str + "a circle";
t := 1;
else
str := str + "an ellipse";
t := 2;
end;
end;
else
if M11 < -EP then
str := str + "a hyperbola";
t := 3;
else
str := str + "a parabola";
t := 4;
end;
end;
end;
print(str);
return t;
end;
// find standard conic equation corresponding to general equation
// by removing rotation and then translation
standard_conic( x, t )
begin
local A := x(1);
local B := x(2);
local C := x(3);
local D := x(4);
local E := x(5);
local F := x(6);
// compute rotation for coordinate frame of standard conic
// and coefficients of rotated conic
local theta;
if A == C then
theta := PI/4;
else
theta := 0.5*atan( B / (A - C) );
end;
local Sa := sin(theta);
local Ca := cos(theta);
local Ar := A*Ca*Ca + B*Sa*Ca + C*Sa*Sa;
local Cr := A*Sa*Sa - B*Sa*Ca + C*Ca*Ca;
local Dr := D*Ca + E*Sa;
local Er := -D*Sa + E*Ca;
local Fr := F;
print("rotation = "+ round((180*theta/PI),4)+ "°");
// compute offset of coordinate frame of standard conic
// and coefficients of translated conic
local Xo, Yo, a2, b2, c, p;
local i;
for i := 1 to 6 do // init coefficients to zero
r(i) := 0;
end;
if (t == 1) or (t == 2) or (t == 3) then
// circle, ellipse, or hyperbola
Xo := -Dr/(2*Ar);
Yo := -Er/(2*Cr);
print("translation = {"+ round(Xo,4)+ ", "+ round(Yo,4)+ "}");
r(1) := Ar/(Ar*Xo*Xo + Cr*Yo*Yo - Fr);
r(3) := Cr/(Ar*Xo*Xo + Cr*Yo*Yo - Fr);
r(6) := -1.0;
if abs(r(1)) < abs(r(3)) then // foci on x-axis
a2 := abs(1 / r(1));
b2 := abs(1 / r(3));
if t == 2 then // ellipse
c := sqrt( a2 - b2 );
print("foci = {±"+ round(c,4)+ ", 0.0000}");
print("eccentricity = "+ c / sqrt(a2));
else
if t == 3 then // hyperbola
c := sqrt( a2 + b2 );
print("foci = {±"+ round(c,4)+ ", 0.0000}");
print("eccentricity = "+ c / sqrt(a2));
end;
end;
else // foci on y-axis
a2 := abs(1 / r(3));
b2 := abs(1 / r(1));
if t == 2 then // ellipse
c := sqrt( a2 - b2 );
print("foci = {0.0000, ±"+ round(c,4)+ "}");
print("eccentricity = "+ c / sqrt(a2));
else
if t == 3 then // hyperbola
c := sqrt( a2 + b2 );
print("foci = {0.0000, ±"+ round(c,4)+ "}");
print("eccentricity = "+ c / sqrt(a2));
end;
end;
end;
else
if t == 4 then // parabola
if abs(Ar) > abs(Cr) then // form: x² - 4·p·y = 0
Xo := Dr/(2*Ar);
Yo := (Fr - Dr*Dr/(4*Ar))/Er;
print("translation = {"+ round(Xo,4)+ ", "+ round(Yo,4)+ "}");
r(1) := 1;
r(5) := Er / Ar;
p := -r(5)/4;
print("focus = {0, "+ round(p,4)+ "}");
end;
else // form: y² - 4·p·x = 0
Xo := (Fr - Er*Er/(4*Cr))/Dr;
Yo := Er/(2*Cr);
print("translation = {"+ round(Xo,4)+ ", "+ round(Yo,4)+ "}");
r(3) := 1;
r(4) := Dr / Cr;
p := -r(4)/4;
print("focus = {"+ round(p,4)+ ", 0}");
end;
end;
end;
// display a general conic equation
put_conic( x )
begin
local i;
local s := false;
local temp := "";
for i := 1 to 6 do
if x(i)==0 then
continue;
end;
if x(i) < 0 then
temp := temp + " - ";
else
if (i > 1) and s <> 0 then
temp := temp +" + ";
end;
end;
s := true;
if (abs(x(i)) <> 1) or (i == 6) then
temp := temp + round(abs(x(i)),4);
if i <> 6 then
temp := temp + "·";
end;
end;
case
if i == 1 then temp := temp + "x²" end;
if i == 2 then temp := temp + "x·y" end;
if i == 3 then temp := temp + "y²" end;
if i == 4 then temp := temp + "x" end;
if i == 5 then temp := temp + "y" end;
end;
end;
print(temp + " = 0");
end;
Search
Member List
Calendar
Help
