1 /*
2 platerev.c:
3
4 Copyright (C) 2006 by Stefan Bilbao
5 2012 John ffitch
6
7 This file is part of Csound.
8
9 The Csound Library is free software; you can redistribute it
10 and/or modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
13
14 Csound is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public
20 License along with Csound; if not, write to the Free Software
21 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
22 02110-1301 USA
23 */
24
25 #include "csdl.h"
26 #include <math.h>
27
28 /* #undef CS_KSMPS */
29 /* #define CS_KSMPS (csound->GetKsmps(csound)) */
30
31 typedef struct {
32 OPDS h;
33 MYFLT *aout[40];
34 MYFLT *tabins;
35 MYFLT *tabout;
36 MYFLT *bndry;
37 MYFLT *asp;
38 MYFLT *stiff;
39 MYFLT *decay;
40 MYFLT *loss;
41 MYFLT *ain[40];
42 // Internals
43 double s00, s10, s01, s11, s20, s02, t00, t01, t10;
44 uint32_t nin, nout, Nx, Ny;
45 double *u, *u1, *u2;
46 AUXCH auxch;
47 double L, dy, dt;
48 MYFLT *in_param, *out_param;
49 double ci[40], si[40], co[40], so[40];
50 } PLATE;
51
52
platerev_init(CSOUND * csound,PLATE * p)53 static int32_t platerev_init(CSOUND *csound, PLATE *p)
54 {
55 FUNC *inp, *outp;
56 double a = *p->asp;
57 double dt = (p->dt = 1.0/csound->GetSr(csound)); /* time step */
58 double sig = (csound->GetSr(csound)+csound->GetSr(csound))*
59 (POWER(10.0, FL(3.0)*dt/(*p->decay))-FL(1.0)); /* loss constant */
60 double b2 = *p->loss;
61 double dxmin = 2.0*sqrt(dt*(b2+hypot(*p->loss, *p->stiff)));
62 uint32_t Nx = (p->Nx = (uint32_t)floor(1.0/dxmin));
63 uint32_t Nx5 = Nx+5;
64 double dx = 1.0/(double)Nx;
65 uint32_t Ny = (p->Ny = (uint32_t)floor(a*Nx));
66 double dy = (p->dy = *p->asp/Ny);
67 double alf = dx/dy;
68 double mu = dt*(*p->stiff)*Nx*Nx;
69 double mu2 = mu*mu;
70 double eta = 1.0/(1.0+sig*dt);
71 double V = 2.0*b2*dt*Nx*Nx;
72 uint32_t qq;
73
74 p->nin = (int32_t) (p->INOCOUNT) - 7; p->nout = (int32_t) (p->OUTOCOUNT);
75 if (UNLIKELY((inp = csound->FTnp2Find(csound,p->tabins)) == NULL ||
76 inp->flen < (uint32_t)3*p->nin)) {
77 return csound->InitError(csound, "%s",
78 Str("Missing input table or too short"));
79 }
80 if (UNLIKELY((outp = csound->FTnp2Find(csound,p->tabout)) == NULL ||
81 outp->flen < (uint32_t)3*p->nout)) {
82 return csound->InitError(csound, "%s",
83 Str("Missing output table or too short"));
84 }
85 p->in_param = inp->ftable;
86 p->out_param = outp->ftable;
87 p->L = (a<1.0 ? a : 1.0);
88 csound->AuxAlloc(csound, 3*Nx5*(Ny+5)*sizeof(double), &p->auxch);
89 p->u = (double*)p->auxch.auxp;
90 p->u1 = p->u+Nx5*(Ny+5);
91 p->u2 = p->u1+Nx5*(Ny+5);
92 p->s00 = 2.0*eta*(1.0-mu2*(3.0+4.0*alf*alf+3.0*alf*alf*alf*alf)-
93 V*(1.0+alf*alf));
94 p->s10 = (4.0*mu2*(1.0+alf*alf)+V)*eta;
95 p->s01 = alf*alf*(4.0*mu2*(1.0+alf*alf)+V)*eta;
96 p->s11 = -2.0*mu2*eta*alf*alf;
97 p->s02 = (p->s20 = -eta*mu2)*alf*alf*alf*alf;
98 p->t00 = (-(1.0-sig*dt)+2.0*V*(1.0+alf*alf))*eta;
99 p->t10 = -V*eta;
100 p->t01 = -V*eta*alf*alf;
101 for (qq=0; qq<p->nin; qq++) {
102 p->ci[qq] = cos((double)p->in_param[3*qq+2]);
103 p->si[qq] = sin((double)p->in_param[3*qq+2]);
104 }
105 for (qq=0; qq<p->nout; qq++) {
106 p->co[qq] = cos((double)p->out_param[3*qq+2]);
107 p->so[qq] = sin((double)p->out_param[3*qq+2]);
108 }
109
110 return OK;
111 }
112
113
platerev(CSOUND * csound,PLATE * p)114 static int32_t platerev(CSOUND *csound, PLATE *p)
115 {
116 IGN(csound);
117 uint32_t offset = p->h.insdshead->ksmps_offset;
118 uint32_t early = p->h.insdshead->ksmps_no_end;
119 uint32_t i, j, nsmps = CS_KSMPS;
120 uint32_t Ny = p->Ny, Nx = p->Nx;
121 uint32_t Nx5 = Nx+5;
122 int32_t bc = (int32_t) MYFLT2LONG(*p->bndry);
123 double *u = p->u, *u1 = p->u1, *u2 = p->u2;
124 double s00 = p->s00, s10 = p->s10, s01 = p->s01,
125 s11 = p->s11, s20 = p->s20, s02 = p->s02,
126 t00 = p->t00, t10 = p->t10, t01 = p->t01;
127 double dt = p->dt, dy = p->dy;
128 uint32_t n, qq;
129 MYFLT *uin;
130 double wi[40], wo[40], sdi[40], cdi[40], sdo[40], cdo[40];
131
132 if (UNLIKELY(early)) nsmps -= early;
133 for (qq=0; qq<(uint32_t)p->nin; qq++) {
134 double delta = TWOPI*(double)p->in_param[3*qq]*dt;
135 cdi[qq] = cos(delta);
136 sdi[qq] = sin(delta);
137 wi[qq] = p->L*0.5*(double)p->in_param[3*qq+1];
138 delta = TWOPI*(double)p->out_param[3*qq]*dt;
139 cdo[qq] = cos(delta);
140 sdo[qq] = sin(delta);
141 wo[qq] = (p->L*0.5)*(double)p->out_param[3*qq+1];
142 if (UNLIKELY(offset)) memset(p->aout[qq], '\0', offset*sizeof(MYFLT));
143 if (UNLIKELY(early)) memset(&p->aout[qq][nsmps], '\0', early*sizeof(MYFLT));
144 }
145 for (n=offset; n<nsmps; n++) {
146 /* interior grid points*/
147 /*u = conv2(u1,S,'same')+conv2(u2,T,'same'); */
148 for (j=2; j<Ny+3; j++) /* Loop from 2,3,...Nf, Nf+1, Nf+2 */
149 for (i=2; i<Nx+3; i++) {
150 int32_t ij = i+Nx5*j;
151 u[ij] = s00*u1[ij]+
152 s10*(u1[ij-Nx5]+u1[ij+Nx5])+
153 s20*(u1[ij-2*Nx5]+u1[ij+2*Nx5])+
154 s01*(u1[ij-1]+u1[ij+1])+
155 s02*(u1[ij-2]+u1[ij+2])+
156 s11*(u1[ij-1-Nx5]+u1[ij+1+Nx5]+
157 u1[ij-1+Nx5]+u1[ij+1-Nx5]);
158 u[ij] += t00*u2[ij]+
159 t10*(u2[ij-Nx5]+u2[ij+Nx5])+
160 t01*(u2[ij-1]+u2[ij+1]);
161 }
162 /* boundary grid points*/
163 if (bc==1) { /* clamped*/
164 int32_t jj = Nx5*j;
165 for (j=0; j<Ny+5; j++) {
166 u[0+jj] = u[2+jj] = u[Nx+2+jj] = u[Nx+4+jj] = 0.0;
167 }
168 for (j=2; j<Ny+3; j++) {
169 u[1+jj] = u[3+jj];
170 u[Nx+3+jj] = u[Nx+1+jj];
171 }
172 for (i=0; i<Nx+5; i++) {
173 u[i+Nx5*0] = u[i+Nx5*2] = u[i+Nx5*(Ny+2)] = u[i+Nx5*(Ny+4)] = 0.0;
174 }
175 for (i=2; i<Nx+3; i++) {
176 u[i+Nx5*1] = u[i+Nx5*3];
177 u[i+Nx5*(Ny+3)] = u[i+Nx5*(Ny+1)];
178 }
179 u[1+Nx5*1] = u[1+Nx5*(Ny+3)] = u[Nx+3+Nx5*1] = u[Nx+3+Nx5*(Ny+3)] = 0.0;
180 }
181 else if (bc==2) { /* pivoting*/
182 int32_t jj = Nx5*j;
183 for (j=0; j<Ny+5; j++) {
184 u[0+jj] = u[2+jj] = u[Nx+2+jj] = u[Nx+4+jj] = 0.0;
185 }
186 for (j=2; j<Ny+3; j++) {
187 u[1+jj] = -u[3+jj];
188 u[Nx+3+jj] = -u[Nx+1+jj];
189 }
190 for (i=0; i<Nx+5; i++) {
191 u[i+Nx5*0] = u[i+Nx5*2] = u[i+Nx5*(Ny+2)] = u[i+Nx5*(Ny+4)] = 0.0;
192 }
193 for (i=2; i<Nx+3; i++) {
194 u[i+Nx5*1] = -u[i+Nx5*3];
195 u[i+Nx5*(Ny+3)] = -u[i+Nx5*(Ny+1)];
196 }
197 }
198 /* else completely free */
199 /* insert excitation*/
200 for (qq=0; qq<p->nin; qq++) {
201 double w = wi[qq];
202 double cv = p->ci[qq]*cdi[qq] - p->si[qq]*sdi[qq];
203 double sv = p->ci[qq]*sdi[qq] + p->si[qq]*cdi[qq];
204 double xid = (0.5+w*cv)*Nx;
205 double yid = ((*p->asp)*0.5+w*sv)/dy;
206 int32_t xi = (int32_t)(floor(xid))+2;
207 int32_t yi = (int32_t)(floor(yid))+2;
208 double xf = xid-(double)(xi-2);
209 double yf = yid-(double)(yi-2);
210 double xyf = xf*yf;
211 uin=p->ain[qq];
212 p->ci[qq] = cv; p->si[qq] = sv;
213 yi = Nx5*yi + xi;
214 u[yi] += (1.0-xf-yf+xyf)*uin[n];
215 u[1+yi] += (xf-xyf)*uin[n];
216 u[1+Nx5+yi] += xyf*uin[n];
217 u[Nx5+yi] += (yf-xyf)*uin[n];
218 }
219 /* %%%% readout */
220 for (qq=0; qq<p->nout; qq++) {
221 double w = wo[qq];
222 double cv = p->co[qq]*cdo[qq] - p->so[qq]*sdo[qq];
223 double sv = p->co[qq]*sdo[qq] + p->so[qq]*cdo[qq];
224 double xod = (0.5+w*cv)*Nx;
225 double yod = (*p->asp*0.5+w*sv)/dy;
226 int32_t xo = (int32_t)(floor(xod))+2;
227 int32_t yo = (int32_t)(floor(yod))+2;
228 double xf = xod-(double)(xo-2);
229 double yf = yod-(double)(yo-2);
230 double xyf = xf*yf;
231 p->co[qq] = cv; p->so[qq] = sv;
232 yo = yo*Nx5 + xo;
233 (p->aout[qq])[n] = (MYFLT)((1.0-xf-yf+xyf)*u[yo]+
234 (yf-xyf)*u[Nx5+yo]+
235 (xf-xyf)*u[1+yo]+
236 xyf*u[1+Nx5+yo])/FL(25.0);
237 }
238 {
239 double *tmp = u2; /* cycle U*/
240 u2 = u1;
241 u1 = u;
242 u = tmp;
243 }
244 }
245 p->u = u; p->u1 = u1; p->u2 = u2;
246 return OK;
247 }
248
249 static OENTRY localops[] =
250 {
251 { "platerev", sizeof(PLATE), 0, 3, "mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm",
252 "iikiiiiy",
253 (SUBR) platerev_init, (SUBR) platerev
254 },
255 };
256
257 LINKAGE
258