Search code examples
gnuplotellipsepolar-coordinates

Gnuplot: Fit an ellipse to a dataset in polar coordinates

For an assignment, I need to plot a inertia ellipsoid of an object. In order to do this, I plot the reciprocal value of the inertia against the angular rotation. After this is done, I need to fit the plot with an ellipse.

But I can't plot my dataset in the deg mode, only in the rad mode and I have no idea why my code doesn't work. Also I need help to plot the ellipse.

Here is my code:

set terminal png
set output 'tisch.png'

set angles radians
set polar
show polar
set parametric
set grid polar
set size square
set trange [0:360]
set rrange [0:0.5]

plot 'A3-daten.txt'

And here is my dataset:

0       0.494012339
30      0.510681467
60      0.461169413
90      0.42190106
120     0.408044505
150     0.442066272
180     0.496961666

Thank you in advance for your help and sorry for my grammatical errors, English is my second/third language and even though I'm able to understand in quite well, I still have sometimes difficulties expressing myself in a understandable way.

Solution

  • In polar mode, gnuplot can fit and plot functions of the form r(t), where t is the angle and r the distance from the origin. This means that we have to express an ellipse using this form.

    Next we consult Wikipedia. Here we have to be careful: The t used by Wikipedia for the parametric representation of an ellipse is not the angle t used by gnuplot, they call it "eccentric anomaly". Using tw instead, the parametric representation is:

    x = a*cos(tw)
y = b*sin(tw)

    We can take these equations and transform them into the polar representation needed for gnuplot (the known transformation from cartesian to polar coordinates):

    r = sqrt( x**2 + y**2 )
  = sqrt( (a*cos(tw))**2 + (b*sin(tw))**2 )

tan(t) = y/x
       = (b*sin(tw)) / (a*cos(tw))

    We need tw, so we solve the second equation:

    sin(tw) / cos(tw) = (a*sin(t)) / (b*cos(t))
          tan(tw) = (a*sin(t)) / (b*cos(t))
               tw = atan2( a*sin(t), b*cos(t))

    The ellipse can be rotated by using t-phi instead of t.

    Now we are done with transforming math equations, and we can start with gnuplot. The script is straight forward:

    datafile = "A3-daten.txt"                                  

set terminal pngcairo                                      
set output "ellipse.png"                                   
set size square                                            

tw(t) = atan2(a*cos(t-phi),b*sin(t-phi))                
r(t) = sqrt( (a*cos(tw(t)))**2 + (b*sin(tw(t)))**2 ) 

set polar
set angles degrees                                         

fit r(t) datafile via a,b,phi                              
plot datafile title datafile ls 7, r(t) title "Fit"

    I arrive at:

    Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 0.408469         +/- 0.001589     (0.3889%)
b               = 0.51369          +/- 0.001782     (0.3469%)
phi             = 21.0673          +/- 0.7251       (3.442%)

    fitted ellipse