Possible "Angle of Incidence" of an RC tail-dragger plane

I am working for an Aero-Design team and we make RC planes that are capable of humanitarian aid. We wished to make a tail dragger and for stability reasons the "Angle of Incidence" (AOI) of the plane has to be between 10 to 20 degrees when it is grounded. Now the AOI depends upon a number of factors such as LHT, LVT, height of landing-gears etc... Thus the number of combination are too large for a human to manually calculate and thus I wrote MATLAB code that is able to consider all the possible combination and permutation and display all the sets of combinations that make the AOI fall within the desired range.

Now the code is extremely long and is filled with too many for loops, this makes slow to process and difficult to interpret. So I wish to know whether or not I could optimize the code, shorten it and make it faster (maybe by pre-defining variables in vectors form)?

Here is the in-depth description of how I developed the logic for this code.

Here is a gist of the logic used:

All the variables that influence the AOI have been marked in the sketch. By using simple trigonometry we can form 3 equations that relate some of these variables with AOI:

Equation 1
$$
\tan(\alpha)= \frac{x}{(V+q)}
$$
Equation 2
$$
\tan(\alpha)= \frac{T}{y}
$$
Equation 3
$$
\tan(\alpha)= \frac{U}{(L+H+q)}
$$

On rearranging the above 3 equations in terms of alpha we get the following equations:

Equation 1 (constraining tail gear):
$$
\sin(\alpha) = \frac{-2V\cdot c + \sqrt{4V^2c^2-4(V^2+x^2)(c^2-x^2)}} {2(V^2+x^2)}
$$

Equation 2 (constraining main landing gear):
$$
\tan(\alpha) = \frac{T}{y}
$$

Equation 3 (constraining the whole system):
$$
\sin(\alpha) = \frac{-(2Hc+2Lc) + \sqrt{(2Hc+2Lc)^2-4(L^2+H^2+2LH+U^2)(C^2- U^2)}}{2(L^2+H^2+2LH+U^2}
$$
On solving these 3 equations simultaneously, we can find the AOI of the plane for given value of the variables.

AoA of a Tail-Dragger

```
clc
clear all
L=100; %Distance between LEw to LEt

In this answer I'll address only one specific improvement you can make: use vectorized operators instead of for loops. This gave me a speedup of about x70. Here is how:

L = 100;         % Distance between LEw to LEt   
H = 20;          % Chord of H-stab  
B = 15;          % Height of the fuselage   
w = -5;          % Height of center of wing integration above the fuselage   
c = 5;           % Clearance of the trailing edge of H-stab from the ground   
k = 0;           % Distance betweem MLG from LEw (+ve if MLG is after LEw)  
x = 7:0.5:10;    % Height of TG  
y = 60:0.5:130;  % Distance btw MLG and TG      
z = 20:0.5:25;   % Height of MLG  
E = 5:0.5:10;    % Height difference of LEw to LEt  

% span each of the vectors along a separate dimension
x_ = x(:);              %  size(x_) = [numel(x),1]
y_ = shiftdim(y(:),-1); %  size(y_) = [1,numel(y)]
z_ = shiftdim(z(:),-2); %  size(z_) = [1,1,numel(z)]
E_ = shiftdim(E(:),-3); %  size(E_) = [1,1,1,numel(E)]

U = bsxfun(@minus,z_,E_)+B+w;  % Height of tail boom above the ground. size(U) = [1,1,numel(z),numel(E)]
T = bsxfun(@minus,U,x_); % Height of TG from the ground.  size(T) = [numel(x),1,numel(z),numel(E)]
V = L+H-y_-k; % Distance btw TG and trailing edge of H-stab.   size(V) = [1,numel(y)]
dis1 = bsxfun(@minus,4*(V*c).^2,4*bsxfun(@times,bsxfun(@plus,V.^2,x_.^2),(c^2-x_.^2))); % size(dis1) = [numel(x),numel(y)]
dis2=(2*H*c+2*L*c)^2-4*(L^2+H^2+2*L*H+U.^2).*(c^2-U.^2);  % size(dis2) = [1,1,numel(z),numel(E)]
inds = bsxfun(@and,dis1>0,dis2>0); % logical indices of where dis1>0 and dis2>0.  size(inds) = [numel(x),numel(y),numel(z),numel(E)]
alpha1 = asin(bsxfun(@plus,-2*V*c,sqrt(dis1))./(2*bsxfun(@plus,V.^2,x_.^2))); % for eqn1.  size(alpha1) = [numel(x),numel(y)]
alpha11 = round(rad2deg(alpha1)); % for eqn1.  size(alpha11) = [numel(x),numel(y)]
alpha2 = atan(bsxfun(@rdivide,T,y_)); % for eqn_2. size(alpha2) = [numel(x),numel(y),numel(z),numel(E)]
alpha22 = round(rad2deg(alpha2)); %for eqn2. size(alpha2) = [numel(x),numel(y),numel(z),numel(E)]  
alpha3 = asin((-(2*H*c+2*L*c)+sqrt(dis2))./(2*(L^2+H^2+2*L*H+U.^2))); %for eqn_3. size(alpha3) = [1,1,numel(z),numel(E)]
alpha3 = round(rad2deg(alpha3)); %for eqn_3. size(alpha3) = [1,1,numel(z),numel(E)]  
alpha = repmat(alpha11,[1,1,numel(z),numel(E)]); % size(alpha) = [numel(x),numel(y),numel(z),numel(E)]
alpha(bsxfun(@ne,alpha11,alpha22) | bsxfun(@ne,alpha22,alpha3) | (inds==0)) = nan;  % this replaces the if statements
inds2 = find(~isnan(alpha) & (alpha>10) & (alpha<20)); % inds2 is a column vector
[inds_x,inds_y,inds_z,inds_E] = ind2sub(size(inds),inds2);

% ther results will be stored in the following matrix, where each row is a result, and the columns are alpha,x,y,z,E
results = [alpha(inds2),x(inds_x)',y(inds_y)',z(inds_z)',E(inds_E)'];

Here is a short explanation: The key idea is to span each of the vectors which you run over in the for loop, along a different direction. We span x along the first dimension, y along the second, z along the third, and E along the fourth dimension. This lets you use the bsxfun function to operate on these vectors giving you all the combinations you need. So for example, in the code above, T contains all combinations of the expression "z+B+w-E-x". Since it depends on the vectors x,z and E, it's dimensions will be [numel(x),1,numel(z),numel(E)] (where numel is the number of elements in the vector). The 1 is there since T does not depend on y.

The if statements in your code are replaced by vectorized conditioning, in the line the starts with alpha(bsxfun(. A drawback of this approach is that I calculated alpha for all combinations of the input, and then I replace with nan the values which do not meet the conditions of your if statements. This means that I waste time on unnecessary calculations. But as I said, I get a big speed improvement on your example so I think we're okay.

Another drawback of the vectorized approach is that it consumes more memory, since all the calculation results are stored at the same time. This may be somewhat mitigated by making sure to clear variables that are not needed any more. But in the case where there is not enough memory to even start the vectorized calculation, a hybrid approach can be taken, where you vectorize (say) x,y and z, but run over E with a for loop.

StackExchange Code Review Q#152221, answer score: 3