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