how to calculate for a range of values for different initial conditions

function [theta]=theta_beta_M(beta,M,gamma)
% return theta for beta-theta-M relationship for oblique shock
beta=linspace(0,(pi/2),90)
M=[1.25, 2, 6, 10)
gamma= 1.4

%cut off at Mach wave angle
if (beta<=asin(1./M)) theta=0; return; end

theta=atan(2.*cot(beta).*((M.*sin(beta)).^2-1)./(M.^2.*(gamma+cos(2.*beta))+2));

i want to be able to calculate theta values for each value of M for the range of beta.

Answer 1

[Walter Roberson]

function [theta]=theta_beta_M(beta,M,gamma)
% return theta for beta-theta-M relationship for oblique shock
if nargin < 1; beta=linspace(0,(pi/2),90); end
if nargin < 2; M=[1.25, 2, 6, 10); end
if nargin < 3; gamma= 1.4; end

[beta,M] = ndgrid(beta, M);

theta = zeros(size(beta));

%cut off at Mach wave angle
mask = beta > asin(1./M);

theta(mask) = atan(2.*cot(beta(mask)).*((M(mask).*sin(beta(mask))).^2-1)./(M(mask).^2.*(gamma+cos(2.*beta(mask)))+2));

The result will be numel(beta) by numel(M) -- so one column for each different M value.

Answer 2

[David Hill]

beta=linspace(0,(pi/2),90);
m=[1.25, 2, 6, 10];
[B,M]=meshgrid(beta,m);
gamma= 1.4;
theta=atan(2.*cot(B).*((M.*sin(B)).^2-1)./(M.^2.*(gamma+cos(2.*B))+2));
theta(B<=asin(1./M))=0

theta = 4×90

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0143 0.0403 0.0644 0.0870 0.1083 0.1287 0.1483 0.1673 0.1856 0.2035 0.2210 0.2381 0.2549 0.2713 0.2875 0.3035 0.3192 0.3347 0.3500 0.3651 0 0 0 0 0 0 0.0092 0.0349 0.0579 0.0790 0.0987 0.1175 0.1356 0.1530 0.1701 0.1867 0.2031 0.2191 0.2350 0.2506 0.2661 0.2814 0.2966 0.3116 0.3265 0.3413 0.3559 0.3705 0.3849 0.3992