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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
\documentclass{article}
\usepackage[margin=0.3in]{geometry}
\usepackage{fontspec}
\usepackage{minted}
\usepackage{amsmath}
\setmainfont{FreeSans}
\setmonofont{FreeMono}
\begin{document}
\begin{minted}[linenos, frame=lines, fontsize=\scriptsize]{cpp}
#include <iostream>
#include "include/praktable.hpp"
using table = prak::table<double>;
using vecarg = const std::vector<f64> &;
using f64p = prak::pvalue<f64>;
table data;
const f64p g = {9.815710602, 0.001};
f64 get(std::string key) { return data[key, 0]; }
f64 getsqrt(vecarg v) { return std::sqrt(std::abs(get("x4_1") - v[0])); }
f64 a_A(vecarg a) { return 2/a[0]/a[0]; }
f64 J_mrga(vecarg a) { return a[0] * a[1] * a[1] * (a[2]/a[3] - 1); }
f64 ξ_x012(vecarg x) { return (x[1] - x[0]) / (2*x[2] - x[0] - x[1]); }
f64 Mfr_mgRξ(vecarg a) { return a[0] * a[1] * a[2] * a[3]; }
f64 J_mrgtx034(vecarg a) { return a[0]*a[1]*a[1]*(a[2]*a[3]*a[3] / 2 / std::pow(std::sqrt(a[4]-a[6]) - std::sqrt(a[4]-a[5]), 2) - 1); }
table ex1(std::string s) {
table ret(s);
const f64 dx03 = 0.04;
f64p m1 = {get("m1"), 0.00001},
R = {get("R2"), 0.00005},
x2 = {get("x2_1"), 0.01};
ret .add_column("x1", std::vector<f64>(ret.rows, NAN))
.add_column("t", std::vector<f64>(ret.rows, NAN))
.add_column("st", std::vector<f64>(ret.rows, NAN))
.add_column("sqrt", std::vector<f64>(ret.rows, NAN))
.add_column("x0", std::vector<f64>(ret.rows, NAN))
.add_column("ξ", std::vector<f64>(ret.rows, NAN))
.apply([dx03](vecarg a){return a[0]+dx03;}, {"x3"}, "x0")
.apply(prak::avg<f64>, {"x11", "x12", "x13"}, "x1")
.apply(prak::avg<f64>, {"t1", "t2", "t3"}, "t")
.apply(prak::stddev<f64>, {"t1", "t2", "t3"}, "st")
.apply(getsqrt, {"x0"}, "sqrt")
.apply([x2](vecarg a){return (a[0]-a[1])/((x2.val-a[0])+(x2.val-a[1]));}, {"x0", "x1"}, "ξ")
.delete_cols({"x11", "x12", "x13", "t1", "t2", "t3"});
f64p x1 = {ret.col_avg("x1"), 0.01 * std::sqrt(ret.rows)};
auto [A, B] = ret.least_squares_linear("sqrt", "t", "st", std::nullopt);
f64p a = prak::function<f64>(a_A, {A});
f64p J = prak::function<f64>(J_mrga, {m1, R, g, a});
f64p ξ = {ret.col_avg("ξ"), ret.col_stddev("ξ")};
f64p Mfr = prak::function<f64>(Mfr_mgRξ, {m1, g, R, ξ});
ret.write_plot("mnk.plot", "sqrt", "t", "st");
std::cout << ret << "\nA = " << A << "\nB = " << B
<< "\na = " << a << "\nJ = " << J
<< "\nξ = " << ξ << "\nМомент трения Mfr = " << Mfr
<< std::endl;
return ret;
}
void ex2(std::string data) {
table table(data);
f64p R0 = {get("R2_0"), 0.00005},
R1 = {get("R2_1"), 0.00005},
m0 = {get("m2_0"), 0.00001},
m1 = {get("m2_1"), 0.00001},
x0 = {get("x0_2"), 0.01},
x3 = {get("x3_2"), 0.01},
x4 = {get("x4_2"), 0.01};
table .add_column("t", std::vector<f64>(table.rows, NAN))
.add_column("st", std::vector<f64>(table.rows, NAN))
.add_column("Ji", std::vector<f64>(table.rows, NAN))
.add_column("sJi", std::vector<f64>(table.rows, NAN))
.add_column("0", std::vector<f64>(table.rows, 0.0))
.apply(prak::avg<f64>, {"t1", "t2", "t3"}, "t")
.apply(prak::stddev<f64>, {"t1", "t2", "t3"}, "st")
.delete_cols({"t1", "t2", "t3"})
;
for (size_t i = 0; i < table.rows; ++i) {
std::vector<f64p> args = {
table["M", i] == 0 ? m0 : m1,
table["R", i] == 0 ? R0 : R1,
g, {table["t", i], table["st", i]},
// я пр--бался и только на первом измерении у меня x4 = 4см, на остальных 10см
x0, x3, table["M", i] == 0 && table["R", i] == 0 ? f64p{0.04, 0.01} : x4,
};
auto [val, err] = prak::function<f64>(J_mrgtx034, args);
table["Ji", i] = val;
table["sJi", i] = err;
}
table.write_plot("ex2_1.plot", "Ji", "0", "sJi");
table.multiply_column("sJi", 1 / std::sqrt(1 - 0.75));
table.write_plot("ex2_2.plot", "Ji", "0", "sJi");
std::cout << table;
}
int main() {
data = table("common");
ex1("data1");
ex2("data2");
return 0;
}\end{minted}
\newpage
\begin{minted}[linenos, frame=lines]{gnuplot}
set term pngcairo size 1000, 800
set title "МНК t от (x_4 - x_0)^{1/2}"
set output "mnk.png"
f(x) = a1*x+b1
fit f(x) 'mnk.plot' using 1:2:3 yerror via a1, b1
set xlabel "(x_4 - x_0)^{1/2}, м^{1/2}"
set ylabel "t, c"
plot 'mnk.plot' using 1:2:3 with yerrorbars notitle lc 0 pt 1 lw 2,\
f(x) title "t" lc rgb "#a889e7" lw 2
reset
set output 'compar.png'
set title "Расположение J_i на оси J (пунктир: погрешность с коэффициентом доверия a=75%)"
set multiplot layout 4,1
set linetype 10 dashtype (5, 10)
unset ytics
unset ylabel
unset xlabel
set xrange[0.012:0.022]
set yrange[-1:1]
plot 'ex2_1.plot' every ::0::0 using 1:2:3 with xerrorbars title "J_1" lc 0 pt 1 lw 2,\
'ex2_2.plot' every ::0::0 using 1:2:3 with xerrorbars title "J_{1_a}" lt 10 lc 0 pt 1 lw 2
unset title
plot 'ex2_1.plot' every ::1::1 using 1:2:3 with xerrorbars title "J_2" lc 0 pt 1 lw 2,\
'ex2_2.plot' every ::1::1 using 1:2:3 with xerrorbars title "J_{2_a}" lt 10 lc 0 pt 1 lw 2
plot 'ex2_1.plot' every ::2::2 using 1:2:3 with xerrorbars title "J_3" lc 0 pt 1 lw 2,\
'ex2_2.plot' every ::2::2 using 1:2:3 with xerrorbars title "J_{3_a}" lt 10 lc 0 pt 1 lw 2
set xlabel "J, кг*м^2"
plot 'ex2_1.plot' every ::3::3 using 1:2:3 with xerrorbars title "J_4" lc 0 pt 1 lw 2,\
'ex2_2.plot' every ::3::3 using 1:2:3 with xerrorbars title "J_{4_a}" lt 10 lc 0 pt 1 lw 2
\end{minted}
\end{document}
|
