# I want to obtain value of "W" for different combination of values of (psi,r) example. "W" for (psi=0.0001,r = 2.1891) ; (psi=0.001,r = 2.1940); (psi=0.01,r = 2.2373) and so on

2 views (last 30 days)
AVINASH SAHU on 22 Jun 2022
Commented: AVINASH SAHU on 22 Jun 2022
psi = 0.001;
alpha = (12*psi)^(1/3);
r = 2.1940;
l1 = (r + alpha)*(r + alpha);
l1a = r*r - r*alpha + alpha*alpha;
L1 = log(l1/l1a);
t1 = (2*r - alpha) / (alpha*sqrt(3));
T1 = atan(t1);
l2 = (1 + alpha)*(1 + alpha);
l2a = (1*1) - (1*alpha) + (alpha*alpha);
L2 = log(l2/l2a);
t2 = (2*1 - alpha) / (alpha*sqrt(3));
T2 = atan(t2);
W1 = 4 * sqrt(3) * (T1 - T2) * log((r+alpha)/(1+alpha));
W2 = (L1 - L2)*log((r*r - r*alpha + alpha*alpha)/(1 - alpha + alpha*alpha));
W3 = 4*(T1 -T2)*(T1 -T2);
W4 = (L1 - L2) + 2 * sqrt(3) * (T1 - T2);
W5 = 3 ./ ((r - 1) * (r - 1));
W = ((W1 + W2 -W3) * W5)/W4
Dyuman Joshi on 22 Jun 2022
Edited: Dyuman Joshi on 22 Jun 2022
Make it a function and call from workspace.
%calling individually
W1=calculateW(0.0001,2.1891)
W1 = 0.1602
W2=calculateW(0.001,2.194)
W2 = 0.1595
W3=calculateW(0.01,2.2373)
W3 = 0.1534
function W = calculateW(psi,r)
alpha = (12*psi)^(1/3);
l1 = (r + alpha)*(r + alpha);
l1a = r*r - r*alpha + alpha*alpha;
L1 = log(l1/l1a);
t1 = (2*r - alpha) / (alpha*sqrt(3));
T1 = atan(t1);
l2 = (1 + alpha)*(1 + alpha);
l2a = (1*1) - (1*alpha) + (alpha*alpha);
L2 = log(l2/l2a);
t2 = (2*1 - alpha) / (alpha*sqrt(3));
T2 = atan(t2);
W1 = 4 * sqrt(3) * (T1 - T2) * log((r+alpha)/(1+alpha));
W2 = (L1 - L2)*log((r*r - r*alpha + alpha*alpha)/(1 - alpha + alpha*alpha));
W3 = 4*(T1 -T2)*(T1 -T2);
W4 = (L1 - L2) + 2 * sqrt(3) * (T1 - T2);
W5 = 3 ./ ((r - 1) * (r - 1));
W = ((W1 + W2 -W3) * W5)/W4;
end
AVINASH SAHU on 22 Jun 2022
thank you for clarification but how can we use loop so that we don't have to enter manually for large dataset?

Dyuman Joshi on 22 Jun 2022
%loop
psi=[0.0001,0.001,0.01];
r=[2.1891,2.194,2.2373];
for i=1:numel(psi)
W(i)=calculateW(psi(i),r(i));
end
W
W = 1×3
0.1602 0.1595 0.1534
%arrayfun
W1=arrayfun(@calculateW, psi, r)
W1 = 1×3
0.1602 0.1595 0.1534
function W = calculateW(psi,r)
alpha = (12*psi)^(1/3);
l1 = (r + alpha)*(r + alpha);
l1a = r*r - r*alpha + alpha*alpha;
L1 = log(l1/l1a);
t1 = (2*r - alpha) / (alpha*sqrt(3));
T1 = atan(t1);
l2 = (1 + alpha)*(1 + alpha);
l2a = (1*1) - (1*alpha) + (alpha*alpha);
L2 = log(l2/l2a);
t2 = (2*1 - alpha) / (alpha*sqrt(3));
T2 = atan(t2);
W1 = 4 * sqrt(3) * (T1 - T2) * log((r+alpha)/(1+alpha));
W2 = (L1 - L2)*log((r*r - r*alpha + alpha*alpha)/(1 - alpha + alpha*alpha));
W3 = 4*(T1 -T2)*(T1 -T2);
W4 = (L1 - L2) + 2 * sqrt(3) * (T1 - T2);
W5 = 3 ./ ((r - 1) * (r - 1));
W = ((W1 + W2 -W3) * W5)/W4;
end
AVINASH SAHU on 22 Jun 2022
Thank you so much