Sketcher: Ellipse support

[marktaff]

Yep, conic sections are a *lot* easier in polar; in fact, they *all* have the same equation! I have an advantage in that I have my orbital mechanics text (Curtiss) sitting right here.

The equation you posted is one form of the 'orbit equation', aka Kepler's First Law. The e is the eccentricity; if it is 0, you have a circle; if it is 1, you have a parabola; if it is > 1, you have a hyperbola; if it is between 0 and 1, you have an ellipse. This means that if we approach the math from a polar perspective, we get solutions for all the conic sections for free (neglecting issues of what snaps, constraints, etc we want).

The 'perifocal' frame is a reference frame centered at the focus (F), where the other focus would be F' (read as "F prime"). The point on the apse line (defined in R^2 as the line defined by the two foci) closest to F is periapsis, and r_p is the position of that point relative to F; r_p is the angle (theta) reference, i.e. 0 degrees (radians), with positive direction being counter-clockwise. N.B. 'r' is the notation we use for 'position', so r(theta) is the position of the point on the conic section that a line from F to the ellipse with the angle of theta (relative to our angle reference). In polar coords, position would be a two-dimensional vector of <r, theta>, which can be converted to a two dimensional vector in <x, y> as required.

Any point or equation can be translated back and forth between coord systems; my approach would be to use whatever coord system was easiest for the task at hand. I will say that I would find constraint(s) on F to be very useful (the focus of my parabolic reflector must be 'here' relative to this feature).

For example, we could use a perifocal frame to do everything we need to with a conic section, and then convert our results back to cartesian coords if we needed to do something that is easier in cartesian coords.

Anyhow, like I said, I still have no knowledge of the FC codebase yet, I was just wondering about the approach being taken with the math here.

[DeepSOIC]

I think we can nail the internal representation of an ellipse now.
The abovementioned polar formula does not stand for parabola case (because a is infinite for parabola, while (1-e^2) is zero, thus the product in nominator is uncertain). So a(vector) and e are not good enough.

ellipse in perifocal coordinates

350px-Ellipse_Polar_svg.png (18.23 KiB) Viewed 1428 times

So, my suggestion is to internally represent ellipse, as well as any conic by:
1. vector to a focus point F1, which defines the position of the conic as a whole
2. vector from F1 to point A (as denoted by d1 on picture, but a vector value). The vector portion of this value sets the orientation of the ellipse, while its length, coupled with...
3. ...the eccentricity e (any positive value) defines the size and shape of the ellipse,
From the user's point of view, this is not acceptable (i.e. if we expose only points A and F1 - no access to the center of ellipse), so UI stuff needs to be tailored separately for each conic case. But at least this approach will hopefully keep the math consistent and unified.
What do you think?

BTW. Shame, but I have never looked into sketcher code. I just help with math here, but for some stupid reason I also insist on ways to implement stuff

[abdullah]

Here is a small patch to fix the constraint incompatibility to the master branch.

Until I saw that diff I could not understand what it was... of course it is the numbering of the constraints, so that the same as in master is kept... I was not giving this much attention, but it is true, I am integrating this in my branch. Additionally I will be numbering them (forcing them to take the right number), so that it is clear that order matters.

[abdullah]

I have a more human-readable way to tell this: addConstraintXXX is supposed to remove exactly one degree of freedom. UI constraint "coincident" has to remove two degrees of freedom, so has two addConstraintXXX calls. Correct?

Conclustion: a tangent constraint must be a single new row into the matrix, since it removes only one degree of freedom.
Or no? More rows can be added if we introduce additional degrees of freedom.

It is correct!! The rank of the gradient matrix is the amount of degrees of freedom removed by the constraints

[abdullah]

Numeric derivatives should be avoided whenever possible. And it is possible to get code that is much more efficient. A quartic equation has algebraic solutions you just have to choose the right one.

I agree.

I have actually done some computations these last couple of days to visualize the solutions. Assuming A is the major axis the following code produce the image of the theta-angle (named t by abdullah I think) of the closest point on the ellipse.

I will look at the code later. Just a clarification t is the "eccentric anomaly" as in: http://en.wikipedia.org/wiki/Eccentric_anomaly. I have not used theta because IMO it induces to error with the elevation angle of the point of the ellipse, which is not. Also in wikipedia, t was used.

Getting the partials from this code is quite straightforward and doesn't require 600kb of code.

Please, folks, stop referring to the 600 kb. That was only a test in full cartesian. When I wrote that, it was only to illustrate how "complex" this can get in full cartesian non-parametric. No part of the code of any branch of my repository has a 600 kb partial, and no gcc compilers were killed in the tests

BTW, I selected the "right" root in that cartesian code by calculating the distance to the points on ellipse and selecting the closest to the midpoint of the line.

Thanks for you code!!, I will look into it later on.

[abdullah]

To avoid the solution jumping from one point on the ellipse and circle to another point on the ellipse and circle I think some kind of extra data is required. An intermediate line is one possibility where you constrain one end of the line to the ellipse and circle curves and add tangency constraints ellipse-line and circle-line.

This tangency point does not necessarily have to be visible, but it can be quite useful if it is.

This is brilliant! Here's how I see it:
A user selects two curves and a point and hits "tangent" button (if no point is supplied, one is created).
The point is constrained onto two objects (two point-on-object constraints, already done). And the final constraint uses the point in calculation of a new, much simpler error function, assuming that the point is the point of tangency. I can take care for making such function, for line, circle and possibly even second ellipse, I feel it's easy! It's just calculating tangent vectors to the curves at the point and F=length of vector product of those two vectors.
If the user deletes the point, it will automatically kill the constraint since it will be associated with the constraint.
What do you think?

I think it is a possibility. If we do not manage to get better, it may be the solution. If we manage to get better, it is an overcomplicated solution...

[abdullah]

Any point or equation can be translated back and forth between coord systems; my approach would be to use whatever coord system was easiest for the task at hand.

It makes a lot of sense.

I will say that I would find constraint(s) on F to be very useful (the focus of my parabolic reflector must be 'here' relative to this feature).

We will have a way to represent the foci in the UI, like, for sure. It won't be basic subelement of the ellipse (like the center), but it will be possible to constrain one or both foci, as needed.

Anyhow, like I said, I still have no knowledge of the FC codebase yet, I was just wondering about the approach being taken with the math here.

Probably having no knowledge of the FC codebase is at this stage an advantage for this discussion, as you are not limited by possible coding constraints. I very much appreciate you introducing us to the perifocal equations. I am quite sure it will change our approach to the maths. Your comments are very much welcome.

[abdullah]

So, my suggestion is to internally represent ellipse, as well as any conic by:
1. vector to a focus point F1, which defines the position of the conic as a whole
2. vector from F1 to point A (as denoted by d1 on picture, but a vector value). The vector portion of this value sets the orientation of the ellipse, while its length, coupled with...
3. ...the eccentricity e (any positive value) defines the size and shape of the ellipse,

So your proposal:
1. So more than a vector, a point having coordinates in cartesian F1_X and F1_Y. (2 DOFs)
2. So a non-unitary vector defining a direction and distance d1 (2 DOFs)
3. Eccentricity (1 DOF)

Sounds plausible to me.

From the user's point of view, this is not acceptable (i.e. if we expose only points A and F1 - no access to the center of ellipse), so UI stuff needs to be tailored separately for each conic case. But at least this approach will hopefully keep the math consistent and unified.
What do you think?

UI staff is a totally separate topic. One way or the other, at the end, the user will have accessible anything that can reasonably be expected. Right now I am most insterested in a solver representation that does not produce convergence problems and that is stable (does not jump to the sides when dragging it). It seems that a source of those problems is when one of the parameters is an angle (I have seen that at the end of the convergence it is not strange to arrive to things like 1000 radians). So it keeps rolling while converging. This is IMHO unwanted behaviour and should be out (and I thank again Ulrich for pointing this out).

BTW. Shame, but I have never looked into sketcher code. I just help with math here, but for some stupid reason I also insist on ways to implement stuff

[/quote]

That is very much ok. If I think it is not possible, I will tell you. I may be right or just show my ignorance, in the latter case either you or others will make me aware of my mistake. So we all win...

[abdullah]

FUTURE STEPS IN ELLIPSE IMPLEMENTATION

I thank all of you for such a fruitful discussion. I am confident that now we are on the right track to have ellipse working soon.

IMO, the main issue now is the solver representation.

1. Solver representation

I think this is the major decision we are taking. The solver representation has to:
1. Offer fast convergence (as opposed to phi parameter).
2. Offer stability.
3. Be convenient when drafting the error functions and partials (allow for simpler representation).

What do we go for as solver representation, Ulrich-DeepSOIC or Mark-DeepSOIC ones?

I would really like to hear any strong opinions or possible drawbacks. To my untrained eye, they look very similar.

Should we fear any (convergence) problem with the fact that e is between 0 and 1 for the ellipse? I think probably not, because it is not periodical as phi is, but I would like to hear any comment you have to this respect.

With your comments I expect that we can decide on this topic soon. Once the solver parameters are fixed...

...I will write a reference framework branch on my github with this solver representation in branch "ellipse_newsolver_master".

I invite anyone here to join to github, and I will give you access to my FreeCAD repository, so that you can create you own branches for your own tests if you want; or you can just fork my ellipse repository and ask me for a pull request if you want to contribute something to "ellipse_newsolver_master". My idea is that the implementation goes growing this way.

Of course, you can help me by throwing equations to this thread, if you prefer to do so. I am very much ok with that.

2. Solver constraint error functions and partials

We need to review the work done so far in the frame of the new solver parameters. It is worth mentioning that while the error function usually can merely be rewritten, the partials will change, when changing the solver parameters. It might also come a situation where with the new representation, an alternative approach is deemed more convenient for the error function.

I propose we continue following Mark's idea of trying to put constraint equations in the easiest form for each problem.

I take any kind of representation of the equations you can offer (as far as I can understand the syntax, if I do not, I ask), python, latex, maxima... I will anyway do the calculations of the partials with sage and post the results so that you can check it, improve them, correct them,...

3. UI

I will work on the interface based on solution 4, as presented before. I will introduce a new UI constraint for general use for any element that wants to implement it called "Align to internal geometry constraint". It will be responsible for fully constraining an element within another element (e.g. a point on a focus, a line on a major radius,...).

4. Ellipse arcs

5. Testing before making it into master

Please, folks, really really really tell me now if you think anything I wrote here is not right or should be further discussed...

[DevJohan]

Ellipse-ellipse tangent

ellipsetoellipse.png (8.62 KiB) Viewed 1371 times