% Define several symbolic variables. 
% Note that x is r/R 
x=sym('x');  
mu=sym('mu'); 
psi = sym('psi'); 
ut=sym('ut'); 
up=sym('up'); 
beta=sym('beta'); 
beta0=sym('beta0'); 
beta1c=sym('beta1c'); 
beta1s=sym('beta1s'); 
theta=sym('theta'); 
theta0=sym('theta0'); 
thetaw=sym('thetaw'); 
y=sym('y'); 
mu=sym('mu'); 
lambda=sym('lambda'); 
sigma = sym('sigma'); 
a=sym('a'); 
alpha = sym('alpha'); 
cl = sym('cl'); 
 
% Express theta in terms of theta0 and thetaw 
 
theta=theta0+thetaw*x; 
 
 
% Define UT already non-dimensionalized by OMEGA * R 
ut=x+mu*sin(psi); 
 
% Define UP , non-dimensionalized by OEMGA * R 
% Recall UP = lambda * (OMEGA*R)+ r * beta_dot + Vinf * beta * cos(psi) 
 
up=lambda+x*(beta1s*cos(psi)-beta1c*sin(psi))... 
   +mu*(beta0+beta1c*cos(psi)+beta1s*sin(psi))*cos(psi); 
 
% set up the integrand for Qi. Multiply by sigma = bc/pi*R already 
% Rho will cancel out because it appears in the integrand, and in the denominator 
% when we non-dimensionalize Qi. 
 
theta=theta0+thetaw*x; 
alpha = theta - up/ut; 
cl = a * alpha; 
y=0.5*sigma* cl* ut * up * x; 
 
% Intgerate from x=0 to 1. Store answer in a temporary variable called part1 
 
part1=sym('part1'); 
part1=int(y,x,0,1); 
 
% Take azimuthal average by integrating with respect to psi, divide by 2*pi 
% recall we have already picked up the pi in defining sigma above. 
 
cqi=sym('cqi'); 
 
cqi= int(part1,psi,0,2*pi)/(2*pi); 
simplify(cqi)