I have code in matlab and Im trying "rewrite" it into Python. The problem is with different values of for loop. Unfortunately I can't find a solution.
Problem is with value in every array: nu theta mu M from i element (32 in matlab and 31 in python).
For both codes n=30
32 element in matlab = 31 element in python
In Python DTOR or RTOD are degrees to radians and vice versa.
-------------------------------------Matlab:-------------------------------------
for t = 1:2
i=n+2;
h = i;
k=0;
z=2;
for j= 1:n-1
while (i<=h+n-k-1)
if (i==h)
if (i==node_number-1)
y(i)=0;
theta(i)=0;
M(i) = Me;
nu(i) =PrandtlMeyer(M(i),gamma);
mu(i)= asin(1/M(i))*180/pi;
mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i))/2);
x(i)= (y(i-n+k)/mI)+x(i-n+k);
if (t==2)
plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
hold on;
end
else
%here is 32 element
y(i)=0;
theta(i)=0;
nu(i) = theta(i-n+k)+nu(i-n+k)+(1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))));
M(i) = InversePrandtlMeyer(nu(i));
mu(i)= asin(1/M(i))*180/pi;
mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i))/2);
x(i)= (y(i-n+k)/mI)+x(i-n+k);
%here ends 32 element
if (t==2)
plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
hold on;
end
end
else if (i==h+1)
if (i==node_number)
mI = tand((theta(i-n+k)+theta(i))/2);
mII =tand((mu(i-1)+theta(i-1)+theta(i)+mu(i))/2);
x(i) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
y(i) = y(i-1)+((x(i)-x(i-1))*mII);
y1(z) = y(i-1)+((x(i)-x(i-1))*mII);
x1(z) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
theta(i) = theta(i-1);
nu(i) = nu(i-1);
M(i) = M(i-1);
mu(i)= asin(1/M(i))*180/pi;
z=z+1;
if (t==2)
plot([x(i-1) x(i)],[y(i-1) y(i)]);
hold on;
plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
hold on;
end
else
%here is 33 element
mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i)-theta(i))/2);
mII = tand((mu(i)+theta(i)+mu(i-1))/2);
x(i) = (y(i-n+k)+(x(i-1)*mII)+(x(i-n+k)*mI))/(mII+mI);
y(i) = (x(i)-x(i-1))*mII;
nu(i) = (theta(i-n+k)+nu(i-n+k)+(nu(i-1)/2)+((1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))))*((y(i)-y(i-n+k))/y(i-n+k))))/1.5;
theta(i)= (nu(i)-nu(i-1))/2;
M(i) = InversePrandtlMeyer(nu(i));
mu(i)= asin(1/M(i))*180/pi;
%here end 33 element
if (t==2)
plot([x(i-1) x(i)],[y(i-1) y(i)]);
hold on;
plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
hold on;
end
end
else if (i==h+n-k-1)
mI = tand((theta(i-n+k)+theta(i))/2);
mII =tand((mu(i-1)+theta(i-1)+theta(i)+mu(i))/2);
x(i) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
y(i) = y(i-1)+((x(i)-x(i-1))*mII);
y1(z)=y(i-1)+((x(i)-x(i-1))*mII);
x1(z) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
theta(i) = theta(i-1);
nu(i) = nu(i-1);
M(i) = M(i-1);
mu(i)= asin(1/M(i))*180/pi;
z=z+1;
if (t==2)
plot([x(i-1) x(i)],[y(i-1) y(i)]);
hold on;
plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
hold on;
end
else
mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i)-theta(i))/2);
mII = tand((mu(i-1)+theta(i-1)+mu(i)+theta(i))/2);
x(i) = (y(i-n+k)-y(i-1)+(x(i-n+k)*mI)+(x(i-1)*mII))/(mII+mI);
y(i) = y(i-1)+((x(i)-x(i-1))*mII);
theta(i) = (theta(i-n+k)+nu(i-n+k)+theta(i-1)-nu(i-1)+((1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))))*((y(i)-y(i-n+k))/y(i-n+k)))-((1/(sqrt((M(i-1)^2)-1)-cotd(theta(i-1))))*((y(i)-y(i-1))/y(i-1))))/2;
nu(i) = (theta(i-n+k)+nu(i-n+k)-theta(i-1)+nu(i-1)+((1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))))*((y(i)-y(i-n+k))/y(i-n+k)))+((1/(sqrt((M(i-1)^2)-1)-cotd(theta(i-1))))*((y(i)-y(i-1))/y(i-1))))/2;
M(i) = InversePrandtlMeyer(nu(i));
mu(i)= asin(1/M(i))*180/pi;
if (t==2)
plot([x(i-1) x(i)],[y(i-1) y(i)]);
hold on;
plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
hold on;
end
end
end
end
i = i+1;
end
k=k+1;
h=i;
end
end
And first few rows of Matlab "theta" values:
0.4397 0.8793 1.3190 1.7587 2.1983 2.6380 3.0776 3.5173
3.9570 4.3966 4.8363 5.2760 5.7156 6.1553 6.5949 7.0346
7.4743 7.9139 8.3536 8.7933 9.2329 9.6726 10.1122 10.5519
10.9916 11.4312 11.8709 12.3106 12.7502 13.1899 13.1899 0
0.3032 0.7453 1.1877 1.6301 2.0725 2.5149 2.9572 3.3995
3.8418 4.2840 4.7262 5.1684 5.6106 6.0528 6.4951 6.9373
7.3795 7.8218 8.2641 8.7065 9.1488 9.5913 10.0337 10.4762
10.9188 11.3614 11.8041 12.2469 12.2469 0 0.2977 0.7407
1.1839 1.6271 2.0702 2.5133 2.9562 3.3991 3.8420 4.2848
4.7276 5.1704 5.6131 6.0559 6.4986 6.9414 7.3842 7.8270
-------------------------------------Python:-------------------------------------
i = n + 1
h = i
k = 0
z = 1
for j in range(n-1):
while i <= h + n - k-1:
if i == h:
if i == node_number -1:
y[i] = 0
theta[i] = 0
M[i] = Me
nu[i] = PM(M[i], gamma)
mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
mI = np.tan((mu[i - n + k] - theta[i - n + k] + mu[i]) / 2 * DTOR)
x[i] = (y[i - n + k] / mI) + x[i - n + k]
else:
#here is 31 element
y[i] = 0
theta[i] = 0
nu[i] = theta[i - n + k] + nu[i - n + k] + (1 / (np.sqrt((M[i - n + k]**2) - 1) - np.cos((theta[i - n + k])*DTOR)))
M[i] = InversePrandtlMeyer(nu[i])
mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
mI = np.tan(((mu[i - n + k] - theta[i - n + k] + mu[i]) / 2) * DTOR)
x[i] = (y[i - n + k] / mI) + x[i - n + k]
#here ends 31 element
elif i == h + 1:
if i == node_number:
mI = np.tan((theta[i - n + k] + theta[i]) / 2 * DTOR)
mII = np.tan((mu[i - 1] + theta[i - 1] + theta[i] + mu[i]) / 2 * DTOR)
x[i] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
y[i] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
y1[z] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
x1[z] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
theta[i] = theta[i - 1]
nu[i] = nu[i - 1]
M[i] = M[i - 1]
mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
z = z + 1
else:
#here is 32 element
mI = np.tan((mu[i - n + k] - theta[i - n + k] + mu[i] - theta[i]) / 2 * DTOR)
mII = np.tan((mu[i] + theta[i] + mu[i - 1]) / 2 * DTOR)
x[i] = (y[i - n + k] + (x[i - 1] * mII) + (x[i - n + k] * mI)) / (mII + mI)
y[i] = (x[i] - x[i - 1]) * mII
nu[i] = (theta[i - n + k] + nu[i - n + k] + (nu[i - 1] / 2) + ((1 / (np.sqrt((M[i - n + k] ** 2) - 1) - np.cos((theta[i - n + k])*DTOR))) * ((y[i] - y[i - n + k]) / y[i - n + k]))) / 1.5
theta[i] = (nu[i] - nu[i - 1]) / 2
M[i] = InversePrandtlMeyer(nu[i])
mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
#here end 32 element
elif i == h + n - k-1:
mI = np.tan((theta[i - n + k] + theta[i]) / 2 * DTOR)
mII = np.tan((mu[i - 1] + theta[i - 1] + theta[i] + mu[i]) / 2 * DTOR)
x[i] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
y[i] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
y1[z] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
x1[z] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
theta[i] = theta[i - 1]
nu[i] = nu[i - 1]
M[i] = M[i - 1]
mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
z = z + 1
else:
mI = np.tan((mu[i - n + k] - theta[i - n + k] + mu[i] - theta[i]) / 2 * DTOR)
mII = np.tan((mu[i - 1] + theta[i - 1] + mu[i] + theta[i]) / 2 * DTOR)
x[i] = (y[i - n + k] - y[i - 1] + (x[i - n + k] * mI) + (x[i - 1] * mII)) / (mII + mI)
y[i] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
theta[i] = (theta[i - n + k] + nu[i - n + k] + theta[i - 1] - nu[i - 1] + ((1 / (np.sqrt((M[i - n + k] ** 2) - 1) - np.cos((theta[i - n + k]) * DTOR))) * ((y[i] - y[i - n + k]) / y[i - n + k])) - ((1 / (np.sqrt((M[i - 1] ** 2) - 1) - np.cos((theta[i - 1])*DTOR))) * ((y[i] - y[i - 1]) / y[i - 1]))) / 2
nu[i] = (theta[i - n + k] + nu[i - n + k] - theta[i - 1] + nu[i - 1] + ((1 / (np.sqrt((M[i - n + k] ** 2) - 1) - np.cos((theta[i - n + k])*DTOR))) * ((y[i] - y[i - n + k]) / y[i - n + k])) + ((1 / (np.sqrt((M[i - 1] ** 2) - 1) - np.cos((theta[i - 1])*DTOR))) * ((y[i] - y[i - 1]) / y[i - 1]))) / 2
M[i] = InversePrandtlMeyer(nu[i])
mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
i = i + 1
k = k + 1
h = i
And first few rows of Python "theta" values:
[ 0.43966268 0.87932536 1.31898804 1.75865072 2.1983134 2.63797608
3.07763876 3.51730144 3.95696412 4.3966268 4.83628948 5.27595216
5.71561484 6.15527752 6.5949402 7.03460288 7.47426556 7.91392824
8.35359092 8.7932536 9.23291628 9.67257896 10.11224165 10.55190433
10.99156701 11.43122969 11.87089237 12.31055505 12.75021773 13.18988041
13.18988041 0. 1.12391865 2.04904177 2.75249233 3.37815396
3.96686921 4.53558282 5.09310744 5.64490943 6.19500187 6.74690367
7.30431746 7.87186621 8.45626406 9.06871626 9.73089387 10.49343203
11.51563061 13.74816848 0.85949953 nan nan nan
nan nan nan nan nan nan
nan 0. 1.48303468 2.8885035 3.74795106 4.46269741
5.11686874 5.74107193 6.35015955 6.95312861 7.55674416 8.16737947
8.79250326 9.44272463 10.13594268 10.90840157 11.85472483 13.39395428
13.37483292 nan nan
More answer...
2019-11-25, 05:20
Matlab for loop, different values in Python
平板模式
