1 /*
2 wterrain2.c:
3
4 Copyright (C) 2002 Matt Gilliard, John ffitch
5 for the original file wave-terrain.c from the csound distribution
6
7 Modifications and enhancements by (C) 2020 Christian Bacher
8
9 This file is part of Csound.
10
11 The Csound Library is free software; you can redistribute it
12 and/or modify it under the terms of the GNU Lesser General Public
13 License as published by the Free Software Foundation; either
14 version 2.1 of the License, or (at your option) any later version.
15
16 Csound is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 GNU Lesser General Public License for more details.
20
21 You should have received a copy of the GNU Lesser General Public
22 License along with Csound; if not, write to the Free Software
23 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
24 02110-1301 USA
25 */
26
27 #include <csdl.h>
28 #include <math.h>
29
30 /* Wave-terrain synthesis opcode
31 *
32 * author: m gilliard
33 * en6mjg@bath.ac.uk
34 *
35 * enhancements and modifications
36 * Christian Bacher docb22@googlemail.com
37 * Changes to the original:
38 * - Added curves: lemniskate (G), limacon with parameter, cornoid with parameter, trisec (Ceva) with parameter, scarabeus with 2 parameters, folium with parameter
39 * - tables are krate
40 * - added k parameter for rotating the curve arround the current x,y
41 */
42 typedef struct {
43
44 OPDS h;
45
46 MYFLT *aout;
47 MYFLT *kamp;
48 MYFLT *kpch;
49 MYFLT *kx;
50 MYFLT *ky;
51 MYFLT *krx;
52 MYFLT *kry;
53 MYFLT *krot; // rotation of the curve
54 MYFLT *ktabx, *ktaby; /* Table numbers */
55 MYFLT *kfunc; // the curve index
56 MYFLT *kparam;
57
58 /* Internals */
59 MYFLT oldfnx; // storage of the current table for k-rate table change
60 MYFLT oldfny; // storage of the current table for k-rate table change
61
62 MYFLT *xarr, *yarr; /* Actual tables */
63
64 MYFLT sizx, sizy;
65 double theta;
66
67 } WAVETER;
68
rotate_point(MYFLT cx,MYFLT cy,MYFLT angle,MYFLT * x,MYFLT * y)69 static void rotate_point(MYFLT cx, MYFLT cy, MYFLT angle, MYFLT *x, MYFLT *y)
70 {
71 if(angle == 0) return;
72 MYFLT s = SIN(angle);
73 MYFLT c = COS(angle);
74
75 *x -= cx;
76 *y -= cy;
77
78 float xnew = *x * c - *y * s;
79 float ynew = *x * s + *y * c;
80
81 *x = xnew + cx;
82 *y = ynew + cy;
83 }
84
85 /* the normal eclipse function with center kx,ky and radius krx and kry */
86
ellipse(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)87 static void ellipse(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
88 double x = t+kparam*SIN(t);
89 *outX = kx + krx * SIN(x);
90 *outY = ky + kry * COS(x);
91 }
92
93 /* the limacon curve parametrized by the kparam value
94 see e.g. http://www.2dcurves.com/roulette/roulettel.html
95 for kparam = 1 we have a cardioid
96 */
97
limacon(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)98 static void limacon(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
99 *outX = kx + krx * SIN(t) * (COS(t) + kparam);
100 *outY = ky + kry * COS(t) * (COS(t) + kparam);
101 }
102
103 /* a simple 8 */
104
lemniskateG(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)105 static void lemniskateG(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
106 double x = t+kparam*SIN(t);
107 *outX = kx + krx * COS(x);
108 *outY = ky + kry * SIN(x)*COS(x);
109 }
110
111 /* the cornoid curve
112 see e.g. http://www.2dcurves.com/sextic/sexticco.html
113 */
cornoid(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)114 static void cornoid(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
115 *outX = kx + krx * COS(t) * COS(2*t);
116 *outY = ky + kry * SIN(t) * (kparam + COS(2*t));
117 }
118
119 /* Chevas trisextix
120 see e.g. http://www.2dcurves.com/sextic/sextict.html
121 */
trisec(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)122 static void trisec(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
123 *outX = kx + krx * COS(t) * (1+kparam*SIN(2*t));
124 *outY = ky + kry * SIN(t) * (1+kparam*SIN(2*t));
125 }
126
127 /* Scarabeus curve see e.g http://www.2dcurves.com/sextic/sexticsc.html
128 */
129
scarabeus(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)130 static void scarabeus(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
131 *outX = kx + krx * COS(t) * (kparam*SIN(2*t)+SIN(t));
132 *outY = ky + kry * SIN(t) * (kparam*SIN(2*t)+SIN(t));
133 }
134 /* folium see http://www.2dcurves.com/quartic/quarticfo.html */
folium(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)135 static void folium(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
136 double sint = SIN(t);
137 double cost = COS(t);
138 *outX = kx + krx * cost * cost * (sint*sint - kparam);
139 *outY = ky + kry * sint * cost * (sint*sint - kparam);
140 }
141
142 /* talbot see http://www.2dcurves.com/trig/trigta.html */
talbot(MYFLT t,MYFLT kx,MYFLT ky,MYFLT krx,MYFLT kry,MYFLT kparam,MYFLT * outX,MYFLT * outY)143 static void talbot(MYFLT t, MYFLT kx, MYFLT ky, MYFLT krx, MYFLT kry, MYFLT kparam, MYFLT *outX, MYFLT *outY ) {
144 double sint = SIN(t);
145 double cost = COS(t);
146 *outX = kx + krx * cost * (1 + kparam * sint*sint);
147 *outY = ky + kry * sint * (1 - kparam - kparam*cost*cost);
148 }
149
150 static void (*ifuncs[8])(MYFLT,MYFLT,MYFLT,MYFLT,MYFLT,MYFLT,MYFLT*,MYFLT*) = { ellipse, lemniskateG, limacon, cornoid, trisec, scarabeus, folium, talbot };
151
wtinit(CSOUND * csound,WAVETER * p)152 static int32_t wtinit(CSOUND *csound, WAVETER *p)
153 {
154 p->xarr = NULL;
155 p->yarr = NULL;
156
157 p->oldfnx = -1;
158 p->oldfny = -1;
159 p->sizx = 0;
160 p->sizy = 0;
161 p->theta = 0.0;
162 return OK;
163 }
164
wtPerf(CSOUND * csound,WAVETER * p)165 static int32_t wtPerf(CSOUND *csound, WAVETER *p)
166 {
167 uint32_t offset = p->h.insdshead->ksmps_offset;
168 uint32_t early = p->h.insdshead->ksmps_no_end;
169 uint32_t i, nsmps = CS_KSMPS;
170 int32_t xloc, yloc;
171 MYFLT xc, yc;
172 MYFLT amp = *p->kamp;
173 MYFLT pch = *p->kpch;
174
175 if (*(p->ktabx) != p->oldfnx || p->xarr == NULL) {
176 p->oldfnx = *(p->ktabx);
177 FUNC *ftp = csound->FTFindP(csound, p->ktabx); /* new table parameters */
178 if (UNLIKELY((ftp == NULL) || ((p->xarr = ftp->ftable) == NULL))) return NOTOK;
179 p->sizx = (MYFLT)ftp->flen;
180 }
181 if (*(p->ktaby) != p->oldfny || p->yarr == NULL) {
182 p->oldfny = *(p->ktaby);
183 FUNC *ftp = csound->FTFindP(csound, p->ktaby); /* new table parameters */
184 if (UNLIKELY((ftp == NULL) || ((p->yarr = ftp->ftable) == NULL))) return NOTOK;
185 p->sizy = (MYFLT)ftp->flen;
186 }
187
188
189 uint32_t kfunc = (uint32_t)*p->kfunc;
190 if(kfunc>7) kfunc = 7;
191 MYFLT period = 1;
192 MYFLT sizx = p->sizx, sizy = p->sizy;
193 MYFLT theta = p->theta;
194 MYFLT *aout = p->aout;
195
196 if (UNLIKELY(offset)) memset(aout, '\0', offset*sizeof(MYFLT));
197 if (UNLIKELY(early)) {
198 nsmps -= early;
199 memset(&aout[nsmps], '\0', early*sizeof(MYFLT));
200 }
201 for (i=offset; i<nsmps; i++) {
202
203 /* COMPUTE LOCATION OF SCANNING POINT */
204 ifuncs[kfunc](theta,*p->kx,*p->ky,*p->krx,*p->kry,*p->kparam,&xc,&yc);
205 rotate_point(*p->kx,*p->ky,*p->krot,&xc,&yc);
206 /* MAP SCANNING POINT TO BE IN UNIT SQUARE */
207 xc = xc-FLOOR(xc);
208 yc = yc-FLOOR(yc);
209
210 /* SCALE TO TABLE-SIZE SQUARE */
211 xloc = (int32_t)(xc * sizx);
212 yloc = (int32_t)(yc * sizy);
213
214 /* OUTPUT AM OF TABLE VALUES * KAMP */
215 aout[i] = p->xarr[xloc] * p->yarr[yloc] * amp;
216
217 /* MOVE SCANNING POINT ROUND THE ELLIPSE */
218 theta += pch*((period*TWOPI_F) / csound->GetSr(csound));
219 }
220
221 p->theta = theta;
222 return OK;
223 }
224
225 #define S(x) sizeof(x)
226
227 static OENTRY localops[] = {
228 { "wterrain2", S(WAVETER), TR, 3, "a", "kkkkkkkkkkk",
229 (SUBR)wtinit, (SUBR)wtPerf },
230 };
231
232 LINKAGE
233