To all those interested in sink rate polars, follows is a C++ program I have written to fit two polar equations through experimetal polar data. Curve fitting is useful in severals ways; it provides a means of smoothing the data, provides the the polar equations which can be used for speed to fly theory, and finnally is a useful technique for quantifying the drag produced by gliders. I will eventually convert the following C++ program into Java and put it on a web site for all to play with...
Enjoy!
--------------------start of C++ file ---------------------------------------
/********************************************************/
/* */
/* PolarFit.C */
/* */
/* A C++ program for calculation sink rate */
/* polar coefficients from an ASCII file of */
/* a speed (mph) against sink rate (fpm). */
/* */
/* The two polar equations supported are :- */
/* */
/* sink rate = a*V^2 + b*V + c */
/* and */
/* sink rate = i/V + p*V^3 */
/* */
/* The coefficients a,b,c & i,p are derived */
/* by least squared algorithms. */
/* */
/* All code and algorithms written by */
/* */
/* Robert Osborn, */
/* June 1996. */
/* */
/* ---------------------------------------------------- */
/* */
/* To compile on an SGI/Irix5.3 I used :- */
/* */
/* DCC PolarFit.C -o polarfit -lm */
/* */
/* To compile on an Ultra Sparc/Solaris2.5.1 :- */
/* */
/* CC PolarFit.C -o polarfit -lm */
/* */
/* ---------------------------------------------------- */
/* */
/* To run the program I used :- */
/* */
/* polarfit < PolarData.AT */
/* */
/* Where PolarData.AT is an ASCII file in form :- */
/* */
/* 21 -184 */
/* 23 -183 */
/* 25 -188 */
/* 27 -199 */
/* */
/* speed sink rate, fpm, */
/* mph +ve upwards */
/* */
/* */
/********************************************************/
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>
const int MAXIUMUM_NO_SAMPLES = 1000;
const double AIRSPEED_UNITS_RATIO = 1.0/2.25;
const double SINKRATE_UNITS_RATIO = 1.0/197;
// parabolic polar equation coefficients.
double a_co;
double b_co;
double c_co;
// parabolic sink rate eqaution.
double sink_rate_abc(double v)
{
return a_co*v*v + b_co*v + c_co;
}
// function for calculating best glide speed in headwind and sink.
double best_glide_abc(double wind,double sink)
{
return -wind + sqrt(wind*wind-(b_co*wind-c_co-sink)/a_co);
}
// induced drag & profile drag polar equation coefficients.
double i_co;
double p_co;
// induced drag & profile drag polar sink rate equation.
double sink_rate_ip(double v)
{
return i_co/v + p_co*v*v*v;
}
// crude iterative function for calculating best glide speed.
double best_glide_ip(double wind, double sink)
{
double v_lower = 1;
double v_upper = 30;
int interations = 0;
while ((v_upper-v_lower)>0.01 && interations<100000) {
double v_mid=(v_lower+v_upper)/2.0;
double ds_dv_lower = (v_lower+wind)
*(-i_co/(v_lower*v_lower)+3*p_co*v_lower*v_lower)
-i_co/v_lower-p_co*v_lower*v_lower*v_lower-sink;
double ds_dv_mid = (v_mid+wind)
*(-i_co/(v_mid*v_mid)+3*p_co*v_mid*v_mid)
-i_co/v_mid-p_co*v_mid*v_mid*v_mid-sink;
double ds_dv_upper = (v_upper+wind)
*(-i_co/(v_upper*v_upper)+3*p_co*v_upper*v_upper)
-i_co/v_upper-p_co*v_upper*v_upper*v_upper-sink;
if (ds_dv_lower*ds_dv_mid<0) v_upper = v_mid;
if (ds_dv_upper*ds_dv_mid<0) v_lower = v_mid;
interations++;
}
return (v_lower+v_upper)/2.0;
}
int main(int noParams, int* paramList)
{
double vList[MAXIUMUM_NO_SAMPLES], sList[MAXIUMUM_NO_SAMPLES];
int noValues = 0;
// read in data samples
while ((!cin.eof()) && (noValues<MAXIUMUM_NO_SAMPLES)) {
float v,s;
cin >> v >> s ;
vList[noValues] = v*AIRSPEED_UNITS_RATIO;
sList[noValues] = fabs(s)*SINKRATE_UNITS_RATIO;
noValues++;
}
noValues--;
// calculate a,b,c.
double M=0;
double N=0;
double O=0;
double P=0;
double Q=0;
double R=0;
double S=0;
double T=0;
for(int i=0;i<noValues;i++) {
double v = vList[i];
double s = sList[i];
M += v*v*v*v;
N += v*v*v;
O += v*v;
P += v;
Q += -s*v*v;
R += -s*v;
S += -s;
T += 1.0;
}
a_co = ((Q*T-S*O)*(N*P-O*O)-(Q*P-R*O)*(N*T-P*O))/
((M*P-N*O)*(N*T-P*O)-(M*T-O*O)*(N*P-O*O));
b_co = ((S*O-Q*T)+a_co*(O*O-M*T))/(N*T-P*O);
c_co = (-S - a_co*O - b_co*P)/T ;
// calculate min sink speed and min sink.
double Vmin_abc = -b_co/(2*a_co);
double Smin_abc = a_co*Vmin_abc*Vmin_abc + b_co*Vmin_abc + c_co;
// calculate max glide speed and max LD.
double Vmax_abc = sqrt(c_co/a_co);
double LDmax_abc = Vmax_abc /
(a_co*Vmax_abc*Vmax_abc + b_co*Vmax_abc + c_co);
// calculate i,p.
double A=0;
double B=0;
double C=0;
double D=0;
double E=0;
double F=0;
for(i=0;i<noValues;i++) {
double v = vList[i];
double s = sList[i];
A += v*v*v*v*v*v;
B += v*v;
C += -s*v*v*v;
D += v*v;
E += 1/(v*v);
F += -s/v;
}
i_co = (F*A-C*D)/(B*D-E*A);
p_co = (F*B-C*E)/(A*E-D*B);
// calculate min sink speed and min sink.
double Vmin_ip = pow(i_co/(3.0*p_co),0.25);
double Smin_ip = fabs(i_co/Vmin_ip + p_co*Vmin_ip*Vmin_ip*Vmin_ip);
// calculate max glide speed and max LD.
double Vmax_ip = pow(i_co/p_co,0.25);
double LDmax_ip = Vmax_ip/(i_co/Vmax_ip + p_co*Vmax_ip*Vmax_ip*Vmax_ip);
// calculate errors;
double abc_error = 0;
double ip_error = 0;
for(i=0;i<noValues;i++) {
double v = vList[i];
double s = sList[i];
double s_abc = a_co*v*v + b_co*v + c_co;
abc_error += fabs((s-s_abc)/s_abc);
double s_ip = p_co*v*v*v + i_co/v;
ip_error += fabs((s-s_ip)/s_ip);
}
abc_error /= (double) noValues;
ip_error /= (double) noValues;
// output results.
cout << "A = " << a_co << " B = " << b_co << " C = " << c_co << endl;
cout << "Average % abc_error = " << abc_error << endl;
cout << "Min sink speed " << Vmin_abc/AIRSPEED_UNITS_RATIO
<< " Min sink rate " << Smin_abc/SINKRATE_UNITS_RATIO << endl;
cout << "Max glide speed " << Vmax_abc/AIRSPEED_UNITS_RATIO
<< " Max Glide " << LDmax_abc << endl;
double v_abc_s400 = best_glide_abc(0,400*SINKRATE_UNITS_RATIO);
cout << "Best glide speed in 400fpm down "
<< v_abc_s400/AIRSPEED_UNITS_RATIO;
cout << " LD_abc "
<< v_abc_s400/(sink_rate_abc(v_abc_s400)+400*SINKRATE_UNITS_RATIO);
cout << " LD_ip "
<< v_abc_s400/(sink_rate_ip(v_abc_s400)+400*SINKRATE_UNITS_RATIO)
<< endl;
double v_abc_l200 = best_glide_abc(0,-200*SINKRATE_UNITS_RATIO);
cout << "Best glide speed in 200fpm lift "
<< v_abc_l200/AIRSPEED_UNITS_RATIO;
cout << " LD_abc "
<< v_abc_l200/(sink_rate_abc(v_abc_l200)-200*SINKRATE_UNITS_RATIO);
cout << " LD_ip "
<< v_abc_l200/(sink_rate_ip(v_abc_l200)-200*SINKRATE_UNITS_RATIO)
<< endl;
cout << endl;
cout << "I = " << i_co << " P = " << p_co << endl;
cout << "Average % ip_error = " << ip_error << endl;
cout << "Min sink speed " << Vmin_ip/AIRSPEED_UNITS_RATIO
<< " Min sink rate " << Smin_ip/SINKRATE_UNITS_RATIO << endl;
cout << "Max glide speed " << Vmax_ip/AIRSPEED_UNITS_RATIO
<< " Max Glide " << LDmax_ip << endl;
double v_ip_s400 = best_glide_ip(0,400*SINKRATE_UNITS_RATIO);
cout << "Best glide speed in 400fpm down "
<< v_ip_s400/AIRSPEED_UNITS_RATIO;
cout << " LD_ip "
<< v_ip_s400/(sink_rate_ip(v_ip_s400)+400*SINKRATE_UNITS_RATIO);
cout << " LD_abc "
<< v_ip_s400/(sink_rate_abc(v_ip_s400)+400*SINKRATE_UNITS_RATIO)
<< endl;
double v_ip_l200 = best_glide_ip(0,-200*SINKRATE_UNITS_RATIO);
cout << "Best glide speed in 200fpm lift "
<< v_ip_l200/AIRSPEED_UNITS_RATIO;
cout << " LD_ip "
<< v_ip_l200/(sink_rate_ip(v_ip_l200)-200*SINKRATE_UNITS_RATIO);
cout << " LD_abc "
<< v_ip_l200/(sink_rate_abc(v_ip_l200)-200*SINKRATE_UNITS_RATIO)
<< endl;
cout << endl;
cout << setw(12) << "Airspeed"
<< setw(12) << "Sink Rate"
<< setw(12) << "Sink_abc"
<< setw(12) << "Sink_ip"
<< endl;
cout << setw(12) << "^^^^^^^^"
<< setw(12) << "^^^^^^^^^"
<< setw(12) << "^^^^^^^^"
<< setw(12) << "^^^^^^^"
<< endl;
for(i=0;i<noValues;i++) {
double v = vList[i];
double s = sList[i];
double s_abc = a_co*v*v + b_co*v + c_co;
double s_ip = p_co*v*v*v + i_co/v;
cout << setw(12) << setprecision(4) << v/AIRSPEED_UNITS_RATIO
<< setw(12) << s/SINKRATE_UNITS_RATIO;
cout << setw(12) << setprecision(4) << s_abc/SINKRATE_UNITS_RATIO
<< setw(12) << s_ip/SINKRATE_UNITS_RATIO << endl;
}
return 0;
}
----------- end of C++ file -------------------------------------------------
--
Robert Osborn, Glasgow, Scotland. robert@mvel.demon.co.uk
Race 2+, Davron, Airwave Klassic 12m, Winglets and cool pair of shades ;)