1 #include <nvector/nvector_serial.h>
2 
3 /* Accessor macro */
4 #define Ith(v,i)    NV_Ith_S(v,i)
5 
6 const int fs = 96000;
7 //const int fs = 8*44100;
8 
9 struct Dim {
10     double lower, upper;
11     int size, stride;
pointDim12     inline double point(int i) { return lower + ((upper-lower) / (size-1)) * i; }
pointDim13     inline double point(double i) { return lower + ((upper-lower) / (size-1)) * i; }
indexDim14     inline void index(double v, int& i, double& w) {
15 	w = (v-lower) / (upper-lower) * (size-1);
16 	i = static_cast<int>(w);
17 	if (i < 0) {
18 	    i = 0;
19 	    w = 0.0;
20 	} else if (i > size-2) {
21 	    i = size-2;
22 	    w = 1.0;
23 	} else {
24 	    w -= i;
25 	}
26     }
27 };
28 
29 
30 class Grid {
31 protected:
32     float *grid;		// ndim-dimensional cube
33     int vals;			// size of cube cell
34     int ndim;			// number of dimensions
35     Dim *dim;			// dimension descriptor
36     double *weight;		// coordinates inside grid cell
37     int *idx;			// index of base point
38     int *idx2;			// running index of neighbor nodes
39     int p2dim;			// 2^ndim
40     double *pts;		// size = p2dim*vals
41     inline float *cell(int *ix);
42     int calc_strides();
43 public:
44     Grid(int d, Dim *v, int n);
45     ~Grid();
46 };
47 
48 
49 class ComponentBase;
50 
51 class Cube: public Grid {
52 private:
53     ComponentBase& cb;
54     int mmap_size;
55     int n_in;                   // length of input vector
56     int fill(int n, int *ix, N_Vector x, N_Vector u0);
57 public:
58     Cube(int d, Dim *v, int n, int *ix, int n_in_, ComponentBase& cb, N_Vector u0, const char* fname=0);
59     ~Cube();
60     void startvalue(N_Vector x, N_Vector s, N_Vector u);
61     int calc(N_Vector x, N_Vector s, N_Vector u);
62 };
63 
64 
65 class ComponentBase {
66 public:
67     int NEQ;
68     int NDIM;
69     int NVALS;
70     int N_IN;
71     int N_OUT;
72     int n_params;
73     realtype *params;
74     const char **param_names;
75     const char **in_names;
76     const char **out_names;
77     const char **state_names;
78     const char **var_names;
79     int *ix;
80     bool verbose;
81     struct UserData {
82 	realtype *inval;
83 	realtype *state;
84     };
85 private:
86     Dim *ranges;
87     Cube *cube;
88     UserData user_data;
89     static int sfunc(N_Vector x, N_Vector u, void *data);
90 protected:
91     N_Vector constraints;
92     realtype& Gco;
93     realtype& Gcf;
94     realtype& mu;
95     realtype& Ex;
96     realtype& Kg1;
97     realtype& Kg2;
98     realtype& Kp;
99     realtype& Kvb;
100     inline realtype Ig(realtype Ugk);
101     inline realtype Ia(realtype Ugk, realtype Uak);
102     inline realtype Iap(realtype Ug1k, realtype Ug2k, realtype Uak);
103     inline realtype Is(realtype Ug1k, realtype Ug2k);
104     static const int param_off = 8;
set_names(const char ** pn,const char ** p,int n)105     inline void set_names(const char **pn, const char **p, int n) {
106 	while (*p) { *pn++ = *p++; n--; }
107 	assert(n == 0);
108     }
109 public:
110     virtual int func(N_Vector x, N_Vector u, UserData *data) = 0;
111     virtual void update(N_Vector y, N_Vector x, N_Vector state) = 0;
112     int findzero(realtype *x, realtype *s, N_Vector u);
113     int set_range(int i, double lower, double upper, int size);
114     void init(N_Vector u0, const char* fname=0);
115     int calc(N_Vector x, N_Vector s, N_Vector u);
116     void startvalue(N_Vector x, N_Vector s, N_Vector u);
117     ComponentBase(int neq, int ndim, int nvals, int n_in, int n_out, int n_params_);
118     virtual ~ComponentBase();
119 };
120 
121 
122 class TriodeCircuit: public ComponentBase {
123 private:
124     realtype& Ck;
125 
126     realtype& G1; // 1/R1
127     realtype& G2; // 1/R2
128     realtype& Gg; // 1/Rg
129     realtype& Gk; // 1/Rk
130     realtype& Ga; // 1/Ra
131     realtype& Gl; // 1/Rl
132 
133     realtype& Un;
134 
135     virtual int func(N_Vector u, N_Vector f, UserData *user_data);
136     virtual void update(N_Vector y, N_Vector x, N_Vector state);
137 public:
138     TriodeCircuit();
139 };
140 
141 
142 class CoupledTriodeCircuit: public ComponentBase {
143 private:
144     realtype& Ck;
145     realtype& Ca;
146     realtype& Un;
147     realtype& G1;
148     realtype& G2;
149     realtype& Gg;
150     realtype& Gk;
151     realtype& Ga;
152     realtype& G3;
153     realtype& Gg2;
154     realtype& Gk2;
155     realtype& Ga2;
156     realtype& Gl;
157 
158     virtual int func(N_Vector u, N_Vector f, UserData *user_data);
159     virtual void update(N_Vector y, N_Vector x, N_Vector state);
160 public:
161     CoupledTriodeCircuit();
162 };
163 
164 
165 class PowerAmpGate: public ComponentBase {
166 private:
167     realtype& C1;
168     realtype& G1; // 1/R1
169     realtype& Gb; // 1/Rb
170     realtype& Gg; // 1/Rg
171     realtype& Ub;
172 
173     virtual int func(N_Vector u, N_Vector f, UserData *user_data);
174     virtual void update(N_Vector y, N_Vector x, N_Vector state);
175 public:
176     PowerAmpGate();
177 };
178 
179 
180 class PowerAmpPlate: public ComponentBase {
181 private:
182     realtype& C2;
183     realtype& Gd; // 1/Rd
184     realtype& Ga; // 1/Ra
185     realtype& Gs; // 1/Rs
186     realtype& Un;
187 
188     virtual int func(N_Vector u, N_Vector f, UserData *user_data);
189     virtual void update(N_Vector y, N_Vector x, N_Vector state);
190 public:
191     PowerAmpPlate();
192 };
193 
194 
195 class PhaseSplitter: public ComponentBase {
196 private:
197     realtype& C1;
198     realtype& C2;
199     realtype& C3;
200     realtype& Un;
201     realtype& G1;
202     realtype& Gg1;
203     realtype& Gg2;
204     realtype& Gk;
205     realtype& G2;
206     realtype& G3;
207     realtype& G4;
208     realtype& G5;
209     realtype& Ga1;
210     realtype& Ga2;
211     realtype& Gl1;
212     realtype& Gl2;
213 
214     virtual int func(N_Vector u, N_Vector f, UserData *user_data);
215     virtual void update(N_Vector y, N_Vector x, N_Vector state);
216 public:
217     PhaseSplitter();
feedback_coeff(realtype * b0,realtype * b1,realtype * a1,realtype pres)218     inline void feedback_coeff(realtype *b0, realtype *b1, realtype *a1, realtype pres) {
219 	realtype R4v = pres/G4;
220 	realtype B0 = 1/G5;
221 	realtype B1 = R4v/G5*C3;
222 	realtype A0 = 1/G3 + 1/G5;
223 	realtype A1 = C3*R4v/G5 + C3/(G3*G5) + C3/G3*R4v;
224 	realtype a = A0 + 2*A1*fs;
225 	*a1 = -(A0 - 2*A1*fs) / a;
226 	*b0 = (B0 + 2*B1*fs) / a;
227 	*b1 = (B0 - 2*B1*fs) / a;
228     }
229 };
230 
231 
232 class PowerAmpPlateTrans: public ComponentBase {
233 private:
234     realtype& C2;
235     realtype& Gd; // 1/Rd
236     realtype& Ga; // 1/Ra
237     realtype& Gs; // 1/Rs
238     realtype& Un;
239 
240     virtual int func(N_Vector u, N_Vector f, UserData *user_data);
241     virtual void update(N_Vector y, N_Vector x, N_Vector state);
242 public:
243     PowerAmpPlateTrans();
244 };
245 
246 
247 class Spline {
248 private:
249     int count;
250     double x0, step;
251     double *c;
252     double *f;
253 public:
254     Spline(double *y, int n, double xstart, double h);
255     ~Spline();
256     double eval(double x);
257 };
258