#! /usr/bin/octave -qf
# Octave script for bem calculation
# Grant Ingram (c) 2007

# this does just one radial element - mainly for demonstration purposes...

# blade details
r = 1; gamma = 66; c = 0.17; r_tip = 4;
printf( "Element: r:%5.2f [m] gamma:%5.2f [deg] c:%5.2f [m] \n",r,gamma,c );

# wind details
V = 10; B = 3; tsr = 4; 
omega = tsr * V / r_tip; # note blade speed is derived
sigma = (B * c) / (2 * pi * r);
tsr_r = omega * r / V;
printf("V: %5.2f [m/s] Omega: %5.2f [rad/s] tsr_r: %5.2f [-]\n",V,omega,tsr_r);

# guess initial value
beta = pi/2 - (2.0/3) *atan(1.0/tsr_r);
axial = 1 + ((4 * cos(beta) *cos(beta))/(sigma * 1* sin(beta) ));
axial = 1.0 / axial;
angular = (1 - 3.0*axial)/(4.0*axial -1);
printf("Initial values: axial: %5.2f [-] angular: %5.2f [-]\n",axial,angular);

# iterate
error_axial = 10; error_angular = 10;
while (abs(error_axial) > 0.000001 || abs(error_angular) > 0.000001)
    old_axial = axial; old_angular = angular;
    beta = atan ( (tsr_r * ( 1.0 + angular)) / ( 1.0 - axial));
    i = gamma - (beta * 180.0 / pi); # convert to degrees here
    CL = 0.084*i; # find new values of CL,CD from the aerofoil curve
    axial = 1.0 / ((1 + (4 * cos(beta)*cos(beta))/(sigma * CL * sin(beta))));
    angular = (( sigma * CL ) / (4* tsr_r * cos(beta)))* (1-axial);
    error_axial = axial - old_axial;
    error_angular = angular - old_angular;            
endwhile

printf("Final values: axial: %5.2f [-] angular: %5.2f [-]\n",axial,angular);
printf("Flow params: i: %5.2f [deg] beta: %5.2f [deg] CL %5.2f [-]\n",i,beta*180/3.14,CL);
