C++

Forum

 
"Completing the the square"

everybodyCallsMeM (16)
Given the general equation for an ellipse

Ax^2+ Bxy + Cy^2 + Dx + Ey + F = 0

How would I code the equation into simplified form to
extract the central point (h,k) and the radius (r)?

so I want to get from the general form to this

(x – h)2 + (y – k)2 = r^2

with the center being at the point (h, k) and the radius being r

I'm having difficulties with getting to this point
Ganado (6814)
The general equation for an ellipse allows for oblique ellipses.

It sounds like you're trying to express the equation of a circle in terms of the general equation for an ellipse?

I wouldn't be able to explain it better than Wikipedia, so let me just link the section that talks about this.
https://en.wikipedia.org/wiki/Ellipse#General_ellipse

You can see here it gives the letters A, B, C, D, E, F in terms of "a" (semi-major axis), "b" (semi-minor axis), "(x_c, y_c)"* (center coordinate), and "theta" the rotation angle.

* in your case: (h, k)

If you're trying to simplify this for a circle,
1. set a = b = r.
2. set theta = 0

So, you get something along the lines of:
A = r^2
B = 0
C = r^2
D = -2 A h
E = - 2 C k
F = A h^2 + C k^2 - 2 r^2
F = A h^2 + C k^2 - r^4

This probably can be further simplified once you expand and combine the terms.

Hoping I didn't make a typo...
r^2 x^2 + r^2 y^2 - 2r^2hx - 2r^2ky + r^2h^2 + r^2k^2 - 2r^2 = 0
r^2 ( x^2 + y^2 - 2hx - 2ky + h^2 + k^2 - 2 ) = 0
...that -2 looks suspicious, however.

r^2 x^2 + r^2 y^2 - 2r^2hx - 2r^2ky + r^2h^2 + r^2k^2 - 2r^2 = 0
r^2 ( x^2 + y^2 - 2hx - 2ky + h^2 + k^2 - r^2 ) = 0
x^2 + y^2 - 2hx - 2ky + h^2 + k^2 - r^2 = 0
x^2 + y^2 - 2hx - 2ky + h^2 + k^2 = r^2


Edit: Fixed typo.
Last edited on
everybodyCallsMeM (16)
thanks
that looks do able!

yeah my bad, but a circle is an ellipse just with the special case when the semi major and semi minor are equal.
The general form of the equation for both is Ax^2 + Bxy + Cy2 + Dx + Ey + F =0.

For the circle A and C must be equal and B must equal 0.
For the ellipse, A and C cannot be equal but must have the same sign. B must equal 0
Ganado (6814)
By the way, I realized what I did wrong. In (F), I "simplified" a^2 * b^2 to be 2 * r^2, when it should be r^2 * r^2, or r^4.
Last edited on
lastchance (6980)
It's a bit sensitive to round-off errors and isn't fully checked ...

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
 
#include <iostream>
#include <cmath>
using namespace std;

const double SMALL = 1.0e-10;
const double PI = 4.0 * atan( 1.0 );

//======================================================================

bool nearZero( double x )
{
   return abs( x ) < SMALL;
}

//======================================================================

// Rotate position vector by theta radians
void rotate( double x, double y, double theta, double &xp, double &yp )
{
   double c = cos( theta );
   double s = sin( theta );
   xp = c * x - s * y;
   yp = s * x + c * y;
}

//======================================================================

// Coefficients in orthonormal axes rotated by theta radians (equivalent to rotating coordinates by -theta )
void rotateCoefficients( double A, double B, double C, double D, double E, double F, double theta,
                         double &AP, double &BP, double &CP, double &DP, double &EP, double &FP )
{
   double c = cos( theta ), csq = c * c;
   double s = sin( theta ), ssq = s * s;
   double sc = s * c;

   AP = csq * A + sc * B + ssq * C;
   BP = ( csq - ssq ) * B + 2 * sc * ( C - A );
   CP = ssq * A - sc * B + csq * C;
   DP = c * D + s * E;
   EP = c * E - s * D;
   FP = F;
}

//======================================================================

int main()
{
   double A, B, C, D, E, F;
   cout << "Conic section Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0\n";
   cout << "Input A, B, C, D, E, F: ";   cin >> A >> B >> C >> D >> E >> F;

   double theta;
   double rhs, xc, yc;

   // Determine angle by which to rotate axes if B != 0
   if ( nearZero( B ) )
   {
      theta = 0.0;
   }
   else if ( nearZero( A - C ) )
   {
      theta = PI / 4;
   }
   else
   {
      theta = 0.5 * atan( B / ( A - C ) );
   }
   // Get coefficients in rotated axes
   if ( !nearZero( B ) ) rotateCoefficients( A, B, C, D, E, F, theta, A, B, C, D, E, F );
// cout << A << ' ' << B << ' ' << C << ' ' << D << ' ' << E << ' ' << F << '\n';
   cout << "Axes rotation: " << theta * 180.0 / PI << " degrees\n";


   // Pick off cases
   if ( nearZero( A ) && nearZero( C ) )
   {
      cout << "Not a conic section\n";
   }

   else if ( nearZero( A ) || nearZero( C ) )
   {
      if ( nearZero( A ) ) rhs = 0.25 * E * E / ( C * C ) - F / C;
      else                 rhs = 0.25 * D * D / ( A * A ) - F / A;
      if ( rhs <= 0.0 )
      {
         cout << "Not a conic section\n";
      }
      else
      {
         cout << "Parabola\n";
      }
   }

   else if ( A * C < 0 )
   {
      rhs = 0.25 * D * D / A + 0.25 * E * E / C - F;
      {
         if ( nearZero( rhs ) )
         {
            cout << "Line pair\n";
         }
         else
         {
            cout << "Hyperbola\n";
         }
      }
   }
   else
   {
      rhs = 0.25 * D * D / A + 0.25 * E * E / C - F;
      {
         if ( rhs <= 0 )
         {
            cout << "Not a conic section\n";
         }
         else
         {
            xc = -0.5 * D / A;
            yc = -0.5 * E / C;
            rotate( xc, yc, theta, xc, yc );
            if ( nearZero( A - C ) )
            {
               cout << "Circle:\n";
               cout << "   centre: " << xc << ", " << yc << '\n';
               cout << "   radius: " << sqrt( rhs / A ) << '\n';
            }
            else
            {
               cout << "Ellipse:\n";
               cout << "   centre: " << xc << ", " << yc << '\n';
               cout << "   semi-axes: " << sqrt( rhs / A ) << " and " << sqrt( rhs / C ) << '\n';
            }
         }
      }
   }
}


Conic section Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
Input A, B, C, D, E, F: 0 2 0 0 0 -4
Axes rotation: 45 degrees
Hyperbola


Conic section Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
Input A, B, C, D, E, F: 9 0 4 -36 8 4
Axes rotation: 0 degrees
Ellipse:
   centre: 2, -1
   semi-axes: 2 and 3


Conic section Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
Input A, B, C, D, E, F: 4 0 4 -16 8 4
Axes rotation: 0 degrees
Circle:
   centre: 2, -1
   radius: 2


Conic section Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
Input A, B, C, D, E, F: 2 -1 -1 0 0 0
Axes rotation: -9.21747 degrees
Line pair

Last edited on
everybodyCallsMeM (16)
TEST INPUT. TEST OUTPUT
1. 1, 0, 1, -2, -6, -26. 1. Circle, (1, 3), 6
2. 9, 0, 16, -72, 64, 64 2. Ellipse (4, -2), 8
3. 16, 0, 16, 0, 0, -64 3. Circle (0, 0), 2
4. 9, 0, 4, 0, 0, -144 4. Ellipse (0, 0), 12
Where the test output for the elliptical geometry is central coordinates and the length of the semi major axis. I'm attempting to find what is going on with the arithmetic because the output I am getting is always off by half of the calculated value.
MikeStgt (466)
the output I am getting
... by using which program? Yours or this one: http://www.cplusplus.com/forum/beginner/252180/#msg1110269
And, what do you mean by "off half of the calculated value"? 50% less, or 100% more? Based on what? How about "the semi major axis" (or major semi axis?) Is this the same semi as in semi-equivalent or semi-conservative or semicolon?
everybodyCallsMeM (16)
My apologies for the vague response

I had essentially stared over with the code you have referenced to closely mirror yours

to which I am thankful for, I truly am.


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
 
#include <iostream>
#include <cmath>

using namespace std;

const double SMALL = 1.0e-10;
const double PI = 4.0 * atan(1.0);

//======================================================================

bool nearZero(double x)
{
    return abs(x) < SMALL;
}

//======================================================================

// Rotate position vector by theta radians
void rotate(double x, double y, double theta, double &xp, double &yp)
{
    double c = cos(theta);
    double s = sin(theta);
    xp = c * x - s * y;
    yp = s * x + c * y;
}

//======================================================================

// Coefficients in orthonormal axes rotated by theta radians (equivalent to rotating coordinates by -theta )
void rotateCoefs(double A, double B, double C, double D, double E, double F, double theta,
    double &AP, double &BP, double &CP, double &DP, double &EP, double &FP)
{
    double c = cos(theta), csq = c * c;
    double s = sin(theta), ssq = s * s;
    double sc = s * c;

    AP = csq * A + sc * B + ssq * C;
    BP = (csq - ssq) * B + 2 * sc * (C - A);
    CP = ssq * A - sc * B + csq * C;
    DP = c * D + s * E;
    EP = c * E - s * D;
    FP = F;
}

//======================================================================

int main()
{

    double theta;
    double rhs, xc, yc;
    double semiMajor, semiMinor;


    // Create an array structure for a six entry array to fill the general eqn.
    struct arrayInputs
    {
        double A, B, C, D, E, F;
    }coef[5];

    // Asking for and inputing values for the coefficents in of the 5 eqns. 
    for (int i = 1; i < 6; ++i)
    {
        cout << "Enter the coefficients for equation number " << i << ": ";
        cin >> coef[i].A >> coef[i].B >> coef[i].C >> coef[i].D >> coef[i].E >> coef[i].F;
        cout << "\n";

        while (! coef[i].B == 0)  // restricting the value for coef[i].B to have value 0
        {
            cout << "The value for the second coefficient in the needs to be zero! \n";
            cout << "Enter the value 0";
            cin >> coef[i].B;
        }

        cout << "General Equation number " << i << " is: " << coef[i].A << "x^2 + "
            << coef[i].B << "xy + " << coef[i].C << "y^2 + " << coef[i].D << "x + " << coef[i].E
            << "y + " << coef[i].F << " = 0 " << endl;

    // Determine angle by which to rotate axes if B != 0
    if (nearZero(coef[i].B))
    {
        theta = 0.0;
    }
    else if (nearZero(coef[i].A - coef[i].C))
    {
        theta = PI / 4;
    }
    else
    {
        theta = 0.5 * atan(coef[i].B / (coef[i].A - coef[i].C));
    }
    // Get coefficients in rotated axes
    if (!nearZero(coef[i].B)) rotateCoefs(coef[i].A, coef[i].B, coef[i].C, coef[i].D, coef[i].E, coef[i].F, theta, coef[i].A, coef[i].B, coef[i].C, coef[i].D, coef[i].E, coef[i].F);
    // cout << A << ' ' << B << ' ' << C << ' ' << D << ' ' << E << ' ' << F << '\n';
    cout << "Axes rotation: " << theta * 180.0 / PI << " degrees\n";


    // Pick off cases
    if (nearZero(coef[i].A) && nearZero(coef[i].C))
    {
        cout << "Not a conic section\n";
    }

    else if (nearZero(coef[i].A) || nearZero(coef[i].C))
    {
        if (nearZero(coef[i].A)) rhs = 0.25 * coef[i].E * coef[i].E / (coef[i].C * coef[i].C) - coef[i].F / coef[i].C;
        else                 rhs = 0.25 * coef[i].D * coef[i].D / (coef[i].A * coef[i].A) - coef[i].F / coef[i].A;
        if (rhs <= 0.0)
        {
            cout << "Not a conic section\n";
        }
        else
        {
            cout << "Parabola\n";
        }
    }

    else if (coef[i].A * coef[i].C < 0)
    {
        rhs = 0.25 * coef[i].D * coef[i].D / coef[i].A + 0.25 * coef[i].E * coef[i].E / coef[i].C - coef[i].F;
        {
            if (nearZero(rhs))
            {
                cout << "Line pair\n";
            }
            else
            {
                cout << "Hyperbola\n";
            }
        }
    }
    else
    {
        rhs = 0.25 * coef[i].D * coef[i].D / coef[i].A + 0.25 * coef[i].E * coef[i].E / coef[i].C - coef[i].F;
        {
            if (rhs <= 0)
            {
                cout << "Not a conic section\n";
            }
            else
            {
                xc = -0.5 * coef[i].D / coef[i].A;
                yc = -0.5 * coef[i].E / coef[i].C;
                rotate(xc, yc, theta, xc, yc);
                if (nearZero(coef[i].A - coef[i].C))
                {
                    cout << "Circle:\n";
                    cout << "   central point: " << "(" << xc << ", " << yc << ")" << '\n';
                    cout << "   radius: " << sqrt(rhs / coef[i].A) << '\n';
                }
                else
                {
                    cout << "Ellipse:\n";
                    cout << "   central point: " << xc << ", " << yc << '\n';
                    
                    // determine the semi major and the semi minor axes 
                    if ((sqrt(rhs / coef[i].A)) > (sqrt(rhs / coef[i].C)))
                    {
                        semiMajor = (sqrt(rhs / coef[i].A));
                        semiMinor = (sqrt(rhs / coef[i].C));
                    }
                    else
                    {
                        semiMajor = (sqrt(rhs / coef[i].C));
                        semiMinor = (sqrt(rhs / coef[i].A));
                    }

                    cout << "   a = " << semiMajor << ", b = " << semiMinor << '\n';
                }
            }
        }
    }

    }
    return 0;
}


with the sample input of: 1, 0, 4, -6, -16, -11

This is the output from the console window:

Enter the coefficients for equation number 2: 1 0 4 -6 -16 -11

General Equation number 2 is: 1x^2 + 0xy + 4y^2 + -6x + -16y + -11 = 0
Axes rotation: 0 degrees
Ellipse:
central point: 3, 2
a = 6, b = 3

When the expected sample output calculates the semi major axis to be 12

I greatly appreciate the help.
MikeStgt (466)
Thank you to publish your latest revision.
I run it and the 'role model' too:
Enter the coefficients for equation number 1: 9 0 4 -36 8 4

General Equation number 1 is: 9x^2 + 0xy + 4y^2 + -36x + 8y + 4 = 0 
Axes rotation: 0 degrees
Ellipse:
   central point: 2, -1
   a = 3, b = 2

and (drum roll...)
Conic section Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
Input A, B, C, D, E, F: 9 0 4 -36 8 4
Axes rotation: 0 degrees
Ellipse:
   centre: 2, -1
   semi-axes: 2 and 3


Once more, another example...
Enter the coefficients for equation number 1: 1 0 4 -6 -16 -11

General Equation number 1 is: 1x^2 + 0xy + 4y^2 + -6x + -16y + -11 = 0 
Axes rotation: 0 degrees
Ellipse:
   central point: 3, 2
   a = 6, b = 3

and...
Conic section Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
Input A, B, C, D, E, F: 1 0 4 -6 -16 -11
Axes rotation: 0 degrees
Ellipse:
   centre: 3, 2
   semi-axes: 6 and 3

So both routines are correct or both are not correct (for the two examples).

When the expected sample output calculates the semi major axis to be 12
Question: major axis or semi major axis? Please make sure.

BTW, if you "sort" the output of semi axis by size, here
1
2
3
4
5
6
7
8
9
10
11
 
                    // determine the semi major and the semi minor axes 
                    if ((sqrt(rhs / coef[i].A)) > (sqrt(rhs / coef[i].C)))
                    {
                        semiMajor = (sqrt(rhs / coef[i].A));
                        semiMinor = (sqrt(rhs / coef[i].C));
                    }
                    else
                    {
                        semiMajor = (sqrt(rhs / coef[i].C));
                        semiMinor = (sqrt(rhs / coef[i].A));
                    }

this is equivalent a rotation by 90°. (In addition, the comparison of coeff_A and coeff_C could be simplified.)
everybodyCallsMeM (16)
thank you. I was looking over everything again, and I feel like a doofus. Everything is correct!

for the most part. some minor tweaks on my part to some of the logic, and increments. And it works for how I would like it to.
Last edited on
MikeStgt (466)
... and I feel like a doofus
Get used to it, it's human to be error-prone. (Disadvantageous, even hurtful for individuals it accelerates evolution.)
Topic archived. No new replies allowed.