Post Reply 
Help Match?
12-01-2022, 03:04 PM
Post: #3
RE: Help Match?
(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;

Tom L
Cui bono?
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
Help Match? - toml_12953 - 12-01-2022, 01:33 AM
RE: Help Match? - roadrunner - 12-01-2022, 02:11 AM
RE: Help Match? - toml_12953 - 12-01-2022 03:04 PM



User(s) browsing this thread: 1 Guest(s)