curve fitting and graphing.

Here’s how i do the curve fits and resulting graphing.

I use the free program ‘gnuplot’ which is pretty powerful, can do curve-fits, and graph everything.

Steps:

1: measure your curve and get some points. I created a file called ‘measured.dat’:

0  .0044
.1 .0123
.2 .0218
.3 .0382
.4 .0636
.5 .098
.6 .157
.7 .221
.8 .300
.9 .404
1 .495

2)  fit the curve. I use a cubic parabola: Ax3+Bx2+Cx+D

shell$ gnuplot
...
gnuplot> f(x)=a*x*x*x+b*x*x+c*x+d
gnuplot> fit f(x) './measured.dat' via a,b,c,d
...
Final set of parameters            Asymptotic Standard Error
=======================            ==========================a               = 0.142852         +/- 0.07749      (54.25%)
b               = 0.385536         +/- 0.1181       (30.63%)
c               = -0.035066        +/- 0.0493       (140.6%)
d               = 0.00833846       +/- 0.005414     (64.93%)
So now you have your coefficients: a,b,c,d

3) plug them into a graphing file. here’s mine:

# let's call this file "fit.gnuplot"
# you can keep the body of the file the same, and just
# modify a,b,c,d each time you run a new curve
a=0.142852
b=0.385536
c=-0.035066
d=0.00833846

#Cubic Parabola
#C=ax3 + bx2 + cx + d
B(x)=a*x**3+b*x**2+c*x+d
# this trick limits the curve to x=0 through x=1.05
C(x)= (x>0 && x<1.05) ? B(x) : 1/0

# tangent slope, derivative of cantilever
# D = B'=3ax2 + 2bx + c
D(x)=3*a*x**2+2*b*x+c

#focus (reflected) line slope (ts = tangent slope)
# y = 1/(TAN(2*(ATAN(1/ts))))
# incidence angle (0 = directly above,
#                  + means low in the sky
#                  - means high in the sky

# radians = degrees x PI / 180
# precalculate one, just showing the formula, not used
ideg= 0;
irad = ideg*pi/180

## slope of reflected ray
F(ts)=1/tan(2*atan(1/ts)-irad)
# focus y-int
G(tp)=C(tp)-(tp*F(D(tp)))

## where will the graph be drawn?
set xrange [-.5:1.2]
set yrange [-.1:1.2]

## setup graphing
set multiplot
# offset (first ray will be ioff+istep)
ioff=-3
istep=1.5
# each ray will offset on x axis by this much
xstep=.1
# arrow number
ar=0
# tangent point
tp=0
ideg=ioff
irad = ideg*pi/180

## how far to the left of the x-axis to extend arrows?
xp = -.35

## increment various counters
ar=ar+1; tp=tp+xstep

## downward ray
set arrow ar from tp,5 to tp,C(tp) lt rgb 'black'

## normal reflected ray
ideg=ideg+istep; irad = ideg*pi/180; ar=ar+1;
set arrow ar from tp,C(tp) to xp,F(D(tp))*xp+G(tp) lt rgb 'red'

## plus 1 step
ideg=ideg+istep; irad = ideg*pi/180; ar=ar+1;
set arrow ar from tp,C(tp) to xp,F(D(tp))*xp+G(tp) lt rgb 'black'

## plus another step
ideg=ideg+istep; irad = ideg*pi/180; ar=ar+1;
set arrow ar from tp,C(tp) to xp,F(D(tp))*xp+G(tp) lt rgb 'green'

## repeat this stanza as many times as you need
## next set
ideg=ioff; irad = ideg*pi/180
ar=ar+1; tp=tp+xstep
set arrow ar from tp,5 to tp,C(tp)  lt rgb 'black'   
ideg=ideg+istep; irad = ideg*pi/180; ar=ar+1;
set arrow ar from tp,C(tp) to xp,F(D(tp))*xp+G(tp) lt rgb 'red'
ideg=ideg+istep; irad = ideg*pi/180; ar=ar+1;
set arrow ar from tp,C(tp) to xp,F(D(tp))*xp+G(tp) lt rgb 'black'
ideg=ideg+istep; irad = ideg*pi/180; ar=ar+1;
set arrow ar from tp,C(tp) to xp,F(D(tp))*xp+G(tp) lt rgb 'green'

set label "  red: 1.5 degrees high\nblack: perfectly aimed\ngreen: 1.5 degrees low" at -.4,.2

set key at 1,.05
plot C(x) lt 2 title 'best fit line for measured sheet'
set key at 1,0
plot 'measured.dat' with points title 'measured points'

unset multiplot
pause -1

Now you have a graph definition that will do all the fancy arrows and such.

4) Show the graph:

shell$ gnuplot fit.gnuplot

Here’s what it makes:

You’ll need to duplicate the ‘repeat this stanza’ section once more for each point without a reflection on it.

About these ads