NAL TR-797 cf. http://airex.tksc.jaxa.jp/pl/dr/NALTR0797000

TR797.m

% -----------
% Optimum Design of Nonplanar wings-Minimum Induced Drag 
% for A Given Lift and Wing Root Bending Moment (NAL TR-797)
% 
% Created by Takahiro Inagawa on 2013-2-28.
% Copyright (c) 2013 Takahiro Inagawa. All rights reserved.
% -----------
clear all
tic;  %Computation time measurement

% ---- Constant ----
Lift = 1000;	% Lift[N]
Uinf = 7.2;	    % Aircraft speed[m/s]
rho = 1.2;	    % Air density[kg/m3]
span = 30;	    % Span width[m]
le = span / 2;	% length from root to tip
N = 100;	    % Partition number
beta = 0.85;	% Coefficient related to the bendin moment
delta_S = linspace(le / N / 2, le / N / 2, N);	% half of partition, Equal distance


% ---- Geometric Condition ----
% yz plane
% y axis is vertical, z axis is horizontal.
% fai is the local dihedral angel.
% y is the center of the partition.
% -----------------------------
y = linspace(delta_S(1), le - delta_S(N), N);
z = linspace(0, 0, N);
fai = linspace(0, 0, N);

% ---- Geometric variable number ----
% Variables required to seek the Q.
for i = 1:N
	for j = 1:N
		ydash(i,j)  =   (y(i) - y(j)) .* cos(fai(j)) + (z(i) - z(j)) .* sin(fai(j));
		zdash(i,j)  = - (y(i) - y(j)) .* sin(fai(j)) + (z(i) - z(j)) .* cos(fai(j));
		y2dash(i,j) =   (y(i) + y(j)) .* cos(fai(j)) - (z(i) - z(j)) .* sin(fai(j));
		z2dash(i,j) =   (y(i) + y(j)) .* sin(fai(j)) + (z(i) - z(j)) .* cos(fai(j));
		
		R2puls(i,j)  = (ydash(i,j) - delta_S(j)).^2 + zdash(i,j).^2;
		R2minus(i,j) = (ydash(i,j) + delta_S(j)).^2 + zdash(i,j).^2;
		Rdash2puls(i,j)  = (y2dash(i,j) + delta_S(j)).^2 + z2dash(i,j).^2;
		Rdash2minus(i,j) = (y2dash(i,j) - delta_S(j)).^2 + z2dash(i,j).^2;
	end
end

for i = 1:N
	for j = 1:N
		Q(i,j) = - 1 / (2 * pi) *(((ydash(i,j) - delta_S(j)) / R2puls(i,j) 
		- (ydash(i,j) + delta_S(j)) / R2minus(i,j)) * cos(fai(i) - fai(j))
		+ (zdash(i,j) / R2puls(i,j) - zdash(i,j) / R2minus(i,j)) * sin(fai(i) - fai(j))
		+ ((y2dash(i,j) - delta_S(j)) / Rdash2minus(i,j) 
		- (y2dash(i,j) + delta_S(j)) / Rdash2puls(i,j)) * cos(fai(i) + fai(j))
		+ (z2dash(i,j) / Rdash2minus(i,j) - z2dash(i,j) / Rdash2puls(i,j))
		* sin(fai(i) + fai(j)));
	end
end

% ---- Normalization ----
% Variables required to seek q.
delta_sigma = delta_S ./ le;
eta = y ./ le;
etadash = ydash ./ le;
eta2dash = y2dash ./ le;
zeta = z ./ le;
zetadash = zdash ./ le;
zeta2dash = z2dash ./ le;
gamma2puls = R2puls / (le^2);
gamma2minus = R2minus / (le^2);
gammadash2puls = Rdash2puls / (le^2);
gammadash2minus = Rdash2minus / (le^2);
for i = 1:N
	for j = 1:N
		q(i,j) = - 1 / (2 * pi) *(((etadash(i,j) - delta_sigma(j)) / gamma2puls(i,j) 
		- (etadash(i,j) + delta_sigma(j)) / gamma2minus(i,j)) * cos(fai(i) - fai(j))
		+ (zetadash(i,j) / gamma2puls(i,j) - zetadash(i,j) / gamma2minus(i,j))
		* sin(fai(i) - fai(j))
		+ ((eta2dash(i,j) - delta_sigma(j)) / gammadash2minus(i,j) 
		- (eta2dash(i,j) + delta_sigma(j)) / gammadash2puls(i,j)) * cos(fai(i) + fai(j))
		+ (zeta2dash(i,j) / gammadash2minus(i,j) - zeta2dash(i,j) / gammadash2puls(i,j))
		* sin(fai(i) + fai(j)));
	end
end

% ---- elliptic loading aerodynamic force ----
% Vn is Induced vertical velocisy.
% Vn is constant when elliptical circulation distribution.
BendingMoment_elpl = 2 / 3 / pi * le * Lift;
InducedDrag_elpl = Lift^2 / (2 * pi * rho * Uinf^2 * le^2);
Vn_elpl = Lift / (2 * pi * rho * Uinf * le^2);

% ---- Creating the optimization equation ----
c = 2 * cos(fai) .* delta_sigma;
b = 3 * pi / 2 *(eta .* cos(fai) + zeta .* sin(fai)) .* delta_sigma;
for i = 1:N
	for j = 1:N
		A(i,j) = pi * q(i,j) .* delta_sigma(i);
	end
end

% ---- solve optimization problem ----
AAA = A + A';
ccc = -c;
bbb = -b;
LeftMat = vertcat([AAA ccc' bbb'], [ccc 0 0], [bbb 0 0]);
RightMat = vertcat(linspace(0,0,N)', -1, -beta);
SolveMat = LeftMat \ RightMat;
g = SolveMat(1:N);
mu = SolveMat(N+1:N+2);	% Lagrange multiplier


% ---- After Solve ----
% efficient : span efficiency
% Gamma : Local circulation
% InducedDrag : Sum of induced drag
% Vn : Local induced vertical velocisy
% Lift0_elpl : Lift at root culcurated by area of the ellipse circulation distribution
% Gamma0_elpl : Circulation at root
% Gamma_elpl : Analytical ellipse circulation distribution @ beta = 1
% ---------------------
efficiency_inverse = g' * A * g;
efficiency = 1 / efficiency_inverse;
Gamma = g * Lift / (2 * le * rho * Uinf);
InducedDrag = efficiency_inverse * InducedDrag_elpl;

Local_Lift = 4 * rho * Uinf * Gamma' .* cos(fai);
Vn = zeros(1,N);
for i = 1:N
	for j = 1:N
		Vn(i) += Q(i,j) .* Gamma(j);
	end
end
Local_InducedDrag = rho * Gamma' .* Vn;

% ---- Aerodynamic force when Elliptical cierculation distribution----
Lift0_elpl = 2 * Lift / pi / le;
Lift_elpl = 4 * Lift0_elpl * sqrt(1 - (y / le).^2);
Gamma0_elpl = Lift0_elpl / (rho * Uinf * cos(fai(1)));
Gamma_elpl = Gamma0_elpl * sqrt(1 - (y / le).^2);
Local_InducedDrag_elpl = 2 * rho * Gamma_elpl * Vn_elpl;

% ---- Bending Moment ----
Local_BendingMoment = zeros(1,N);
Local_BendingMoment_elpl = zeros(1,N);
for i = 1:N
	tmp1 = 0;
	tmp2 = 0;
	for j = i:N
		tmp1 += Local_Lift(j) * (y(j) - y(i));
		tmp2 += Lift_elpl(j) * (y(j) - y(i));
	end
	Local_BendingMoment(i) = tmp1;
	Local_BendingMoment_elpl(i) = tmp2;
end


% ---- Plot ----
% figure 1 : Circulation
%        2 : Lift
%        3 : Benging Moment
%        4 : Induced vertical velocity
%        5 : Induced Drag
figure(1)
plot(y, Gamma, y, Gamma_elpl);
xlabel('position [m]');ylabel('Gamma');
str = sprintf('\beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Circulation');
xlim([0 ,max(y) * 1.05]);
grid on;
print("Circulation.jpg","-djpg","-S800,500","-r300");

figure(2)
plot(y, Local_Lift, y, Lift_elpl);
xlabel('position [m]'); ylabel('Lift[N]');
str = sprintf('\beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Lift');
xlim([0 ,max(y) * 1.05]);
grid on;
print("Lift.jpg","-djpg","-S800,500","-r300");

figure(3)
plot(y, Local_BendingMoment, y, Local_BendingMoment_elpl);
xlabel('position [m]'); ylabel('Bending Moment [Nm]');
str = sprintf('beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Bending Moment');
xlim([0, max(y) * 1.05]);
grid on;
print("Bending Moment.jpg","-djpg","-S800,500","-r300");

figure(4)
plot(y, linspace(Vn_elpl, Vn_elpl, N), y, Vn);
xlabel('position [m]'); ylabel('Induced vertical velocity [m/s]');
str = sprintf('beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Induced vertical velocity');
xlim([0, max(y) * 1.05]);
grid on;
print("Induced vertical velocity.jpg","-djpg","-S800,500","-r300");

figure(5)
plot(y, Local_InducedDrag, y, Local_InducedDrag_elpl);
xlabel('position [m]'); ylabel('Induced Drag [N]');
str = sprintf('beta = %.2f', beta);
legend(str, 'ellipse circulation distribution');
title('Induced Drag');
xlim([0, max(y) * 1.05]);
grid on;
print("Induced Drag.jpg","-djpg","-S800,500","-r300");

% ---- Display Input and Output ----
disp('---- Input ----')
str = sprintf('Lift             : %.1f', Lift);
disp(str)
str = sprintf('Aircraft speed   : %.2f', Uinf);
disp(str)
str = sprintf('Wing span width  : %.2f', span);
disp(str)
str = sprintf('Partition number : %d', N);
disp(str)
str = sprintf('beta             : %.3f', beta);
disp(str)
disp('---- Output ----')
str = sprintf('Induced Drag     : %.4f', InducedDrag);
disp(str)
str = sprintf('efficiency       : %.4f', efficiency);
disp(str)
str = sprintf('');
disp(str)
disp('---- End ----')
toc

https://twitter.com/kikatyan/status/314421508313317379

・改行しても式が終わらないところは ... をつける．
・ダブルクォーテション　"　　->　シングルクオーテション '
・matlabで+=演算子はつかえない？ので a = a+xとなおす．