1 /***************************************************************************
2 * fractal.cpp is part of Math Graphic Library
3 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru> *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU Lesser General Public License as *
7 * published by the Free Software Foundation; either version 3 of the *
8 * License, or (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU Lesser General Public *
16 * License along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 #include "mgl2/other.h"
21 #include "mgl2/data.h"
22 MGL_NO_EXPORT char *mgl_read_gz(gzFile fp);
23 //-----------------------------------------------------------------------------
24 //
25 // IFS series
26 //
27 //-----------------------------------------------------------------------------
mgl_ifs_2d_point(HCDT A,mreal & x,mreal & y,mreal amax)28 void static mgl_ifs_2d_point(HCDT A, mreal& x, mreal& y, mreal amax)
29 {
30 long i, n=A->GetNy();
31 mreal r = amax*mgl_rnd(), sum_prob = 0, x1;
32 for(i=0;i<n;i++)
33 {
34 sum_prob += A->v(6,i);
35 if(r<sum_prob) break;
36 }
37 x1= A->v(0,i)*x + A->v(1,i)*y + A->v(4,i);
38 y = A->v(2,i)*x + A->v(3,i)*y + A->v(5,i); x = x1;
39 }
mgl_data_ifs_2d(HCDT A,long n,long skip)40 HMDT MGL_EXPORT mgl_data_ifs_2d(HCDT A, long n, long skip)
41 {
42 if(!A || A->GetNx()<7 || n<1) return 0; // incompatible dimensions
43 mreal amax=0;
44 for(long i=0; i<A->GetNy(); i++) amax += A->v(6,i);
45 if(amax<=0) return 0;
46
47 mglData *f = new mglData(2,n);
48 mreal x = 0, y = 0;
49 for(long i=0; i<skip; i++) mgl_ifs_2d_point(A, x, y,amax);
50 for(long i=0; i<n; i++)
51 {
52 mgl_ifs_2d_point(A, x, y, amax);
53 f->a[2*i] = x; f->a[2*i+1] = y;
54 }
55 return f;
56 }
mgl_data_ifs_2d_(uintptr_t * d,long * n,long * skip)57 uintptr_t MGL_EXPORT mgl_data_ifs_2d_(uintptr_t *d, long *n, long *skip)
58 { return uintptr_t(mgl_data_ifs_2d(_DT_,*n,*skip)); }
59 //-----------------------------------------------------------------------------
mgl_ifs_3d_point(HCDT A,mreal & x,mreal & y,mreal & z,mreal amax)60 void static mgl_ifs_3d_point(HCDT A, mreal& x, mreal& y, mreal& z, mreal amax)
61 {
62 int i, n=A->GetNy();
63 mreal r = amax*mgl_rnd(), sum_prob = 0, x1, y1;
64 for (i=0; i<n; i++)
65 {
66 sum_prob += A->v(12,i);
67 if(r<sum_prob) break;
68 }
69 x1= A->v(0,i)*x + A->v(1,i)*y + A->v(2,i)*z + A->v(9,i);
70 y1= A->v(3,i)*x + A->v(4,i)*y + A->v(5,i)*z + A->v(10,i);
71 z = A->v(6,i)*x + A->v(7,i)*y + A->v(8,i)*z + A->v(11,i);
72 x = x1; y = y1;
73 }
mgl_data_ifs_3d(HCDT A,long n,long skip)74 HMDT MGL_EXPORT mgl_data_ifs_3d(HCDT A, long n, long skip)
75 {
76 if(!A || A->GetNx()<13 || n<1) return 0; // incompatible dimensions
77 mreal amax = 0;
78 for(int i=0; i<A->GetNy(); i++) amax += A->v(12,i);
79 if(amax <= 0) return 0;
80
81 mglData *f = new mglData(3,n);
82 mreal x = 0, y = 0, z = 0;
83 for(long i=0; i<skip; i++) mgl_ifs_3d_point(A, x, y, z, amax);
84 for(long i=0; i<n; i++)
85 {
86 mgl_ifs_3d_point(A, x, y, z, amax);
87 f->a[3*i] = x; f->a[3*i+1] = y; f->a[3*i+2] = z;
88 }
89 return f;
90 }
mgl_data_ifs_3d_(uintptr_t * d,long * n,long * skip)91 uintptr_t MGL_EXPORT mgl_data_ifs_3d_(uintptr_t *d, long *n, long *skip)
92 { return uintptr_t(mgl_data_ifs_3d(_DT_,*n,*skip)); }
93 //-----------------------------------------------------------------------------
mgl_data_ifs_file(const char * fname,const char * name,long n,long skip)94 HMDT MGL_EXPORT mgl_data_ifs_file(const char *fname, const char *name, long n, long skip)
95 {
96 gzFile fp = gzopen(fname,"r");
97 if(!fp) return 0; // Couldn't open file file
98 char *buf = mgl_read_gz(fp); gzclose(fp);
99 char *s = strstr(buf,name);
100 if(!s) return 0; // No data for fractal 'name' in the file
101
102 char *p = strchr(s,'{'), *e;
103 if(!p) return 0; // Wrong data format for fractal 'name' in the file
104 bool ext3d = false;
105 e = strstr(s,"(3D)"); if(e && e<p) ext3d = true;
106 e = strstr(s,"(3d)"); if(e && e<p) ext3d = true;
107 e = strchr(p,'}');
108
109 std::vector<mreal> nums;
110 for(size_t i=0;p[i] && p+i<e;i++)
111 {
112 while(p[i]<=' ') i++;
113 if(p[i]==';' || p[i]=='#') while(p[i] && p[i]!='\n') i++;
114 if(strchr("0123456789.+-",p[i])) // this is number
115 {
116 nums.push_back(atof(p+i));
117 while(p[i]>' ') i++;
118 }
119 }
120 HMDT dat = new mglData, res;
121 if(ext3d)
122 {
123 dat->Set(&(nums[0]), 13, nums.size()/13, 1);
124 res = mgl_data_ifs_3d(dat, n, skip);
125 }
126 else
127 {
128 dat->Set(&(nums[0]), 7, nums.size()/7, 1);
129 res = mgl_data_ifs_2d(dat, n, skip);
130 }
131 delete dat; free(buf); return res;
132 }
mgl_data_ifs_file_(const char * fname,const char * name,long * n,long * skip,int l,int m)133 uintptr_t mgl_data_ifs_file_(const char *fname, const char *name, long *n, long *skip,int l,int m)
134 { char *s=new char[l+1]; memcpy(s,fname,l); s[l]=0;
135 char *t=new char[m+1]; memcpy(t,name,m); t[m]=0;
136 uintptr_t r = uintptr_t(mgl_data_ifs_file(s,t,*n,*skip));
137 delete []s; delete []t; return r; }
138 //-----------------------------------------------------------------------------
139 //
140 // Functions for flame fractal
141 //
142 //-----------------------------------------------------------------------------
mgl_linear_var0(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)143 void static mgl_linear_var0(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
144 { xNew += par[0]*x; yNew += par[0]*y; }
mgl_sinusoidal_var1(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)145 void static mgl_sinusoidal_var1(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
146 { xNew += par[0]*sin(x); yNew += par[0]*sin(y); }
mgl_spherical_var2(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)147 void static mgl_spherical_var2(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
148 { mreal c1 = par[0]/(x*x+y*y); xNew += c1*x; yNew += c1*y; }
mgl_swirl_var3(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)149 void static mgl_swirl_var3(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
150 {
151 mreal r2=x*x+y*y, c1=sin(r2), c2=cos(r2);
152 xNew += par[0]*(x*c1 - y*c2);
153 yNew += par[0]*(x*c2 + y*c1);
154 }
mgl_horseshoe_var4(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)155 void static mgl_horseshoe_var4(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
156 {
157 mreal c1 = par[0]/hypot(x,y);
158 xNew += c1*(x*x-y*y);
159 yNew += 2*c1*x*y;
160 }
mgl_polar_var5(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)161 void static mgl_polar_var5(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
162 {
163 xNew += par[0]*atan2(x,y)/M_PI;
164 yNew += par[0]*(hypot(x,y)-1);
165 }
mgl_handkerchief_var6(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)166 void static mgl_handkerchief_var6(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
167 {
168 mreal r=hypot(x,y), t=atan2(x,y), c1=par[0]*r;
169 xNew += c1*sin(t+r); yNew += c1*cos(t-r);
170 }
mgl_heart_var7(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)171 void static mgl_heart_var7(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
172 {
173 mreal r=hypot(x,y), c1=par[0]*r, c2=atan2(x,y)*r;
174 xNew += c1*sin(c2); yNew -= c1*cos(c2);
175 }
mgl_disc_var8(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)176 void static mgl_disc_var8(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
177 {
178 mreal c1=par[0]*atan2(x,y)/M_PI, c2=M_PI*hypot(x,y);
179 xNew += c1*sin(c2); yNew += c1*cos(c2);
180 }
mgl_spiral_var9(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)181 void static mgl_spiral_var9(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
182 {
183 mreal r=hypot(x,y), t=atan2(x,y), c1=par[0]/r;
184 xNew += c1*(cos(t)+sin(r));
185 yNew += c1*(sin(t)-cos(r));
186 }
mgl_hyperbolic_var10(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)187 void static mgl_hyperbolic_var10(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
188 {
189 mreal r=hypot(x,y), t=atan2(x,y);
190 xNew += par[0]*sin(t)/r;
191 yNew += par[0]*r*cos(t);
192 }
mgl_diamond_var11(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)193 void static mgl_diamond_var11(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
194 {
195 mreal r=hypot(x,y), t=atan2(x,y);
196 xNew += par[0]*sin(t)*cos(r);
197 yNew += par[0]*cos(t)*sin(r);
198 }
mgl_ex_var12(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)199 void static mgl_ex_var12(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
200 {
201 mreal r=hypot(x,y), t=atan2(x,y), c1=par[0]*r;
202 mreal c2=mgl_ipow(sin(t+r),3), c3 = mgl_ipow(cos(t-r), 3);
203 xNew += c1*(c2 + c3); yNew += c1*(c2 - c3);
204 }
mgl_julia_var13(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)205 void static mgl_julia_var13(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
206 {
207 mreal c1=par[0]*sqrt(hypot(x,y)), c2=atan2(x,y)/2, c3=(rand()%2)*M_PI;
208 xNew += c1*cos(c2+c3); yNew += c1*sin(c2+c3);
209 }
mgl_bent_var14(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)210 void static mgl_bent_var14(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
211 {
212 if(x>=0 && y>=0)
213 { xNew += par[0]*x; yNew += par[0]*y; }
214 else if(x<0 && y>=0)
215 { xNew += par[0]*2*x; yNew += par[0]*y; }
216 else if(x>=0 && y<0)
217 { xNew += par[0]*x; yNew += par[0]*y/2; }
218 else
219 { xNew += par[0]*2*x; yNew += par[0]*y/2; }
220 }
mgl_waves_var15(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)221 void static mgl_waves_var15(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
222 { // NOTE: par[1]=b[i], par[2]=1/c[i]^2, par[3]=e[i], par[4]=1/f[i]^2
223 xNew += par[0]*(x + par[1]*sin(y*par[2]));
224 yNew += par[0]*(y + par[3]*sin(x*par[4]));
225 }
mgl_fisheye_var16(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)226 void static mgl_fisheye_var16(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
227 {
228 mreal c1 = par[0]*2/(hypot(x,y) + 1);
229 xNew += c1*y; yNew += c1*x;
230 }
mgl_popcorn_var17(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)231 void static mgl_popcorn_var17(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
232 { // NOTE: par[1]=c[i], par[2]=f[i]
233 xNew += par[0]*(x + par[1]*sin(tan(3*y)));
234 yNew += par[0]*(y + par[2]*sin(tan(3*x)));
235 }
mgl_exponential_var18(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)236 void static mgl_exponential_var18(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
237 {
238 mreal c1=par[0]*exp(x-1);
239 xNew += c1*cos(M_PI*y); yNew += c1*sin(M_PI*y);
240 }
mgl_power_var19(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)241 void static mgl_power_var19(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
242 {
243 mreal t=atan2(x,y), c1=par[0]*pow(hypot(x,y), sin(t));
244 xNew += c1*cos(t); yNew += c1*sin(t);
245 }
mgl_cosine_var20(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)246 void static mgl_cosine_var20(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
247 {
248 xNew += par[0]*cos(M_PI*x)*cosh(y);
249 yNew -= par[0]*sin(M_PI*x)*sinh(y);
250 }
mgl_rings_var21(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)251 void static mgl_rings_var21(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
252 { // NOTE: par[1]=c[i]^2
253 mreal t=atan2(x,y), r=hypot(x,y), c1=par[0]*(fmod(r+par[1],2*par[1])-par[1]+r*(1-par[1])); // convert to int?
254 xNew += c1*cos(t); yNew += c1*sin(t);
255 }
mgl_fan_var22(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)256 void static mgl_fan_var22(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
257 { // NOTE: par[1]=c[i]^2, par[2]=f[i]
258 mreal t=atan2(x,y), c1=par[0]*hypot(x,y), c2;
259 c2 = fmod(t+par[2], M_PI*par[1]); // convert to int?
260 if(c2>M_PI_2*par[1]) c2 = t - M_PI_2*par[1];
261 else c2 += t;
262 xNew += c1*cos(c2); yNew += c1*sin(c2);
263 }
mgl_blob_var23(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)264 void static mgl_blob_var23(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
265 {
266 mreal t=atan2(x,y), c1=par[0]*hypot(x,y)*(par[2]+(par[1]-par[2])/2*(sin(par[3]*t)));
267 xNew += c1*cos(t); yNew += c1*sin(t);
268 }
mgl_pdj_var24(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)269 void static mgl_pdj_var24(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
270 {
271 xNew += par[0]*(sin(par[1]*y) - cos(par[2]*x));
272 yNew += par[0]*(sin(par[3]*x) - cos(par[4]*y));
273 }
mgl_fan2_var25(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)274 void static mgl_fan2_var25(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
275 {
276 mreal t=atan2(x,y), c1, c2;
277 c1 = M_PI*par[1]*par[1];
278 c2 = t + par[2] - c1*int(2*t*par[2]/c1);
279 c1 /= 2; c2 = c2>c1?t-c1:t+c1;
280 c1 = par[0]*hypot(x,y);
281 xNew += c1*sin(c2); yNew += c1*cos(c2);
282 }
mgl_rings2_var26(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)283 void static mgl_rings2_var26(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
284 {
285 mreal r=hypot(x,y), t=atan2(x,y), c1=par[1]*par[1];
286 c1 = par[0]*(r - 2*c1*int((r+c1)/(2*c1)) + r*(1-c1));
287 xNew += c1*cos(t); yNew += c1*sin(t);
288 }
mgl_eyefish_var27(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)289 void static mgl_eyefish_var27(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
290 {
291 mreal c1 = par[0]*2/(hypot(x,y)+1);
292 xNew += c1*x; yNew += c1*y;
293 }
mgl_bubble_var28(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)294 void static mgl_bubble_var28(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
295 {
296 mreal c1 = par[0]*4/(x*x+y*y+4);
297 xNew += c1*x; yNew += c1*y;
298 }
mgl_cylinder_var29(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)299 void static mgl_cylinder_var29(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
300 { xNew += par[0]*sin(x); yNew += par[0]*y; }
mgl_perspective_var30(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)301 void static mgl_perspective_var30(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
302 {
303 mreal c1 = par[0]*par[2]/(par[2]-y*sin(par[1]));
304 xNew += c1*x; yNew += c1*y*cos(par[1]);
305 }
mgl_noise_var31(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)306 void static mgl_noise_var31(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
307 {
308 mreal c1=par[0]*mgl_rnd(), c2=2*M_PI*mgl_rnd();
309 xNew += c1*x*cos(c2); yNew += c1*y*sin(c2);
310 }
mgl_juliaN_var32(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)311 void static mgl_juliaN_var32(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
312 {
313 mreal c1=int(fabs(par[1])*mgl_rnd()), c2;
314 c2 = (atan2(y,x) + 2*M_PI*c1)/par[1];
315 c1 = par[0]*pow(hypot(x,y), par[2]/par[1]);
316 xNew += c1*cos(c2); yNew += c1*sin(c2);
317 }
mgl_juliaScope_var33(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)318 void static mgl_juliaScope_var33(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
319 {
320 mreal c1=int(fabs(par[1])*mgl_rnd()), c2;
321 c2 = ((2*(rand()%2)-1)*atan2(y,x) + 2*M_PI*c1)/par[1];
322 c1 = par[0]*pow(hypot(x,y), par[2]/par[1]);
323 xNew += c1*cos(c2); yNew += c1*sin(c2);
324 }
mgl_blur_var34(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)325 void static mgl_blur_var34(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
326 {
327 mreal c1=par[0]*mgl_rnd(), c2=2*M_PI*mgl_rnd();
328 xNew += c1*cos(c2); yNew += c1*sin(c2);
329 }
mgl_gaussian_var35(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)330 void static mgl_gaussian_var35(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
331 {
332 mreal c1=par[0]*(4*mgl_rnd()-2), c2=2*M_PI*mgl_rnd();
333 xNew += c1*cos(c2); yNew += c1*sin(c2);
334 }
mgl_radialBlur_var36(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)335 void static mgl_radialBlur_var36(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
336 {
337 mreal r=hypot(x,y), c1=par[1]*M_PI_2, c2=par[0]*(4*mgl_rnd()-2), c3;
338 c3 = c2*cos(c1) - 1; c2 = atan2(y,x) + c2 *sin(c1);
339 xNew += r*cos(c2) + c3*x; yNew += r*sin(c2) + c3*y;
340 }
mgl_pie_var37(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)341 void static mgl_pie_var37(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
342 {
343 mreal c1=int(mgl_rnd()*par[1] + 0.5), c2;
344 c2 = par[2] + 2*M_PI/par[1]*(c1 + mgl_rnd()*par[3]);
345 c1 = par[0]*mgl_rnd();
346 xNew += c1*cos(c2); yNew += c1*sin(c2);
347 }
mgl_ngon_var38(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)348 void static mgl_ngon_var38(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
349 {
350 mreal c1=2*M_PI/par[2], c2;
351 c2 = atan2(y,x) - c1*floor(atan2(y,x)/c1);
352 if(c2 <= c1/2) c2 -= c1;
353 c1 = par[0]*(par[3]*(1/cos(c2) - 1) + par[4])/pow(hypot(x,y), par[1]);
354 xNew += c1*x; yNew += c1*y;
355 }
mgl_curl_var39(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)356 void static mgl_curl_var39(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
357 {
358 mreal c1=1 + par[1]*x + par[2]*(x*x - y*y);
359 mreal c2 = par[1]*y + 2*par[2]*x*y;
360 mreal c3 = par[0]/(c1*c1 + c2*c2);
361 xNew += c3*(c1*x + c2*y); yNew += c3*(c1*x - c2*y);
362 }
mgl_rectangles_var40(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)363 void static mgl_rectangles_var40(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
364 {
365 xNew += par[0]*((2*floor(x/par[1]) + 1)*par[1] - x);
366 yNew += par[0]*((2*floor(y/par[2]) + 1)*par[2] - y);
367 }
mgl_arch_var41(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)368 void static mgl_arch_var41(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
369 {
370 mreal c1=mgl_rnd()*M_PI*par[0], c2=sin(c1);
371 xNew += par[0]*c2; yNew += par[0]*c2*c2/cos(c1);
372 }
mgl_tangent_var42(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)373 void static mgl_tangent_var42(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
374 { xNew += par[0]*sin(x)/cos(y); yNew += par[0]*tan(y); }
mgl_square_var43(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)375 void static mgl_square_var43(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
376 {
377 xNew += par[0]*(mgl_rnd() - 0.5);
378 yNew += par[0]*(mgl_rnd() - 0.5);
379 }
mgl_blade_var44(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)380 void static mgl_blade_var44(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
381 {
382 mreal c1=par[0]*x, c2=mgl_rnd()*hypot(x,y)*par[0];
383 xNew += c1*(cos(c2) + sin(c2)); // TODO check use of c2
384 yNew += c1*(cos(c2) - sin(c2));
385 }
mgl_secant_var45(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)386 void static mgl_secant_var45(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
387 { xNew += par[0]*x; yNew += 1/cos(par[0]*hypot(x,y)); }
mgl_rays_var46(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)388 void static mgl_rays_var46(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
389 {
390 mreal c1 = par[0]*par[0]*tan(mgl_rnd()*M_PI*par[0])/(x*x+y*y);
391 xNew += c1*cos(x); yNew += c1*sin(y);
392 }
mgl_twintrian_var47(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)393 void static mgl_twintrian_var47(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
394 {
395 mreal c1=par[0]*x, c2, c3;
396 c2 = mgl_rnd()*hypot(x,y)*par[0];
397 c3 = log10(sin(c2)*sin(c2)) + cos(c2);
398 xNew += c1*c3; yNew += c1*(c3 - M_PI*sin(c2));
399 }
mgl_cross_var48(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)400 void static mgl_cross_var48(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
401 {
402 mreal c1 = par[0]/fabs(x*x - y*y);
403 xNew += c1*x; yNew += c1*y;
404 }
mgl_disc2_var49(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)405 void static mgl_disc2_var49(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
406 {
407 mreal c1 = cos(par[2])-1, c2 = sin(par[2]), c3;
408 if(par[2]>2*M_PI)
409 {
410 c3 = 1+par[2]-2*M_PI;
411 c1 *= c3; c2 *= c3;
412 }
413 else if(par[2]<-2*M_PI)
414 {
415 c3 = 1+par[2]+2*M_PI;
416 c1 *= c3; c2 *= c3;
417 }
418 c3 = par[1]*M_PI*(x+y);
419 mreal a = par[0]*atan2(x,y)/M_PI;
420 xNew += a*(c1+sin(c3)); yNew += a*(c2+cos(c3));
421 }
mgl_supershape_var50(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)422 void static mgl_supershape_var50(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
423 {
424 mreal c1 = par[2]/4*atan2(y,x)+M_PI_4, r = hypot(x,y);
425 c1 = pow(fabs(sin(c1)), par[5])+pow(fabs(cos(c1)), par[4]);
426 c1 = par[0]*((par[1]*mgl_rnd()+(1-par[1])*r)-par[6])*pow(c1, -1.0/par[3])/r;
427 xNew += c1*x; yNew += c1*y;
428 }
mgl_flower_var51(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)429 void static mgl_flower_var51(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
430 {
431 mreal c1 = par[0]*(mgl_rnd()-par[2])*cos(par[1]*atan2(y,x))/hypot(x,y);
432 xNew += c1*x; yNew += c1*y;
433 }
mgl_conic_var52(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)434 void static mgl_conic_var52(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
435 {
436 mreal c1 = x/hypot(x,y);
437 c1 = par[0]*(mgl_rnd()-par[2])*par[1]/(1+par[1]*c1)/hypot(x,y);
438 xNew += c1*x; yNew += c1*y;
439 }
mgl_parabola_var53(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)440 void static mgl_parabola_var53(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
441 {
442 mreal c1 = hypot(x,y), c2;
443 c2 = cos(c1); c1 = sin(c1);
444 xNew += par[0]*par[1]*c1*c1 *mgl_rnd();
445 yNew += par[0]*par[2]*c2*mgl_rnd();
446 }
mgl_bent2_var54(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)447 void static mgl_bent2_var54(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
448 {
449 mreal c1 = par[0]*x, c2 = par[0]*y;
450 if(x<0.0) c1 *= par[1];
451 if(y<0.0) c2 *= par[2];
452 xNew += c1; yNew += c2;
453 }
mgl_bipolar_var55(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)454 void static mgl_bipolar_var55(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
455 {
456 mreal c1 = 0.5*atan2(2*y, x*x+y*y-1)-M_PI_2*par[1], c2, c3, c4;
457 if(c1>M_PI_2) c1 = -M_PI_2+fmod(c1+M_PI_2, M_PI);
458 else if(y<-M_PI_2) c1 = M_PI_2-fmod(M_PI_2-y, M_PI);
459 c2 = par[0]*M_2_PI;
460 c3 = x*x+y*y+1; c4 = 2*x;
461 xNew += 0.25*c2*log((c3+c4)/(c3-c4)); yNew += c1*c2;
462 }
mgl_boarders_var56(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)463 void static mgl_boarders_var56(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
464 {
465 mreal c1 = round(x), c2 = round(y), c3, c4;
466 c3 = x-c1; c4 = y-c2;
467 if(mgl_rnd()>=0.75)
468 {
469 xNew += par[0]*(c1+0.5*c3);
470 yNew += par[0]*(c2+0.5*c4);
471 }
472 else
473 {
474 if(fabs(c3)>=fabs(c4))
475 {
476 if(c3>=0)
477 {
478 xNew += par[0]*(c1+0.5*c3+0.25);
479 yNew += par[0]*(c2+0.5*c4+0.25*c4/c3);
480 }
481 else
482 {
483 xNew += par[0]*(c1+0.5*c3-0.25);
484 yNew += par[0]*(c2+0.5*c4-0.25*c4/c3);
485 }
486 }
487 else
488 {
489 if(c4>=0)
490 {
491 xNew += par[0]*(c2+0.5*c4+0.25);
492 yNew += par[0]*(c1+0.5*c3+0.25*c3/c4);
493 }
494 else
495 {
496 xNew += par[0]*(c2+0.5*c4-0.25);
497 yNew += par[0]*(c1+0.5*c3-0.25*c3/c4);
498 }
499 }
500 }
501 }
mgl_butterfly_var57(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)502 void static mgl_butterfly_var57(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
503 {
504 mreal c1 = 4/sqrt(3*M_PI)*par[0], c2 = 2*y; // replace 4/sqrt(3*M_PI) for the result?
505 c1 *= sqrt(fabs(x*y)/(x*x+c2*c2));
506 xNew += c1*x; yNew += c1*c2;
507 }
mgl_cell_var58(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)508 void static mgl_cell_var58(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
509 {
510 mreal c1=floor(x/par[1]), c2=floor(y/par[1]), c3=x-c1*par[1], c4=y-c2*par[1];
511 if(c2>=0)
512 {
513 if(c1>=0) { c1 *= 2; c2 *= 2; }
514 else { c1=-2*c1-1; c2 *= 2; }
515 }
516 else
517 {
518 if(c1>=0) { c1 *= 2; c2=-2*c2-1; }
519 else { c1=-2*c1-1; c2=-2*c2-1; }
520 }
521 xNew += par[0]*(par[1]*c1+c3);
522 yNew += par[0]*(par[1]*c2+c4);
523 }
mgl_cpow_var59(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)524 void static mgl_cpow_var59(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
525 {
526 mreal c1 = atan2(y, x), c2, c3, c4, c5;
527 c2 = par[1]/par[3];
528 c3 = par[2]/par[3];
529 c4 = 0.5*log(x*x+y*y);
530 c5 = c1*c2+c3*c4+2*M_PI/par[3]*floor(par[3]*mgl_rnd());
531 c1 = par[0]*exp(c2*c4-c1*c3);
532 c2 = cos(c5);
533 c3 = sin(c5);
534 xNew += c1*c2; yNew += c1*c3;
535 }
mgl_curve_var60(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)536 void static mgl_curve_var60(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
537 {
538 xNew += par[0]*(x+par[1]*exp(-(y*y)/(par[3]*par[3])));
539 yNew += par[0]*(y+par[2]*exp(-(x*x)/(par[4]*par[4])));
540 }
mgl_edisc_var61(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)541 void static mgl_edisc_var61(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
542 {
543 mreal c1 = x*x+y*y+1, c2,c3,c4,c5;
544 c2 = 2*x;
545 c3 = sqrt(c1-c2);
546 c1 = sqrt(c1+c2);
547 c1 = (c1+c3)*0.5;
548 c2 = log(c1+sqrt(c1-1));
549 c3 = -acos(x/c1);
550 c1 = par[0]/11.57034632;
551 c4 = cos(c2)*cosh(c3);
552 c5 = sin(c2)*sinh(c3);
553 xNew += c1*c4; yNew += c1*c5;
554 }
mgl_elliptic_var62(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)555 void static mgl_elliptic_var62(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
556 {
557 mreal c1 = x*x+y*y+1, c2, c3, c4, c5;
558 c2 = 2*x;
559 c2 = 0.5*(sqrt(c1+c2)+sqrt(c1-c2));
560 c3 = x/c2;
561 c4 = 1-c3*c3;
562 c5 = c2-1;
563 c1 = par[0]/M_PI_2;
564 if(c4<0) c4 = 0;
565 else c4 = sqrt(c4);
566 if(c5<0) c5 = 0;
567 else c5 = sqrt(c5);
568 xNew += c1*atan2(c3, c4);
569 if(y>0) yNew += c1*log(c2+c5);
570 else yNew -= c1*log(c2+c5);
571 }
mgl_escher_var63(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)572 void static mgl_escher_var63(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
573 {
574 mreal c1 = 0.5*(1+cos(par[1])), c2 = 0.5*sin(par[1]), c3, c4, c5;
575 c3 = 0.5*log(x*x+y*y); c4 = atan2(y, x);
576 c5 = c1*c4+c2*c3; c1 = par[0]*exp(c1*c3-c2*c4);
577 xNew += c1*cos(c5); yNew += c1*sin(c5);
578 }
mgl_foci_var64(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)579 void static mgl_foci_var64(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
580 {
581 mreal c1=0.5*exp(x), c2=0.25/c1, c3=par[0]/(c1+c2-cos(y)); // TODO Check this!!!
582 xNew += c1*(c2-c3); yNew += c1*sin(y);
583 }
mgl_lazySusan_var65(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)584 void static mgl_lazySusan_var65(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
585 {
586 mreal c1 = hypot(x, y), c2;
587 if(c1<par[0])
588 {
589 c2 = atan2(y+par[5], x-par[4])+par[1]+par[3]*(par[0]-c1);
590 c1 *= par[0];
591 xNew += c1*cos(c2)+par[4];
592 yNew += c1*sin(c2)-par[5];
593 }
594 else
595 {
596 c1 = par[0]*(1+par[2]/c1);
597 xNew += c1*(x-par[4])+par[4];
598 yNew += c1*(y+par[5])-par[5];
599 }
600 }
mgl_loonie_var66(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)601 void static mgl_loonie_var66(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
602 {
603 mreal c1 = x*x+y*y, c2 = par[0]*par[0];
604 if(c1<c2)
605 {
606 c1 = par[0]*sqrt(c2/c1-1);
607 xNew += c1*x; yNew += c1*y;
608 }
609 else { xNew += par[0]*x; yNew += par[0]*y; }
610 }
mgl_preBlur_var67(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)611 void static mgl_preBlur_var67(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
612 {
613 mreal c1 = par[0]*(mgl_rnd()+mgl_rnd()+mgl_rnd()+mgl_rnd()-2), c2;
614 c2 = 2*mgl_rnd()*M_PI;
615 x += c1*cos(c2); y += c1*sin(c2); // NOTE: This changes the original coordinates, not the new ones
616 }
mgl_modulus_var68(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)617 void static mgl_modulus_var68(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
618 {
619 if(x>par[1]) xNew += par[0]*(-par[1]+fmod(x-par[1], 2*par[1]));
620 else if(x<par[1]) xNew += par[0]*(par[1]-fmod(par[1]-x, 2*par[1]));
621 else xNew += par[0]*x;
622 if(y>par[2]) yNew += par[0]*(-par[2]+fmod(y+par[2], 2*par[2]));
623 else if(y<par[2]) yNew += par[0]*(par[2]-fmod(par[2]-y, 2*par[2]));
624 else yNew += par[0]*y;
625 }
mgl_oscope_var69(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)626 void static mgl_oscope_var69(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
627 {
628 mreal c1 = par[3]*exp(-fabs(x)*par[4])*cos(2*M_PI*par[2]*x)+par[1];
629 if(fabs(y) <= c1)
630 { xNew += par[0]*x; yNew -= par[0]*y; }
631 else
632 { xNew += par[0]*x; yNew += par[0]*y; }
633 }
mgl_polar2_var70(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)634 void static mgl_polar2_var70(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
635 {
636 mreal t = atan2(x, y);
637 xNew += par[0]*t*t; yNew += par[0]*t/2*log(x*x+y*y);
638 }
mgl_popcorn2_var71(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)639 void static mgl_popcorn2_var71(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
640 {
641 xNew += par[0]*(x+par[1]*sin(tan(par[3]*y)));
642 yNew += par[0]*(y+par[2]*sin(tan(par[3]*x)));
643 }
mgl_scry_var72(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)644 void static mgl_scry_var72(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
645 {
646 mreal c1 = 1/(hypot(x, y)*(x*x+y*y+1/par[0]));
647 xNew += c1*x; yNew += c1*y;
648 }
mgl_separation_var73(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)649 void static mgl_separation_var73(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
650 {
651 if(x>0) xNew += par[0]*(sqrt(x*x+par[1]*par[1])-x*par[2]);
652 else xNew -= par[0]*(sqrt(x*x+par[1]*par[1])+x*par[2]);
653 if(y>0) yNew += par[0]*(sqrt(y*y+par[3]*par[3])-y*par[4]);
654 else yNew -= par[0]*(sqrt(y*y+par[3]*par[3])+y*par[4]);
655 }
mgl_split_var74(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)656 void static mgl_split_var74(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
657 {
658 if(cos(M_PI*x*par[1])>=0) xNew += par[0]*y;
659 else xNew -= par[0]*y;
660 if(cos(M_PI*y*par[2])>=0) yNew += par[0]*x;
661 else yNew -= par[0]*x;
662 }
mgl_splits_var75(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)663 void static mgl_splits_var75(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
664 {
665 if(x>=0) xNew+= par[0]*(x+par[1]);
666 else xNew += par[0]*(x-par[1]);
667 if(y>=0) yNew += par[0]*(y+par[2]);
668 else yNew += par[0]*(y-par[2]);
669 }
mgl_stripes_var76(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)670 void static mgl_stripes_var76(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
671 {
672 mreal c1 = floor(x+0.5), c2 = x-c1;
673 xNew += par[0]*(c2*(1-par[1])+c1);
674 yNew += par[0]*(y+c2*c2*par[2]);
675 }
mgl_wedge_var77(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)676 void static mgl_wedge_var77(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
677 {
678 mreal c1 = hypot(x, y), c2, c3, c4;
679 c2 = atan2(y, x)+par[4]*c1;
680 c3 = 1-0.5*M_1_PI*par[1]*par[3];
681 c4 = floor(0.5*M_1_PI*(par[2]*c2+M_PI));
682 c2 = c2*c3+c4*par[1];
683 c1 = par[0]*(c1+par[2]);
684 xNew += c1*cos(c2); yNew += c1*sin(c2);
685 }
mgl_wedgeJulia_var78(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)686 void static mgl_wedgeJulia_var78(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
687 {
688 mreal c1 = int(fabs(par[3])*mgl_rnd()), c2;
689 c1 = (atan2(y, x)+2*M_PI*c1)/par[3];
690 c2 = floor(0.5*M_1_PI*(par[2]*c1+M_PI));
691 c2 = c1*(1-0.5*M_1_PI*par[1]*par[2]+c1*par[1]); // TODO Check this!!!
692 c1 = par[0]*pow(x*x +y*y, par[4]/(2*par[3]));
693 xNew += c1*cos(c2); yNew += c1*sin(c2);
694 }
mgl_wedgeSph_var79(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)695 void static mgl_wedgeSph_var79(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
696 {
697 mreal c1 = 1/hypot(x, y), c2, c3, c4;
698 c2 = atan2(y, x)+par[4]*c1;
699 c3 = 1-0.5*M_1_PI*par[1]*par[2];
700 c4 = floor(0.5*M_1_PI*(par[2]*c2+M_PI));
701 c2 = c2*c3+c4*par[1];
702 c1 = par[0]*(c1+par[3]);
703 xNew += c1*cos(c2); yNew += c1*sin(c2);
704 }
mgl_whorl_var80(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)705 void static mgl_whorl_var80(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
706 {
707 mreal c1 = hypot(x, y), c2;
708 if(c1<par[0]) c2 = atan2(y, x)+par[1]/(par[0]-c1);
709 else c2 = atan2(y, x)+par[2]/(par[0]-c1);
710 c1 *= par[0];
711 xNew += c1*cos(c2); yNew += c1*sin(c2);
712 }
mgl_waves2_var81(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)713 void static mgl_waves2_var81(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
714 {
715 xNew += par[0]*(x+par[2]*sin(y*par[1]));
716 yNew += par[0]*(y+par[4]*sin(x*par[3]));
717 }
mgl_exp_var82(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)718 void static mgl_exp_var82(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
719 {
720 mreal c1 = par[0]*exp(x);
721 xNew += c1*cos(y); yNew += c1*sin(y);
722 }
mgl_log_var83(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)723 void static mgl_log_var83(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
724 {
725 xNew += par[0]*0.5*log(x*x+y*y); yNew += par[0]*atan2(y, x);
726 }
mgl_sin_var84(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)727 void static mgl_sin_var84(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
728 {
729 xNew += par[0]*sin(x)*cosh(y); yNew += par[0]*cos(x)*sinh(y);
730 }
mgl_cos_var85(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)731 void static mgl_cos_var85(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
732 {
733 xNew += par[0]*cos(x)*cosh(y); yNew -= par[0]*sin(x)*sinh(y);
734 }
mgl_tan_var86(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)735 void static mgl_tan_var86(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
736 {
737 mreal c1 = par[0]/(cos(2*x)+cosh(2*y));
738 xNew += c1*sin(2*x); yNew += c1*sinh(2*y);
739 }
mgl_sec_var87(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)740 void static mgl_sec_var87(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
741 {
742 mreal c1 = 2*par[0]/(cos(2*x)+cosh(2*y));
743 xNew += c1*cos(x)*cosh(y); yNew += c1*sin(x)*sinh(y);
744 }
mgl_csc_var88(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)745 void static mgl_csc_var88(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
746 {
747 mreal c1 = 2*par[0]/(cosh(2*y)-cos(2*x));
748 xNew += c1*sin(x)*cosh(y); yNew -= c1*cos(x)*sinh(y);
749 }
mgl_cot_var89(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)750 void static mgl_cot_var89(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
751 {
752 mreal c1 = par[0]/(cosh(2*y)-cos(2*x));
753 xNew += c1*sin(2*x); yNew -= c1*sinh(2*y);
754 }
mgl_sinh_var90(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)755 void static mgl_sinh_var90(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
756 { xNew += par[0]*sinh(x)*cos(y); yNew += par[0]*cosh(y)*sin(y); }
mgl_cosh_var91(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)757 void static mgl_cosh_var91(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
758 { xNew += par[0]*cosh(x)*cos(y); yNew += par[0]*sinh(x)*sin(y); }
mgl_tanh_var92(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)759 void static mgl_tanh_var92(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
760 {
761 mreal c1 = par[0]/(cos(2*y)+cosh(2*x));
762 xNew += c1*sinh(2*x); yNew += c1*sin(2*y);
763 }
mgl_sech_var93(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)764 void static mgl_sech_var93(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
765 {
766 mreal c1 = 2*par[0]/(cos(2*y)+cosh(2*x));
767 xNew += c1*cos(y)*cosh(x); yNew -= c1*sin(y)*sinh(x);
768 }
mgl_csch_var94(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)769 void static mgl_csch_var94(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
770 {
771 mreal c1 = 2*par[0]/(cosh(2*x)-cos(2*y));
772 xNew += c1*sinh(x)*cos(y); yNew -= c1*cosh(x)*sin(y);
773 }
mgl_coth_var95(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)774 void static mgl_coth_var95(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
775 {
776 mreal c1 = par[0]/(cosh(2*x)-cos(2*y));
777 xNew += c1*sinh(2*x); yNew += c1*sin(2*y);
778 }
mgl_auger_var96(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)779 void static mgl_auger_var96(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
780 {
781 mreal c1 = sin(par[3]*y), c2 = sin(par[3]*x);
782 c1 = x+par[2]*(par[4]*c1/2+fabs(x)*c1);
783 c2 = y+par[2]*(par[4]*c2/2+fabs(y)*c2);
784 xNew += par[0]*(x+par[1]*(c1-x)); yNew += par[0]*c2;
785 }
mgl_flux_var97(mreal & xNew,mreal & yNew,mreal x,mreal y,const mreal * par)786 void static mgl_flux_var97(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par)
787 {
788 mreal c1 = x+par[0], c2 = x-par[0], c3;
789 c3 = 0.5*(atan2(y, c2)-atan2(y, c1));
790 c1 = par[0]*(2+par[1])*sqrt(sqrt(y*y+c1*c1)/sqrt(y*y-c2*c2));
791 xNew += c1*cos(c3); yNew += c1*sin(c3);
792 }
793 //-----------------------------------------------------------------------------
794 typedef void (*flame_func)(mreal &xNew, mreal &yNew, mreal x, mreal y, const mreal *par);
795 static flame_func ffunc[mglFlame2dLAST] = {
796 mgl_linear_var0, mgl_sinusoidal_var1, mgl_spherical_var2, mgl_swirl_var3, mgl_horseshoe_var4,
797 mgl_polar_var5, mgl_handkerchief_var6, mgl_heart_var7, mgl_disc_var8, mgl_spiral_var9,
798 mgl_hyperbolic_var10, mgl_diamond_var11, mgl_ex_var12, mgl_julia_var13, mgl_bent_var14,
799 mgl_waves_var15, mgl_fisheye_var16, mgl_popcorn_var17, mgl_exponential_var18, mgl_power_var19,
800 mgl_cosine_var20, mgl_rings_var21, mgl_fan_var22, mgl_blob_var23, mgl_pdj_var24,
801 mgl_fan2_var25, mgl_rings2_var26, mgl_eyefish_var27, mgl_bubble_var28, mgl_cylinder_var29,
802 mgl_perspective_var30, mgl_noise_var31, mgl_juliaN_var32, mgl_juliaScope_var33, mgl_blur_var34,
803 mgl_gaussian_var35, mgl_radialBlur_var36, mgl_pie_var37, mgl_ngon_var38, mgl_curl_var39,
804 mgl_rectangles_var40, mgl_arch_var41, mgl_tangent_var42, mgl_square_var43, mgl_blade_var44,
805 mgl_secant_var45, mgl_rays_var46, mgl_twintrian_var47,mgl_cross_var48, mgl_disc2_var49,
806 mgl_supershape_var50, mgl_flower_var51, mgl_conic_var52, mgl_parabola_var53, mgl_bent2_var54,
807 mgl_bipolar_var55, mgl_boarders_var56, mgl_butterfly_var57,mgl_cell_var58, mgl_cpow_var59,
808 mgl_curve_var60, mgl_edisc_var61, mgl_elliptic_var62, mgl_escher_var63, mgl_foci_var64,
809 mgl_lazySusan_var65, mgl_loonie_var66, mgl_preBlur_var67, mgl_modulus_var68, mgl_oscope_var69,
810 mgl_polar2_var70, mgl_popcorn2_var71, mgl_scry_var72, mgl_separation_var73, mgl_split_var74,
811 mgl_splits_var75, mgl_stripes_var76, mgl_wedge_var77, mgl_wedgeJulia_var78, mgl_wedgeSph_var79,
812 mgl_whorl_var80, mgl_waves2_var81, mgl_exp_var82, mgl_log_var83, mgl_sin_var84,
813 mgl_cos_var85, mgl_tan_var86, mgl_sec_var87, mgl_csc_var88, mgl_cot_var89,
814 mgl_sinh_var90, mgl_cosh_var91, mgl_tanh_var92, mgl_sech_var93, mgl_csch_var94,
815 mgl_coth_var95, mgl_auger_var96, mgl_flux_var97};
816 //-----------------------------------------------------------------------------
mgl_flame_2d_point(HCDT A,HCDT F,mreal & x,mreal & y,mreal amax)817 long static mgl_flame_2d_point(HCDT A, HCDT F, mreal& x, mreal& y, mreal amax)
818 {
819 long i, n=A->GetNy(), m=F->GetNy(), last_func=0, l=F->GetNx();
820 l = l>6?6:l;
821 mreal r = amax*mgl_rnd(), sum_prob = 0, x1, y1;
822 for(i=0;i<n;i++)
823 {
824 sum_prob += A->v(6,i);
825 if(r<sum_prob) break;
826 }
827 x1 = A->v(0,i)*x+A->v(1,i)*y+A->v(4,i);
828 y1 = A->v(2,i)*x+A->v(3,i)*y+A->v(5,i);
829 x = y = 0;
830 for(long j=0;j<m;j++)
831 {
832 int v=int(F->v(0,j,i)+0.5);
833 mreal par[5] = {F->v(1,j,i),0,0,0,0};
834 for(int k=2;k<l;k++) par[k-1]=F->v(k,j,i);
835 if(v<0 || v>=mglFlame2dLAST) { v=0; par[0]=1; }
836 ffunc[v](x,y,x1,y1,par); last_func=v;
837 }
838 return last_func;
839 }
mgl_data_flame_2d(HCDT A,HCDT F,long n,long skip)840 HMDT MGL_EXPORT mgl_data_flame_2d(HCDT A, HCDT F, long n, long skip)
841 {
842 if(!A || A->GetNx()<7 || n<1) return 0; // incompatible dimensions
843 if(!F || F->GetNx()<2 || F->GetNz()!=A->GetNy()) return 0; // incompatible dimensions
844 mreal amax=0;
845 for(long i=0; i<A->GetNy(); i++) amax += A->v(6,i);
846 if(amax<=0) return 0;
847
848 mglData *f = new mglData(3,n);
849 mreal x = 0, y = 0;
850 for(long i=0; i<skip; i++) mgl_flame_2d_point(A, F, x, y,amax);
851 for(long i=0; i<n; i++)
852 {
853 f->a[3*i+2] = mgl_flame_2d_point(A, F, x, y, amax); // TODO color information ?!!
854 f->a[3*i] = x; f->a[3*i+1] = y;
855 }
856 return f;
857 }
mgl_data_flame_2d_(uintptr_t * d,uintptr_t * f,long * n,long * skip)858 uintptr_t MGL_EXPORT mgl_data_flame_2d_(uintptr_t *d, uintptr_t *f, long *n, long *skip)
859 { return uintptr_t(mgl_data_flame_2d(_DT_,_DA_(f),*n,*skip)); }
860 //-----------------------------------------------------------------------------
861