Sketcher: Ellipse support

Here's the place for discussion related to coding in FreeCAD, C++ or Python. Design, interfaces and structures.
abdullah
Posts: 3584
Joined: Sun May 04, 2014 3:16 pm

Re: Sketcher: Ellipse support

Postby abdullah » Sat Sep 13, 2014 5:09 am

Hey guys!! I have a question for you. Do you know of any other case where the sketcher reports a constraint as redundant when you would not expect it to?
ulrich1a
Posts: 1958
Joined: Sun Jul 07, 2013 12:08 pm

Re: Sketcher: Ellipse support

Postby ulrich1a » Sat Sep 13, 2014 10:45 am

abdullah wrote:Hey guys!! I have a question for you. Do you know of any other case where the sketcher reports a constraint as redundant when you would not expect it to?
In all the cases I investigated so far, I could find a reason for the complaint about redundant constraints. It is more the other way around: In several cases there are redundant constraints, that are not reported.

My hypothesis is, you are using the parameter t together with xp and yp. But this will introduce a redundancy into the equation system by definition. Only two parameters are needed to define a point on the ellipse. The third is more or less redundant. If only xp and yp are used, then there are two different yp-values possible for a given xp. But this is a problem of the ellipse, as it has two valid yp-values for each valid xp.
It may be possible to formulate the equation system to be solved only with xp and t. The needed value yp is just calculated from the actual xp and t. This may have an advantage, but I think a comparison will be needed to decide which formulation is working better.
If my hypothesis is true, only in a few case the code works, due to the weak detection of redundancies in the sketcher solver.
(I have to admit, I never looked into your code. I am not a c++ programmer. The one time, I had a look into the sketcher code, I was totally lost.) Please correct me, if I am wrong.

The ellipse-simulation-sketch does not have the t-parameter. It based on the pin and string method for an ellipse construction. And it works.

By the way: do you know wxMaxima? This program may be helpful in order to generate the needed expressions.

Ulrich
abdullah
Posts: 3584
Joined: Sun May 04, 2014 3:16 pm

Re: Sketcher: Ellipse support

Postby abdullah » Sat Sep 13, 2014 3:44 pm

In all the cases I investigated so far, I could find a reason for the complaint about redundant constraints. It is more the other way around: In several cases there are redundant constraints, that are not reported.
Well I am learning how it works, so maybe in the future we can take a look to those cases.

The redundant check is done in diagnose(). A matrix of the system, number of constraints times number of parameters (x of point, y of point, phi, a, b, ...) is built by iterativelly using the gradients provided by the constraints for each element. Once created the QR descomposition of this matrix is calculated. The rank of the system is calculated and compared to the number of constraints, if lower, then obviously one constraint is redundant (the row is not independent).
My problem is related to the point on ellipse not constraining the y coordinate of the point on ellipse, only the x. Then putting a vertical constraint, which is an equality constraint (x of p1=x of p2), does not affect the y of that point, so no conteain actualky constraints the y of that point and the rank is deficient. I have to look into the point on ellipse constraint to see why this is happening (today I had very little time available.
My hypothesis is, you are using the parameter t together with xp and yp
t is an arrangement to simplify the equations. It is only used to "decide" from which point on the ellipse the distance is to be calculated, is not a parameter to the matrix above. The equation of the distance from point to ellipse is a funtion of t, t is solved numerically using Newton-Raphson, fixing the point of the ellipse closest to our point (and this works fine, tested). Then the gradient of the distance with respect to all the intervening parameters (not t, t is not a parameter of the equation) are calculated and returned. The solver builds the system matrix based on those gradients as above.
(I have to admit, I never looked into your code. I am not a c++ programmer. The one time, I had a look into the sketcher code, I was totally lost.) Please correct me, if I am wrong.
Your comments are very logical. I think it is not the problem, but I have to find the problem first to know.
The ellipse-simulation-sketch does not have the t-parameter. It based on the pin and string method for an ellipse construction. And it works.
Yep, but I would not know how to get the gradient from there...
By the way: do you know wxMaxima? This program may be helpful in order to generate the needed expressions.
Yes. I use sage, which uses a hundred of opensource packages including maxima... take a look it is totally worth it

https://www.flickr.com/photos/sagescree ... 532209437/
http://www.sagemath.org/links-components.html

Thanks for replying!!
wmayer
Site Admin
Posts: 16665
Joined: Thu Feb 19, 2009 10:32 am

Re: Sketcher: Ellipse support

Postby wmayer » Sat Sep 13, 2014 5:54 pm

A few things I have noticed with the ellipse:
* when trying to drag the ellipse on the perimeter it jumps like crazy
* once you succeed to drag the shape to enlarge the major axis it simply stops after a short distance. No way to enlarge it more.
* trying to enlarge the minor axis doesn't seem to work at all, however shortening is possible
* if you shortened the minor axis and then drag the ellipse (by moving its center) the old shape is restored

* once you set the "angle" constraint to the ellipse you can enlarge as much as you want along the major axis but with a lot of errors in the output window
* enlarging along the minor axis works too but you cannot make it longer than the major axis
=> as said before the ellipse editing shouldn't limit the user. If the user wants to make the minor axis larger then internally the representation of minor and major axis should silently swap
=> having this angle constraint seems unnatural to me and limits the user. IMO this should be solved and handled differently

Limitations:
* in case I want to make a line parallel to the the minor or major axis how would I do that? After trying a couple of minutes I have no way found to do that.
wmayer
Site Admin
Posts: 16665
Joined: Thu Feb 19, 2009 10:32 am

Re: Sketcher: Ellipse support

Postby wmayer » Sat Sep 13, 2014 7:31 pm

I do not blame you. It is too long To me it makes sense, but as a generic functionality that would generate desirable construction lines and points of a geometry constrained to the actual geometry... for those who had technical drawing in high school: Do you know a graphical way of obtaining the focuses of an ellipse, given the ellipse and its center?
For an ellipse with the focus points F1, F2 and an arbitrary point S on the shape:
dist(F1, S) + dist(F2,S) = 2 * a = const (with "a" as the length of the major axis)

So, this means that the distance from the focus to the endpoint on the minor axis is equal to the length of the major axis (because of the symmetry). To construct the focus points make a circle around the endpoint on the minor axis with the radius equal to length of the major axis. The intersection points of this circle and the major axis gives you the focus points.
abdullah
Posts: 3584
Joined: Sun May 04, 2014 3:16 pm

Re: Sketcher: Ellipse support

Postby abdullah » Sun Sep 14, 2014 6:00 am

Welcome back Werner!
A few things I have noticed with the ellipse:
* when trying to drag the ellipse on the perimeter it jumps like crazy
* once you succeed to drag the shape to enlarge the major axis it simply stops after a short distance. No way to enlarge it more.
* trying to enlarge the minor axis doesn't seem to work at all, however shortening is possible
* if you shortened the minor axis and then drag the ellipse (by moving its center) the old shape is restored
The code is very alpha state, meaning that these features are only half implemented, most of them have a placeholder code inhereted from circle modified so that it compiles. The real implementation will come with time, but I have given lower priority to it. First working with whistles, then the whistles...
* once you set the "angle" constraint to the ellipse you can enlarge as much as you want along the major axis but with a lot of errors in the output window
The ellipse angle constraint, being an specific constraint, is pending a major design decision, i.e. how do you continue with the design (specific, general constraints on constructions lines, general constraints in new elements on the ellipse, i.e. Inventor material, quadrants as Normand indicated in this thread). Since then, the angle constraint is there to help me fix the ellipse to see if other constraints (like pointonellipse or line tangent to ellipse converge under angle constraint). More or less the same seems to apply to the separate two major and minor radii specific constraints. This we can discuss separately, not to mix things, in reply to
IMO this should be solved and handled differently
.
enlarging along the minor axis works too but you cannot make it longer than the major axis
=> as said before the ellipse editing shouldn't limit the user. If the user wants to make the minor axis larger then internally the representation of minor and major axis should silently swap
The theory is extremely easy, the practice might not be so easy, as the ellipse might be constrained with other constraints. In the process of "swapping" (call it create a new one or assign the magnitudes new major to major, last major to minor, rotate by 90º), all the other constraints shall be "updated" (the angle with respect to major axis will change by 90 degrees in any way because the major is now at 90 degrees with respect to before). Conflicting constraints might arise, Redundant constraints might arise, the solver will be force to numerically solve all those new constraints at the same time, not one after the other, which might not converge easily or not at all... The user might be left with a situation he did not know how he ran into... Anyway, for me this is an "advanced feature", meaning yes we can look into it, but it is going to be with the lowest prio in the Ellipse support.
Limitations:
* in case I want to make a line parallel to the the minor or major axis how would I do that? After trying a couple of minutes I have no way found to do that.
Reply: currently unimplemented. Constraints... parallel, perpendicular, symmetric, angle, ... are on hold until the design decision is made, as the design decission bears an incredible influence on how the will be implemented. Those were the constraints that exposed (to me) the rigidity of the current approach. By only having selectable a point (center point) and an edge, these requirements could difficultly be met.
IMO this should be solved and handled differently
At the most (Jim and Ulrich ideas) of selecting the edge and supposing this means apply the e.g. parallel constraint to the major axis. But then it would not be possible to have a similar way for the minor axis (though one could do a construction line perpendicular to the major axis through center point, point on ellipse on both extremes and the do a parallel to this line).

Then it comes Normand's clever observation: If necessary do new general constraints but not specific ones, bringing forward the Catia way with quadrants being aligned to lines (so new constraint, "align") and new construction elements for all the forms, also circles ("quadrants").

Then Ulrich, Jim and Normand suggested using construction lines "locked" into the ellipse's major and minor axes, and using those to apply constraints to the ellipse.

As it is not clear to me how Catia exactly works, I ignore which advantages could arise from working with quadrants as opposed to working with two construction lines representing the axes. Aren't those quadrants two construction lines and the alignment constraint a 0 angle or parallel angle constraint??!!

I had myself already thought about construction lines on the axes and the two "construction points" (i.e. points) for the foci. But I did not realise of the huge possibilities of the construction lines on the axes until their user input. Then I thought about making a new toolbar button "expose internal geometry" that would draw all this helper construction lines and relevant points constrained in any Geometry implementing it.

With either of the two (Catia or construction lines), it seems feasible to implement the "angle constraint" as a normal angle (line to line or line to axis). Though it seems to be that Jim's proposal on a "exterior angle" constraint having only a definite positive angle that would avoid flipping should be heard and considered together with it if we decide this way (why? flipping in the sketcher must be removed. It is a great source of "for free" frustration. I would not like to have flipping of the ellipse).

In the context of this construction lines, the Major and Minor radii could also be removed, as the changing of the axes can be effected by changing the length of those construction lines.

From where I stand, the most reasonable approach is the button to make the construction lines and possibly ellipse focuses appear (or maybe, more futuristic, a modal dialog listing the available construction elements and asking the user which ones should be exposed)... phi angle will be constrained as the user wants (horizontal, parallel, angle, vertical,...) by constraining one of the axes. However, I think you are the one that has to make that call. I would not like to code all that and then be told that that was not the right approach. If you need more input from me or other users, just continue the discussion...
For an ellipse with the focus points F1, F2 and an arbitrary point S on the shape:
dist(F1, S) + dist(F2,S) = 2 * a = const (with "a" as the length of the major axis)

So, this means that the distance from the focus to the endpoint on the minor axis is equal to the length of the major axis (because of the symmetry). To construct the focus points make a circle around the endpoint on the minor axis with the radius equal to length of the major axis. The intersection points of this circle and the major axis gives you the focus points.
Yup, I figured it out from Ulrich's "ellipse simulation". Now seems obvious...

You will have a question: If everything up does not have a priority or is on hold, what has it? The answer, tangent and point on object constraints and its convergence.

It is good to have you back...
wmayer
Site Admin
Posts: 16665
Joined: Thu Feb 19, 2009 10:32 am

Re: Sketcher: Ellipse support

Postby wmayer » Sun Sep 14, 2014 12:05 pm

Thanks for the overview of this thread. IMO in order to make a decision it would be good if the users who use a proprietary CAD system can give some more details on how to work with ellipses there.
ulrich1a
Posts: 1958
Joined: Sun Jul 07, 2013 12:08 pm

Re: Sketcher: Ellipse support

Postby ulrich1a » Mon Sep 15, 2014 9:16 pm

As the current approach to define the ellipse equation for the solver seems to have some problems, here is a proposal for a different parameter set to define the ellipse in the solver:

Code: Select all

F1, F2: Focus points
a: major radius
xF1, yF1, xF2, yF2: coordinates of the Focus-points
The major radius and the four coordinates of the focus points do
fully define an ellipse.
These parameters are easier to handle in the solver equations.
Ellipse equation for a point in the ellipse curve: xp, yp

2*a=sqrt((yF2−yp)^2+(xF2−xp)^2)+sqrt((yF1−yp)^2+(xF1−xp)^2)
The equation is the mathematical description of the pins and string method to draw an ellipse. The same as the ellipse-simulation sketch.

The current method to draw an ellipse the first time may stay the same. But then the new parameters needs to be calculated in order to set up the solver equations.

I did attach an zip-archive of an wXMaxima-file with the equation. I tried to make the gradients, but I am not shure, I did it right.

Ulrich
Attachments
Ellipse_Focus_point_equation.wxm.zip
(620 Bytes) Downloaded 20 times
abdullah
Posts: 3584
Joined: Sun May 04, 2014 3:16 pm

Re: Sketcher: Ellipse support

Postby abdullah » Tue Sep 16, 2014 1:47 pm

Hi Ulrich!

I have just downloaded your file.

I see what you have done, but I think it can not be done like this (or I do not see how).

The numerical solver works by minimizing a function, so you need to define a function, that when it gets minimized, a point P is on the Ellipse edge.

I have chosen to use the distance from an Ellipse to a point, but this is just my choice, any other function would do it.

Then you calculate the partials of this function with respect to the 5 parameters of the ellipse and the 2 parameters of point P. This gives you the gradient.
Afterwards, you still have to define a function giving you the error. In my case, the error is the distance itself. So in my case this is very easy.

To give you a hint, in the case of point on line, the distance is used, but defined as the Area of the triangle formed by the two extremes of the line and the point, divided by the length of the line. The partials are calculated, for example, as part(d,P1x)=Part(A,P1x).1/length+A*Part(1/length,P1x)... depending on the variable, some of these terms might be zero (for example the length varies when moving the coordinates of one of the points of the line, but does not when moving the point not belonging to the line, so this latest term is zero).

These days I have not been stopped. I have calculated the partials of the tangent, not using the parametrization of t. This gives me partials that are 600kB of calculation each, which is simply not viable, so I was also looking into going towards a focus based calculation (you can see my-without-t-parameter calculations) here:
http://www.sagenb.org/home/pub/5049

However, in the meantime I have started thinking, that maybe the original implementation (that one where the solver complained about redundant constraints) might not be wrong. I am starting to consider the possibility that the calculation of the rank of the system based on the gradients, might have problems in points where the partials are zero because the function (distance from point to ellipse as a function of the y coordinate of the point) has a minimum at the current position, i.e. the partial is zero because the slope is zero, not because the system does not depend on that parameter...
... this would support why it does work (see Jims example with the line from center to edge+vertical constraint, the line is at the correct position after applying the vertical constraint, however the solver complains about redundant constraint)

... however at this moment these are only conjetures, as I also ask myself why in other cases this does not seem to happen...

I keep looking. Let me know if you manage to get to some point in the meantime...
abdullah
Posts: 3584
Joined: Sun May 04, 2014 3:16 pm

Re: Sketcher: Ellipse support

Postby abdullah » Tue Sep 16, 2014 3:31 pm

Hi Ulrich,

If you want to take a look, this is the (math) work behind the "still-not-properly-working" Point on Ellipse:

http://www.sagenb.org/home/pub/5049

One of the nice things of sage, is that you get the c code for free (as far as you do not use abs function, because is not supported)...

This is the gradient implementation:

Code: Select all

double ConstraintEllipseTangentLine::grad(double *param)
{      
    double deriv=0.;
    if (param == p1x() || param == p1y() ||
        param == p2x() || param == p2y() ||
        param == e1x() || param == e1y() ||
        param == rmaj() || param == rmin() ||
        param == phi()) {
        
        // So:
        // 1. Starting from the parametric equations, solve for t that has the same tangent slope as the line
        //          E'(t)=m => two roots t0, t1
        // 2. Calculate distance from E(t0) and E(t1) to line => select the one (ts) shortest distance to line 
        //    (this is implicit selection by choosing an appropriate initial guess for Newton-Raphson)
        // 3. Calculate the partials of this distance assuming constant t of value ts.
        
        // so first is to select the point (X_0,Y_0) in line to calculate
        double X_1 = *p1x();
        double Y_1 = *p1y();
        double X_2 = *p2x();
        double Y_2 = *p2y();
        double X_c = *e1x();
        double Y_c = *e1y();
        double a = *rmaj();
        double b = *rmin();
        double phip = *phi();
        
        // so first get t for our case (the Eccentric anomally corresponding to the closest point of the ellipse)
        int i=0, maxmitr=100;
        double h, t1, f, df;
        // TODO: ellipse check this error
        double error=1e-20;
        // TODO: ellipse check this for convergence
        
        double X_0 = ( X_1 + X_2 ) / 2;
        double Y_0 = ( Y_1 + Y_2 ) / 2;
        
        double t=atan2(a*(-(X_0-X_c)*sin(phip)+(Y_0-Y_c)*cos(phip)),b*((X_0-X_c)*cos(phip)+(Y_0-Y_c)*sin(phip)));
        
        // Newton-Raphson to get the t0 that corresponds to d.
        for (i=1; i<=maxmitr; i++){
            f=-(-a*sin(phip)*sin(t) +
                b*cos(phip)*cos(t))/(a*sin(t)*cos(phip) + b*sin(phip)*cos(t)) - (Y_1 -
                Y_2)/(X_1 - X_2);
            df=(-a*sin(phip)*sin(t) + b*cos(phip)*cos(t))*(a*cos(phip)*cos(t)
                - b*sin(phip)*sin(t))/pow(a*sin(t)*cos(phip) + b*sin(phip)*cos(t), 2) +
                (a*sin(phip)*cos(t) + b*sin(t)*cos(phip))/(a*sin(t)*cos(phip) +
                b*sin(phip)*cos(t));
            h=f/df;
            t1=t-h;
            if (fabs(h) < error)
                break;
            t=t1;
        }
        
        // t is the solution for the distance at this point
        // TODO: Ellipse implementation: Insert the partials
        
        if (param == p1x()) 
            deriv += ((Y_1 - Y_2)/(X_1 - X_2) + (Y_1 - Y_2)*(-X_1 + X_c +
                a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/pow(X_1 - X_2, 2))*(-Y_1 + Y_c
                + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c +
                a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2))/(sqrt(1 + pow(Y_1
                - Y_2, 2)/pow(X_1 - X_2, 2))*sqrt(pow(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2), 2))) + pow(Y_1 - Y_2, 2)*sqrt(pow(-Y_1
                + Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 +
                X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2), 2))/(pow(1 +
                pow(Y_1 - Y_2, 2)/pow(X_1 - X_2, 2), 3.0L/2.0L)*pow(X_1 - X_2, 3));
        if (param == p1y()) 
            deriv += -(1 + (-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2))*(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2))/(sqrt(1 + pow(Y_1 - Y_2, 2)/pow(X_1 -
                X_2, 2))*sqrt(pow(-Y_1 + Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) -
                (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1
                - X_2), 2))) - (Y_1 - Y_2)*sqrt(pow(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2), 2))/(pow(1 + pow(Y_1 - Y_2, 2)/pow(X_1
                - X_2, 2), 3.0L/2.0L)*pow(X_1 - X_2, 2));
       if (param == p2x()) 
            deriv += -(Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))*(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2))/(sqrt(1 + pow(Y_1 - Y_2, 2)/pow(X_1 -
                X_2, 2))*pow(X_1 - X_2, 2)*sqrt(pow(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2), 2))) - pow(Y_1 - Y_2, 2)*sqrt(pow(-Y_1
                + Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 +
                X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2), 2))/(pow(1 +
                pow(Y_1 - Y_2, 2)/pow(X_1 - X_2, 2), 3.0L/2.0L)*pow(X_1 - X_2, 3));
        if (param == p2y()) 
            deriv += (-X_1 + X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))*(-Y_1 +
                Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c
                + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2))/(sqrt(1 +
                pow(Y_1 - Y_2, 2)/pow(X_1 - X_2, 2))*(X_1 - X_2)*sqrt(pow(-Y_1 + Y_c +
                a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c +
                a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2), 2))) + (Y_1 -
                Y_2)*sqrt(pow(-Y_1 + Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) -
                (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1
                - X_2), 2))/(pow(1 + pow(Y_1 - Y_2, 2)/pow(X_1 - X_2, 2),
                3.0L/2.0L)*pow(X_1 - X_2, 2));
        if (param == e1x()) 
            deriv += -(Y_1 - Y_2)*(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2))/(sqrt(1 + pow(Y_1 - Y_2, 2)/pow(X_1 -
                X_2, 2))*(X_1 - X_2)*sqrt(pow(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2), 2)));
        if (param == e1y()) 
            deriv += (-Y_1 + Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 -
                Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 -
                X_2))/(sqrt(1 + pow(Y_1 - Y_2, 2)/pow(X_1 - X_2, 2))*sqrt(pow(-Y_1 + Y_c
                + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c +
                a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2), 2)));
        if (param == rmaj()) 
            deriv += -(-sin(phip)*cos(t) + (Y_1 - Y_2)*cos(phip)*cos(t)/(X_1 -
                X_2))*(-Y_1 + Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 -
                Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 -
                X_2))/(sqrt(1 + pow(Y_1 - Y_2, 2)/pow(X_1 - X_2, 2))*sqrt(pow(-Y_1 + Y_c
                + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c +
                a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2), 2)));
        if (param == rmin()) 
            deriv += (sin(t)*cos(phip) + (Y_1 - Y_2)*sin(phip)*sin(t)/(X_1 -
                X_2))*(-Y_1 + Y_c + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 -
                Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 -
                X_2))/(sqrt(1 + pow(Y_1 - Y_2, 2)/pow(X_1 - X_2, 2))*sqrt(pow(-Y_1 + Y_c
                + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c +
                a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2), 2)));
        if (param == phi()) 
            deriv += (a*cos(phip)*cos(t) - b*sin(phip)*sin(t) + (Y_1 -
                Y_2)*(a*sin(phip)*cos(t) + b*sin(t)*cos(phip))/(X_1 - X_2))*(-Y_1 + Y_c
                + a*sin(phip)*cos(t) + b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c +
                a*cos(phip)*cos(t) - b*sin(phip)*sin(t))/(X_1 - X_2))/(sqrt(1 + pow(Y_1
                - Y_2, 2)/pow(X_1 - X_2, 2))*sqrt(pow(-Y_1 + Y_c + a*sin(phip)*cos(t) +
                b*sin(t)*cos(phip) - (Y_1 - Y_2)*(-X_1 + X_c + a*cos(phip)*cos(t) -
                b*sin(phip)*sin(t))/(X_1 - X_2), 2)));
    }
    return scale * deriv;
}