Compass/LCR_Oscillations.hpp at main · thecompassproject/Compass · GitHub

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
152
153
154
155
156
157
158
159
class LCR
{
private:
    double _L, _C, _R, _h, _method_choice;
    double _initialTime, _finalTime, _I, _Q;
    int _n;

public:
    string input(double method_choice, double L, double C, double R)
    {
        _L = 1e-3 * L;
        _C = 1e-6 * C;
        _R = R;
        _method_choice = method_choice;

        if (pow(_R, 2) > (4 * _L) / _C)
        {
            // std::cout << "\nThe oscillator is over-damped.";
            return "Message: The oscillator is over-damped.";
        }

        else if (pow(_R, 2) < (4 * _L) / _C)
        {
            return "Message: The oscillator is under-damped.";
        }

        else
        {
            return "Message: The oscillator is critically-damped.";
        }
    }

    // Wrapper function to match the required signature
    double f_wrapper(double A, double B, double Q1, double I1)
    {
        return -(((Q1 / _C) + (_R * I1)) / _L);
    }

    double cha_wrapper(double A, double B, double Q1, double I1)
    {
        return I1;
    }

    void kart(double initialTime, double finalTime, double I, double Q, int n)
    {
        _finalTime = 1e-6 * finalTime;
        _I = 1e-3 * I;
        _Q = 1e-3 * Q;
        _n = n;
        _initialTime = 0.0;
        double h = (_finalTime - _initialTime) / _n;

        char *file_name = "lcr.dat";
        if (_method_choice == 1.0)
        {
            file_name = "lcr_rk4.dat";
        }
        else if (_method_choice == 2.0)
        {
            file_name = "lcr_euler.dat";
        }
        else if (_method_choice == 3.0)
        {
            file_name = "lcr_rals.dat";
        }

        std::ofstream out(file_name);

        // A and B are dummy variables for the sake of compatibility
        double A = 0;
        double B = 0;

        for (double t = _initialTime; t <= _finalTime; t += h)
        {
            // steps calculating is replaced by function call from included header file
            if (_method_choice == 1.0)
            {
                calc1(std::bind(&LCR::cha_wrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4),
                      std::bind(&LCR::f_wrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4),
                      A, B, _Q, _I, h);
            }

            else if (_method_choice == 2.0)
            {
                calc4(std::bind(&LCR::cha_wrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4),
                      std::bind(&LCR::f_wrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4),
                      A, B, _Q, _I, h);
            }

            else if (_method_choice == 3.0)
            {
                calc(std::bind(&LCR::cha_wrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4),
                     std::bind(&LCR::f_wrapper, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4),
                     A, B, _Q, _I, h);
            }

            else
            {
                std::cout << "ERROR !!Not Any Method!!\n";
                break;
            }

            out << t << "  " << _I << "  " << _Q << "\n";
        }
        out.close();

        // Important if building for both Unix and Windows OS
        // Plot the data using gnuplot
#ifdef _WIN32
        FILE *gnuplotPipe = _popen("gnuplot -persistent", "w");
#else
        FILE *gnuplotPipe = popen("gnuplot -persistent", "w");
#endif

        if (!gnuplotPipe)
        {
            std::cerr << "Error opening pipe to gnuplot" << std::endl;
            return;
        }

        // Create plot command for gnuplot
        std::string s;
        const char *ss = s.append("plot '").append(file_name).append("' using 1:2 with lines title 'Current'\n").c_str();

        if (pow(_R, 2) > (4 * _L) / _C)
        {

            fprintf(gnuplotPipe, "set term wxt\n");
            fprintf(gnuplotPipe, "set title 'Over Damped'\n");
            fprintf(gnuplotPipe, "set xlabel 'Time'\n");
            fprintf(gnuplotPipe, "set ylabel 'Current'\n");
            fprintf(gnuplotPipe, ss);
        }
        else if (pow(_R, 2) < (4 * _L) / _C)
        {
            fprintf(gnuplotPipe, "set term wxt\n");
            fprintf(gnuplotPipe, "set title 'Under Damped'\n");
            fprintf(gnuplotPipe, "set xlabel 'Time'\n");
            fprintf(gnuplotPipe, "set ylabel 'Current'\n");
            fprintf(gnuplotPipe, ss);
        }
        else
        {
            fprintf(gnuplotPipe, "set term wxt\n");
            fprintf(gnuplotPipe, "set title 'Critically Damped'\n");
            fprintf(gnuplotPipe, "set xlabel 'Time'\n");
            fprintf(gnuplotPipe, "set ylabel 'Current'\n");
            fprintf(gnuplotPipe, ss);
        }

#ifdef _WIN32
        _pclose(gnuplotPipe);
#else
        pclose(gnuplotPipe);
#endif

        return;
    }
};