Oswald efficiency estimation function

This function lets you use different methods to estimate Oswald efficiency factor, e, for purposes of approximating a simple parabolic drag polar:
C_D = C_D0 + C_L^2/(pi*e*AR), where C_L is lift coefficient.

Based on some or all of the following:
- Wing aspect ratio
- Wing sweep
- Parasite drag coefficient (drag independent of lift)
- Ratio of fuselage diameter (or width) to wing span
- Planform efficiency - usually .98 < u < 1

You can also choose between two common methods as proposed by Raymer ("Aicraft Design: A Conceptual Approach") or Shevell ("Fundamentals of Flight"). OSWALDFACTOR can also calculate using both methods and return the most pessimistic or optimistic result or a weighted average of the two results.

Example: Replicate Shevell's Figure 11.8.
u = .99; df_b = 0.114;

ARx = linspace(0,12,25);
[CD0,AR] = meshgrid([.01 .015 .02 .025],ARx);
e = oswaldfactor(AR,0,'shevell',CD0,df_b,u);
plot(ARx,e); axis([0 12 .6 1]); grid on
xlabel('Aspect Ratio'); ylabel('Efficiency factor, e')
legend('C_{D,0} = 0.01','0.015','0.02','0.025')
text(8.5,.63,'u = 0.99 s = 0.975')

sweepx = linspace(0,40,41);
[AR,sweep] = meshgrid([4;8;12],sweepx*pi/180);
e = oswaldfactor(AR,sweep,'shevell',0,df_b,u)./...
oswaldfactor(AR, 0 ,'shevell',0,df_b,u);
axes('pos',[0.25 0.25 0.3 0.25]);
plot(sweepx,e); axis([0 40 .9 1.02]); grid on
xlabel('Sweep, \Lambda (deg)'); ylabel('e_\Lambda/e_{\Lambda=0}')
legend('AR = 4','8','12','Location','SouthWest')