Код матлаба на C#
11.10.2021, 18:24. Показов 319. Ответов 0
Как правильно положить этот код на c#.
Matlab M | 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
| e=1.6021892*10^(-19);% заряд электрона [Кл] *
me=9.109543*10^(-31);% масса электрона [кг]*
c=2.99792458*10^(8);% скорость света в вакууме [м/с]
e0=8.85418782*10^(-12); % постоянная электрическая [Ф/м]
mag=12.5663706144*10^(-7); % постоянная магнитная [Гн/м]*
w=10^15:10^13:10^18
aem=1.6605655*10^(-27); % атомная еденица массы [кг]
r0=0.529177*10^(-10); % боровский радиус [м]
d=0.430; % оптимизированный экранирующий вклад
ro=3580; n1=1; n2=1;
z1=12; m1=24.305*aem; k1=3;%
z2=8; m2=15.9994*aem; k2=2;
n=ro/(n1*m1+n2*m2);
% Mg
%Слэйтор
z1s=z1-2*0.85-1*0.35;
z1p1=z1-2*0.85-3*0.35;
z1p2=z1-2*0.85-5*0.35;
z1p3=z1-2*0.85-7*0.35;
%Бор
r1s=r0*(k1-1)^2/z1s;
r1p1=r0*(k1-1)^2/z1p1;
r1p2=r0*(k1-1)^2/z1p2;
r1p3=r0*(k1-1)^2/z1p3;
w01s=(z1s*e^2/(4*pi*e0*me*r1s^3))^0.5;
w01p1=(z1p1*e^2/(4*pi*e0*me*r1p1^3))^0.5;
w01p2=(z1p2*e^2/(4*pi*e0*me*r1p2^3))^0.5;
w01p3=(z1p3*e^2/(4*pi*e0*me*r1p3^3))^0.5;
b1s=(mag*e^2*w01s^2)/(6*pi*me*c);
b1p1=(mag*e^2*w01p1^2)/(6*pi*me*c);
b1p2=(mag*e^2*w01p2^2)/(6*pi*me*c);
b1p3=(mag*e^2*w01p3^2)/(6*pi*me*c);
num1e=[2*e^2/me*n1];
den1s=[1 2*b1s w01s^2];
den1p1=[1 2*b1p1 w01p1^2];
den1p2=[1 2*b1p2 w01p2^2];
den1p3=[1 2*b1p3 w01p3^2];
[num1,den1]=parallel(num1e,den1s,num1e,den1p1);
[num1,den1]=parallel(num1,den1,num1e,den1p2);
[num1,den1]=parallel(num1,den1,num1e,den1p3);
% O
z2s=z2-2*1.00-1*d;
z2p1=z2-2*1.00-3*d;
z2p2=z2-2*1.00-5*d;
z2p3=z2-2*1.00-7*d;
r2s=r0*(k2+1)^2/z2s;
r2p1=r0*(k2+1)^2/z2p1;
r2p2=r0*(k2+1)^2/z2p2;
r2p3=r0*(k2+1)^2/z2p3;
w02s=(z2s*e^2/(4*pi*e0*me*r2s^3))^0.5;
w02p1=(z2p1*e^2/(4*pi*e0*me*r2p1^3))^0.5;
w02p2=(z2p2*e^2/(4*pi*e0*me*r2p2^3))^0.5;
w02p3=(z2p3*e^2/(4*pi*e0*me*r2p3^3))^0.5;
b2s=(mag*e^2*w02s^2)/(6*pi*me*c);
b2p1=(mag*e^2*w02p1^2)/(6*pi*me*c);
b2p2=(mag*e^2*w02p2^2)/(6*pi*me*c);
b2p3=(mag*e^2*w02p3^2)/(6*pi*me*c);
num2e=[2*e^2/me*n2];
den2s=[1 2*b2s w02s^2];
den2p1=[1 2*b2p1 w02p1^2];
den2p2=[1 2*b2p2 w02p2^2];
den2p3=[1 2*b2p3 w02p3^2];
[num2,den2]=parallel(num2e,den2s,num2e,den2p1);
[num2,den2]=parallel(num2,den2,num2e,den2p2);
[num2,den2]=parallel(num2,den2,num2e,den2p3);
[num,den]=parallel(num2,den2,num1,den1);
[re,im]=nyquist(num,den,w);
sre=1+2/(3*e0)*n*re;
semilogx(w,sre), axis([10^15 10^18 -5 15])
set(0,'DefaultAxesFontSize',14,'DefaultAxesFontName','Times New Roman');
title('Mg0');
xlabel('\omega (рад/с)');
ylabel('\epsilon`','Rotation', 0) |
|
Вот что у меня получилось, но график не совпадает с матлабовским. Что я упустил?
C# | 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
| // double l_min = 1e+15, l_max = 1e+18, h = 1e+13, x, y;
int ro = 3580, n1 = 1, n2 = 1;
int z1 = 12; double m1 = 24.305 * aem; int k1 = 3;
int z2 = 8; double m2 = 15.9994 * aem; int k2 = 2;
double n = ro / (n1 * m1 + n2 * m2);
//для ускорения расчетов
double alpha__const = 2 * qe * qe / me;
double omega0sqr_const = (qe * qe) / (4 * Math.PI * eps0 * me);
double betax2_const = (2 * qe * qe * mag) / (6 * Math.PI * c * me);
//Mg
//Слэйтeр
double z1s = z1 - 2 * 0.85 - 1 * 0.35;
double z1p1 = z1 - 2 * 0.85 - 3 * 0.35;
double z1p2 = z1 - 2 * 0.85 - 5 * 0.35;
double z1p3 = z1 - 2 * 0.85 - 7 * 0.35;
double[] slaytormassmg = new double[] { z1s, z1p1, z1p2, z1p3 };
//Бор
double[] bormassmg = new double[4];
double[] omegamg = new double[4];
double[] kofzatmg = new double[4];
double[] alfamg = new double[4];
double[] alfamg2 = new double[4];
for (int i = 0; i < bormassmg.Length; i++)
{
bormassmg[i] = r0 * Math.Pow(k1 - 1, 2) / slaytormassmg[i];
omegamg[i] = Math.Sqrt(slaytormassmg[i] * Math.Pow(qe, 2) / (4 * Math.PI * e0 * me * Math.Pow(bormassmg[i], 3)));
kofzatmg[i] = (mag * Math.Pow(qe, 2) * Math.Pow(omegamg[i], 2)) / (6 * Math.PI * me * c);
alfamg[i] += alpha__const * (omegamg[i] - (wa * wa)) /
(Math.Pow((omegamg[i] - (wa * wa)), 2.0) + Math.Pow(kofzatmg[i] * wa, 2.0)) * n1;
alfamg2[i] += alpha__const * kofzatmg[i] * wa /
(Math.Pow((omegamg[i] - (wa * wa)), 2.0) + Math.Pow(kofzatmg[i] * wa, 2.0)) * n1;
}
// O
double z2s = z2 - 2 * 1.00 - 1 * d;
double z2p1 = z2 - 2 * 1.00 - 3 * d;
double z2p2 = z2 - 2 * 1.00 - 5 * d;
double z2p3 = z2 - 2 * 1.00 - 7 * d;
double[] slaytormasso = new double[] { z2s, z2p1, z2p2, z2p3 };
//Бор
double[] bormasso = new double[4];
double[] omegao = new double[4];
double[] kofzato = new double[4];
double[] alfao = new double[4];
double[] alfao2 = new double[4];
for (int i = 0; i < bormasso.Length; i++)
{
bormasso[i] = r0 * Math.Pow(k2 - 1, 2) / slaytormasso[i];
omegao[i] = slaytormasso[i] * Math.Pow(qe, 2) / Math.Sqrt(4 * Math.PI * e0 * me * Math.Pow(bormasso[i], 3));
kofzato[i] = (mag * Math.Pow(qe, 2) * Math.Pow(omegao[i], 2)) / (6 * Math.PI * me * c);
alfao[i] += alpha__const * (omegao[i] - (wa * wa)) /
(Math.Pow((omegao[i] - (wa * wa)), 2.0) + Math.Pow(kofzato[i] * wa, 2.0)) * n2;
alfao2[i] += alpha__const * kofzato[i] * wa /
(Math.Pow((omegao[i] - (wa * wa)), 2.0) + Math.Pow(kofzato[i] * wa, 2.0)) * n2;
}
//Объединение
double[] alfasum = alfamg.Concat(alfao).ToArray();
double[] alfasum2 = alfamg2.Concat(alfao2).ToArray();
double[] omega = omegamg.Concat(omegao).ToArray();
double[] epsilon = new double[8];
double[] epsilon2 = new double[8];
double[] optic = new double[8];
for (int i = 0; i < alfasum.Length; i++)
{
epsilon[i] = 1.0 + 2.0 / (3.0 * eps0)*n * alfasum[i];
epsilon2[i] = 2.0 / (3.0 * eps0) *n* alfasum2[i];
optic[i] = Math.Sqrt((Math.Abs(epsilon[i]) + Math.Sqrt(Math.Pow(epsilon[i], 2.0) + Math.Pow(epsilon2[i], 2.0))) / 2.0);
}
for (int i = 0; i < omega.Length; i++)
{
//y = optic[i];
chart1.Series[0].Points.AddXY(omega[i], optic[i]);
chart1.ChartAreas[0].Axes[0].IsLogarithmic = true;
//chart1.ChartAreas[0].Axes[1].IsLogarithmic = true;
} |
|
0
|