Ocean Sim modifier patch
[blender.git] / source / blender / blenkernel / intern / ocean.c
1 /* 
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * Contributors: Matt Ebb, Hamed Zaghaghi
22  * Based on original code by Drew Whitehouse / Houdini Ocean Toolkit
23  * OpenMP hints by Christian Schnellhammer
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28
29 #include <math.h>
30 #include <stdlib.h>
31
32 #include <string.h>
33
34 #include "MEM_guardedalloc.h"
35
36 #include "DNA_scene_types.h"
37
38 #include "BKE_image.h"
39 #include "BKE_ocean.h"
40 #include "BKE_utildefines.h"
41
42 #include "BKE_global.h" // XXX TESTING
43
44 #include "BLI_math_base.h"
45 #include "BLI_math_inline.h"
46 #include "BLI_rand.h"
47 #include "BLI_string.h"
48 #include "BLI_threads.h"
49 #include "BLI_utildefines.h"
50
51 #include "IMB_imbuf.h"
52 #include "IMB_imbuf_types.h"
53
54 #include "RE_render_ext.h"
55
56 #ifdef WITH_OCEANSIM
57
58 // Ocean code
59 #include "fftw3.h"
60
61 #define GRAVITY  9.81f
62
63 typedef struct Ocean {
64         /* ********* input parameters to the sim ********* */
65     float _V;
66     float _l;
67     float _w;
68     float _A;
69     float _damp_reflections;
70     float _wind_alignment;
71     float _depth;
72         
73     float _wx;
74     float _wz;
75         
76     float _L;
77         
78     /* dimensions of computational grid */
79     int _M;
80     int _N;
81         
82     /* spatial size of computational grid */
83     float _Lx;
84     float _Lz;
85         
86     float normalize_factor;                                     // init w
87     float time;
88
89     short _do_disp_y;
90     short _do_normals;
91     short _do_chop;
92     short _do_jacobian;
93         
94         /* mutex for threaded texture access */
95         ThreadRWMutex oceanmutex;
96         
97         /* ********* sim data arrays ********* */
98         
99     /* two dimensional arrays of complex */
100     fftw_complex *_fft_in;                      // init w       sim w
101         fftw_complex *_fft_in_x;                        // init w       sim w
102         fftw_complex *_fft_in_z;                        // init w       sim w
103         fftw_complex *_fft_in_jxx;                      // init w       sim w
104         fftw_complex *_fft_in_jzz;                      // init w       sim w
105         fftw_complex *_fft_in_jxz;                      // init w       sim w
106         fftw_complex *_fft_in_nx;                       // init w       sim w
107         fftw_complex *_fft_in_nz;                       // init w       sim w
108     fftw_complex *_htilda;                      // init w       sim w (only once)
109         
110     /* fftw "plans" */
111     fftw_plan _disp_y_plan;                     // init w       sim r
112     fftw_plan _disp_x_plan;                     // init w       sim r
113     fftw_plan _disp_z_plan;                     // init w       sim r
114     fftw_plan _N_x_plan;                        // init w       sim r
115     fftw_plan _N_z_plan;                        // init w       sim r
116     fftw_plan _Jxx_plan;                        // init w       sim r
117     fftw_plan _Jxz_plan;                        // init w       sim r
118     fftw_plan _Jzz_plan;                        // init w       sim r
119         
120     /* two dimensional arrays of float */
121     double *  _disp_y;                          // init w       sim w via plan?
122     double * _N_x;                                      // init w       sim w via plan?
123     /*float * _N_y; all member of this array has same values, so convert this array to a float to reduce memory usage (MEM01)*/
124     double _N_y;                                        //                      sim w ********* can be rearranged?
125     double * _N_z;                                      // init w       sim w via plan?
126     double * _disp_x;                           // init w       sim w via plan?
127     double * _disp_z;                           // init w       sim w via plan?
128         
129     /* two dimensional arrays of float */
130     /* Jacobian and minimum eigenvalue */
131     double * _Jxx;                                      // init w       sim w
132     double * _Jzz;                                      // init w       sim w
133     double * _Jxz;                                      // init w       sim w
134         
135     /* one dimensional float array */
136     float * _kx;                                        // init w       sim r
137     float * _kz;                                        // init w       sim r
138         
139     /* two dimensional complex array */
140     fftw_complex * _h0;                         // init w       sim r
141     fftw_complex * _h0_minus;           // init w       sim r
142     
143     /* two dimensional float array */
144     float * _k;                                         // init w       sim r
145 } Ocean;
146
147
148
149 static float nextfr(float min, float max)
150 {
151         return BLI_frand()*(min-max)+max; 
152 }
153
154 static float gaussRand (void)
155 {
156     float x;            // Note: to avoid numerical problems with very small
157     float y;            // numbers, we make these variables singe-precision
158     float length2;      // floats, but later we call the double-precision log()
159         // and sqrt() functions instead of logf() and sqrtf().
160     do
161     {
162                 x = (float) (nextfr (-1, 1));
163                 y = (float)(nextfr (-1, 1));
164                 length2 = x * x + y * y;
165     }
166     while (length2 >= 1 || length2 == 0);
167         
168     return x * sqrt (-2 * log (length2) / length2);
169 }
170
171 /**
172  * Som usefull functions
173  * */
174 MINLINE float lerp(float a,float b,float f) 
175 {
176         return a + (b-a)*f; 
177 }
178
179 MINLINE float catrom(float p0,float p1,float p2,float p3,float f) 
180
181     return 0.5 *((2 * p1) +
182                  (-p0 + p2) * f +
183                  (2*p0 - 5*p1 + 4*p2 - p3) * f*f +
184                  (-p0 + 3*p1- 3*p2 + p3) * f*f*f);
185 }
186
187 MINLINE float omega(float k, float depth) 
188
189     return sqrt(GRAVITY*k * tanh(k*depth));
190 }
191
192 // modified Phillips spectrum 
193 static float Ph(struct Ocean* o, float kx,float kz ) 
194 {
195         float tmp;
196     float k2 = kx*kx + kz*kz;
197         
198     if (k2 == 0.0)
199     {
200         return 0.0; // no DC component
201     }
202         
203     // damp out the waves going in the direction opposite the wind
204     tmp = (o->_wx * kx  + o->_wz * kz)/sqrt(k2);
205     if (tmp < 0) 
206     {
207         tmp *= o->_damp_reflections;
208     }
209
210     return o->_A * exp( -1.0f / (k2*(o->_L*o->_L))) * exp(-k2 * (o->_l*o->_l)) * pow(fabs(tmp),o->_wind_alignment) / (k2*k2);
211 }
212
213 static void compute_eigenstuff(struct OceanResult *ocr, float jxx,float jzz,float jxz)
214 {
215     float a,b,qplus,qminus;
216     a = jxx + jzz; 
217     b = sqrt((jxx - jzz)*(jxx - jzz) + 4 * jxz * jxz);
218         
219     ocr->Jminus = 0.5*(a-b);
220     ocr->Jplus  = 0.5*(a+b);
221         
222     qplus  = (ocr->Jplus  - jxx)/jxz;
223     qminus = (ocr->Jminus - jxx)/jxz;
224         
225     a = sqrt(1 + qplus*qplus);
226     b = sqrt(1 + qminus*qminus);
227         
228     ocr->Eplus[0] = 1.0/ a;
229     ocr->Eplus[1] = 0.0;
230     ocr->Eplus[2] = qplus/a;
231         
232     ocr->Eminus[0] = 1.0/b;
233     ocr->Eminus[1] = 0.0;
234     ocr->Eminus[2] = qminus/b;  
235 }
236
237 /*
238  * instead of Complex.h 
239  * in fftw.h "fftw_complex" typedefed as double[2]
240  * below you can see functions are needed to work with such complex numbers.
241  * */
242 static void init_complex(fftw_complex cmpl, float real, float image)
243 {
244         cmpl[0] = real;
245         cmpl[1] = image;
246 }
247
248 #if 0   // unused
249 static void add_complex_f(fftw_complex res, fftw_complex cmpl, float f)
250 {
251         res[0] = cmpl[0] + f;
252         res[1] = cmpl[1];
253 }
254 #endif
255
256 static void add_comlex_c(fftw_complex res, fftw_complex cmpl1, fftw_complex cmpl2)
257 {
258         res[0] = cmpl1[0] + cmpl2[0];
259         res[1] = cmpl1[1] + cmpl2[1];
260 }
261
262 static void mul_complex_f(fftw_complex res, fftw_complex cmpl, float f)
263 {
264         res[0] = cmpl[0]*f;
265         res[1] = cmpl[1]*f;
266 }
267
268 static void mul_complex_c(fftw_complex res, fftw_complex cmpl1, fftw_complex cmpl2)
269 {
270         fftwf_complex temp;
271         temp[0] = cmpl1[0]*cmpl2[0]-cmpl1[1]*cmpl2[1];
272         temp[1] = cmpl1[0]*cmpl2[1]+cmpl1[1]*cmpl2[0];
273         res[0] = temp[0];
274         res[1] = temp[1];
275 }
276
277 static float real_c(fftw_complex cmpl)
278 {
279         return cmpl[0];
280 }
281
282 static float image_c(fftw_complex cmpl)
283 {
284         return cmpl[1];
285 }
286
287 static void conj_complex(fftw_complex res, fftw_complex cmpl1)
288 {
289         res[0] = cmpl1[0];
290         res[1] = -cmpl1[1];
291 }
292
293 static void exp_complex(fftw_complex res, fftw_complex cmpl)
294 {
295         float r = expf(cmpl[0]);
296         
297         res[0] = cos(cmpl[1])*r;
298         res[1] = sin(cmpl[1])*r;
299 }
300
301 float BKE_ocean_jminus_to_foam(float jminus, float coverage) {
302         float foam = jminus * -0.005 + coverage;
303         CLAMP(foam, 0.0, 1.0);
304         return foam*foam;
305 }
306
307 void BKE_ocean_eval_uv(struct Ocean *oc, struct OceanResult *ocr, float u,float v)
308 {
309     int i0,i1,j0,j1;
310     float frac_x,frac_z;
311     float uu,vv;
312         
313     // first wrap the texture so 0 <= (u,v) < 1
314     u = fmod(u,1.0f);
315     v = fmod(v,1.0f);
316         
317     if (u < 0) u += 1.0f;
318     if (v < 0) v += 1.0f;
319         
320     BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
321         
322         uu = u * oc->_M;
323     vv = v * oc->_N;
324         
325     i0 = (int)floor(uu);
326     j0 = (int)floor(vv);
327         
328     i1 = (i0 + 1);
329     j1 = (j0 + 1);
330         
331     frac_x = uu - i0;
332     frac_z = vv - j0;
333         
334     i0 = i0 % oc->_M;
335     j0 = j0 % oc->_N;
336         
337     i1 = i1 % oc->_M;
338     j1 = j1 % oc->_N;
339         
340         
341 #define BILERP(m) (lerp(lerp(m[i0*oc->_N+j0],m[i1*oc->_N+j0],frac_x),lerp(m[i0*oc->_N+j1],m[i1*oc->_N+j1],frac_x),frac_z))
342     {
343         if (oc->_do_disp_y) {
344                 ocr->disp[1] = BILERP(oc->_disp_y);
345         }
346                 
347         if (oc->_do_normals) {
348                 ocr->normal[0] = BILERP(oc->_N_x);
349                 ocr->normal[1] = oc->_N_y/*BILERP(oc->_N_y) (MEM01)*/;
350                 ocr->normal[2] = BILERP(oc->_N_z);
351         }
352                 
353         if (oc->_do_chop) {
354                 ocr->disp[0] = BILERP(oc->_disp_x);
355                 ocr->disp[2] = BILERP(oc->_disp_z); 
356         } else {
357                 ocr->disp[0] = 0.0;
358                 ocr->disp[2] = 0.0;
359         }
360                 
361         if (oc->_do_jacobian) {
362             compute_eigenstuff(ocr, BILERP(oc->_Jxx),BILERP(oc->_Jzz),BILERP(oc->_Jxz));
363         }                      
364     }
365 #undef BILERP
366         
367         BLI_rw_mutex_unlock(&oc->oceanmutex);
368 }
369
370 // use catmullrom interpolation rather than linear
371 void BKE_ocean_eval_uv_catrom(struct Ocean *oc, struct OceanResult *ocr, float u,float v)
372 {
373     int i0,i1,i2,i3,j0,j1,j2,j3;
374     float frac_x,frac_z;
375     float uu,vv;
376         
377     // first wrap the texture so 0 <= (u,v) < 1
378     u = fmod(u,1.0f);
379     v = fmod(v,1.0f);
380         
381     if (u < 0) u += 1.0f;
382     if (v < 0) v += 1.0f;
383         
384         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
385         
386     uu = u * oc->_M;
387     vv = v * oc->_N;
388         
389     i1 = (int)floor(uu);
390     j1 = (int)floor(vv);
391         
392     i2 = (i1 + 1);
393     j2 = (j1 + 1);
394         
395     frac_x = uu - i1;
396     frac_z = vv - j1;
397         
398     i1 = i1 % oc->_M;
399     j1 = j1 % oc->_N;
400         
401     i2 = i2 % oc->_M;
402     j2 = j2 % oc->_N;
403         
404     i0 = (i1-1);
405     i3 = (i2+1);
406     i0 = i0 <   0 ? i0 + oc->_M : i0;
407     i3 = i3 >= oc->_M ? i3 - oc->_M : i3;
408         
409     j0 = (j1-1);
410     j3 = (j2+1);
411     j0 = j0 <   0 ? j0 + oc->_N : j0;
412     j3 = j3 >= oc->_N ? j3 - oc->_N : j3;
413         
414 #define INTERP(m) catrom(catrom(m[i0*oc->_N+j0],m[i1*oc->_N+j0],m[i2*oc->_N+j0],m[i3*oc->_N+j0],frac_x),\
415 catrom(m[i0*oc->_N+j1],m[i1*oc->_N+j1],m[i2*oc->_N+j1],m[i3*oc->_N+j1],frac_x),\
416 catrom(m[i0*oc->_N+j2],m[i1*oc->_N+j2],m[i2*oc->_N+j2],m[i3*oc->_N+j2],frac_x),\
417 catrom(m[i0*oc->_N+j3],m[i1*oc->_N+j3],m[i2*oc->_N+j3],m[i3*oc->_N+j3],frac_x),\
418 frac_z)
419         
420     {
421         if (oc->_do_disp_y)
422         {
423                 ocr->disp[1] = INTERP(oc->_disp_y) ;
424         }
425         if (oc->_do_normals)
426         {
427                 ocr->normal[0] = INTERP(oc->_N_x);
428                 ocr->normal[1] = oc->_N_y/*INTERP(oc->_N_y) (MEM01)*/;
429                 ocr->normal[2] = INTERP(oc->_N_z);
430         }
431         if (oc->_do_chop) 
432         {
433                 ocr->disp[0] = INTERP(oc->_disp_x);
434                 ocr->disp[2] = INTERP(oc->_disp_z); 
435         }
436         else
437         {
438                 ocr->disp[0] = 0.0;
439                 ocr->disp[2] = 0.0;
440         }
441                 
442         if (oc->_do_jacobian)
443         {
444             compute_eigenstuff(ocr, INTERP(oc->_Jxx),INTERP(oc->_Jzz),INTERP(oc->_Jxz));
445         }                      
446     }
447 #undef INTERP
448         
449         BLI_rw_mutex_unlock(&oc->oceanmutex);
450         
451 }
452
453 void BKE_ocean_eval_xz(struct Ocean *oc, struct OceanResult *ocr, float x,float z)
454 {
455     BKE_ocean_eval_uv(oc, ocr, x/oc->_Lx,z/oc->_Lz);
456 }
457
458 void BKE_ocean_eval_xz_catrom(struct Ocean *oc, struct OceanResult *ocr, float x,float z)
459 {
460     BKE_ocean_eval_uv_catrom(oc, ocr, x/oc->_Lx,z/oc->_Lz);
461 }
462
463 // note that this doesn't wrap properly for i,j < 0, but its
464 // not really meant for that being just a way to get the raw data out
465 // to save in some image format.
466 void BKE_ocean_eval_ij(struct Ocean *oc, struct OceanResult *ocr, int i,int j)
467 {
468         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
469         
470     i = abs(i) % oc->_M;
471     j = abs(j) % oc->_N;
472         
473     ocr->disp[1] = oc->_do_disp_y ? oc->_disp_y[i*oc->_N+j] : 0.0f;
474         
475     if (oc->_do_chop)
476     {
477         ocr->disp[0] = oc->_disp_x[i*oc->_N+j];
478         ocr->disp[2] = oc->_disp_z[i*oc->_N+j];
479     }
480     else
481     {
482         ocr->disp[0] = 0.0f;
483         ocr->disp[2] = 0.0f;
484     }
485         
486     if (oc->_do_normals)
487     {
488         ocr->normal[0] = oc->_N_x[i*oc->_N+j];
489         ocr->normal[1] = oc->_N_y/*oc->_N_y[i*oc->_N+j] (MEM01)*/;
490         ocr->normal[2] = oc->_N_z[i*oc->_N+j];
491     }
492         
493     if (oc->_do_jacobian)
494     {
495                 compute_eigenstuff(ocr, oc->_Jxx[i*oc->_N+j],oc->_Jzz[i*oc->_N+j],oc->_Jxz[i*oc->_N+j]);
496     }
497         
498         BLI_rw_mutex_unlock(&oc->oceanmutex);
499 }
500
501 void BKE_simulate_ocean(struct Ocean *o, float t, float scale, float chop_amount)
502 {
503         int i, j;
504         
505         scale *= o->normalize_factor;
506                 
507         BLI_rw_mutex_lock(&o->oceanmutex, THREAD_LOCK_WRITE);
508         
509         // compute a new htilda
510         #pragma omp parallel for private(i, j)
511     for (i = 0 ; i  < o->_M ; ++i)
512     {
513         // note the <= _N/2 here, see the fftw doco about
514         // the mechanics of the complex->real fft storage
515         for ( j  = 0 ; j  <= o->_N / 2 ; ++j)
516         {
517                 fftw_complex exp_param1;
518                 fftw_complex exp_param2;
519                 fftw_complex conj_param;
520                 
521                 
522                 init_complex(exp_param1, 0.0, omega(o->_k[i*(1+o->_N/2)+j],o->_depth)*t);
523                 init_complex(exp_param2, 0.0, -omega(o->_k[i*(1+o->_N/2)+j],o->_depth)*t);
524                 exp_complex(exp_param1, exp_param1);
525                 exp_complex(exp_param2, exp_param2);
526                 conj_complex(conj_param, o->_h0_minus[i*o->_N+j]);
527                 
528                 mul_complex_c(exp_param1, o->_h0[i*o->_N+j], exp_param1);
529                 mul_complex_c(exp_param2, conj_param, exp_param2);
530                         
531                 add_comlex_c(o->_htilda[i*(1+o->_N/2)+j], exp_param1, exp_param2);
532                 mul_complex_f(o->_fft_in[i*(1+o->_N/2)+j], o->_htilda[i*(1+o->_N/2)+j], scale);
533         }
534     }
535         
536         #pragma omp parallel sections private(i, j)
537         {
538                 
539         #pragma omp section
540         {
541                 if (o->_do_disp_y)
542                 {
543                         // y displacement
544                         fftw_execute(o->_disp_y_plan);
545                 }
546         } // section 1
547                 
548         #pragma omp section
549         {
550                 if (o->_do_chop)
551                 {
552                         // x displacement
553                         for ( i = 0 ; i  < o->_M ; ++i)
554                         {   
555                                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j)
556                                 {     
557                                         fftw_complex mul_param;
558                                         fftw_complex minus_i;
559                                         
560                                         init_complex(minus_i, 0.0, -1.0);
561                                         init_complex(mul_param, -scale, 0);
562                                         mul_complex_f(mul_param, mul_param, chop_amount);
563                                         mul_complex_c(mul_param, mul_param, minus_i);
564                                         mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
565                                         mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kx[i] / o->_k[i*(1+o->_N/2)+j]));
566                                         init_complex(o->_fft_in_x[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
567                                 }
568                         }
569                         fftw_execute(o->_disp_x_plan);
570                 }
571         } //section 2
572                 
573         #pragma omp section
574         {
575                 if (o->_do_chop)
576                 {       
577                         // z displacement
578                         for ( i = 0 ; i  < o->_M ; ++i)
579                         {   
580                                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j) 
581                                 {               
582                                         fftw_complex mul_param;
583                                         fftw_complex minus_i;
584                                         
585                                         init_complex(minus_i, 0.0, -1.0);
586                                         init_complex(mul_param, -scale, 0);
587                                         mul_complex_f(mul_param, mul_param, chop_amount);
588                                         mul_complex_c(mul_param, mul_param, minus_i);
589                                         mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
590                                         mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kz[j] / o->_k[i*(1+o->_N/2)+j]));
591                                         init_complex(o->_fft_in_z[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));                             
592                                 }
593                         } 
594                         fftw_execute(o->_disp_z_plan);
595                 }
596         } // section 3
597
598         #pragma omp section
599         {
600                 if (o->_do_jacobian)
601                 {
602                         // Jxx
603                         for ( i = 0 ; i  < o->_M ; ++i)
604                         {
605                                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j)
606                                 {    
607                                         fftw_complex mul_param;
608                                         
609                                         //init_complex(mul_param, -scale, 0);
610                                         init_complex(mul_param, -1, 0);
611                                         
612                                         mul_complex_f(mul_param, mul_param, chop_amount);
613                                         mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
614                                         mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kx[i]*o->_kx[i] / o->_k[i*(1+o->_N/2)+j]));
615                                         init_complex(o->_fft_in_jxx[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
616                                 }
617                         }
618                         fftw_execute(o->_Jxx_plan);
619                         
620                         for ( i = 0 ; i  < o->_M ; ++i)
621                         {   
622                                 for ( j  = 0 ; j  < o->_N ; ++j) 
623                                 {
624                                         o->_Jxx[i*o->_N+j] += 1.0;
625                                 }
626                         }
627                 }
628         } // section 4
629         
630         #pragma omp section
631         {
632                 if (o->_do_jacobian)
633                 {
634                         // Jzz
635                         for ( i = 0 ; i  < o->_M ; ++i)
636                         {
637                                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j)
638                                 {   
639                                         fftw_complex mul_param;
640                                         
641                                         //init_complex(mul_param, -scale, 0);
642                                         init_complex(mul_param, -1, 0);
643                                         
644                                         mul_complex_f(mul_param, mul_param, chop_amount);
645                                         mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
646                                         mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kz[j]*o->_kz[j] / o->_k[i*(1+o->_N/2)+j]));
647                                         init_complex(o->_fft_in_jzz[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
648                                 }
649                         }
650                         fftw_execute(o->_Jzz_plan);
651                         for ( i = 0 ; i  < o->_M ; ++i)
652                         {   
653                                 for ( j  = 0 ; j  < o->_N ; ++j) 
654                                 {
655                                         o->_Jzz[i*o->_N+j] += 1.0;
656                                 }
657                         }
658                 }
659         } // section 5
660                 
661         #pragma omp section
662         {       
663                 if (o->_do_jacobian)
664                 {
665                         // Jxz
666                         for ( i = 0 ; i  < o->_M ; ++i)
667                         {
668                                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j)
669                                 {     
670                                         fftw_complex mul_param;
671                                         
672                                         //init_complex(mul_param, -scale, 0);
673                                         init_complex(mul_param, -1, 0);
674
675                                         mul_complex_f(mul_param, mul_param, chop_amount);
676                                         mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
677                                         mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kx[i]*o->_kz[j] / o->_k[i*(1+o->_N/2)+j]));
678                                         init_complex(o->_fft_in_jxz[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
679                                 }
680                         }
681                         fftw_execute(o->_Jxz_plan);
682                 }
683         } // section 6
684                 
685         #pragma omp section
686         {
687                 // fft normals
688                 if (o->_do_normals)
689                 {
690                         for ( i = 0 ; i  < o->_M ; ++i)
691                         {
692                                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j)
693                                 {     
694                                         fftw_complex mul_param;
695                                         
696                                         init_complex(mul_param, 0.0, -1.0);
697                                         mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
698                                         mul_complex_f(mul_param, mul_param, o->_kx[i]);
699                                         init_complex(o->_fft_in_nx[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
700                                 }
701                         }
702                         fftw_execute(o->_N_x_plan);
703                 
704                 }
705         } // section 7
706         
707         #pragma omp section
708         {
709                 if (o->_do_normals)
710                 {       
711                         for ( i = 0 ; i  < o->_M ; ++i)
712                         {   
713                                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j) 
714                                 {       
715                                         fftw_complex mul_param;
716                                         
717                                         init_complex(mul_param, 0.0, -1.0);
718                                         mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
719                                         mul_complex_f(mul_param, mul_param, o->_kz[i]);
720                                         init_complex(o->_fft_in_nz[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
721                                 }
722                         } 
723                         fftw_execute(o->_N_z_plan);
724                         
725                         /*for ( i = 0 ; i  < o->_M ; ++i)
726                          {   
727                          for ( j  = 0 ; j  < o->_N ; ++j) 
728                          {
729                          o->_N_y[i*o->_N+j] = 1.0f/scale;
730                          }
731                          }
732                          (MEM01)*/
733                         o->_N_y = 1.0f/scale;        
734                 }
735         } // section 8
736         
737         } // omp sections
738         
739         BLI_rw_mutex_unlock(&o->oceanmutex);
740 }
741
742 static void set_height_normalize_factor(struct Ocean *oc) 
743 {  
744     float res = 1.0;
745     float max_h = 0.0;
746         
747         int i,j;
748         
749         if (!oc->_do_disp_y) return;
750         
751         oc->normalize_factor = 1.0;
752         
753         BKE_simulate_ocean(oc, 0.0, 1.0, 0);
754         
755         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
756         
757     for (i = 0; i < oc->_M; ++i)
758     {
759         for (j = 0; j < oc->_N; ++j)
760         {   
761             if( max_h < fabsf(oc->_disp_y[i*oc->_N+j]))
762             {
763                 max_h = fabsf(oc->_disp_y[i*oc->_N+j]);
764             }
765         } 
766     }
767         
768         BLI_rw_mutex_unlock(&oc->oceanmutex);
769         
770     if (max_h == 0.0) max_h = 0.00001f; // just in case ...
771         
772     res = 1.0f / (max_h);
773
774         oc->normalize_factor = res;
775 }
776
777 struct Ocean *BKE_add_ocean(void)
778 {
779         Ocean *oc = MEM_callocN(sizeof(Ocean), "ocean sim data");
780
781         BLI_rw_mutex_init(&oc->oceanmutex);
782         
783         return oc;
784 }
785
786 void BKE_init_ocean(struct Ocean* o, int M,int N, float Lx, float Lz, float V, float l, float A, float w, float damp, 
787                                            float alignment, float depth, float time, short do_height_field, short do_chop, short do_normals, short do_jacobian, int seed)
788 {
789         int i,j,ii;
790         
791         BLI_rw_mutex_lock(&o->oceanmutex, THREAD_LOCK_WRITE);
792         
793         o->_M = M;
794     o->_N = N;
795     o->_V = V;
796     o->_l = l;
797     o->_A = A;
798     o->_w = w;
799     o->_damp_reflections = 1.0 - damp;
800     o->_wind_alignment = alignment;
801     o->_depth = depth;
802     o->_Lx = Lx;
803     o->_Lz = Lz;
804     o->_wx = cos(w);
805     o->_wz = -sin(w); // wave direction
806     o->_L = V*V / GRAVITY;  // largest wave for a given velocity V
807     o->time = time;
808         
809     o->_do_disp_y = do_height_field;
810     o->_do_normals = do_normals;
811     o->_do_chop = do_chop;
812     o->_do_jacobian = do_jacobian;
813     
814     o->_k = (float*) MEM_mallocN(M * (1+N/2) * sizeof(float), "ocean_k");
815     o->_h0 = (fftw_complex*) MEM_mallocN(M * N * sizeof(fftw_complex), "ocean_h0");
816     o->_h0_minus = (fftw_complex*) MEM_mallocN(M * N * sizeof(fftw_complex), "ocean_h0_minus");
817     o->_kx = (float*) MEM_mallocN(o->_M * sizeof(float), "ocean_kx");
818     o->_kz = (float*) MEM_mallocN(o->_N * sizeof(float), "ocean_kz");
819         
820     // make this robust in the face of erroneous usage
821     if (o->_Lx == 0.0)
822         o->_Lx = 0.001;
823         
824     if (o->_Lz == 0.0)
825         o->_Lz = 0.001;
826         
827     // the +ve components and DC
828     for (i = 0 ; i <= o->_M/2 ; ++i)
829         o->_kx[i] = 2.0f * M_PI * i / o->_Lx;
830         
831     // the -ve components
832     for (i = o->_M-1,ii=0 ; i > o->_M/2 ; --i,++ii)
833         o->_kx[i] = -2.0f * M_PI * ii / o->_Lx;
834         
835     // the +ve components and DC
836     for (i = 0 ; i <= o->_N/2 ; ++i)
837         o->_kz[i] = 2.0f * M_PI * i / o->_Lz;
838         
839     // the -ve components
840     for (i = o->_N-1,ii=0 ; i > o->_N/2 ; --i,++ii)
841         o->_kz[i] = -2.0f * M_PI * ii / o->_Lz;
842         
843     // pre-calculate the k matrix
844     for (i = 0 ; i  < o->_M ; ++i)
845         for (j  = 0 ; j  <= o->_N / 2 ; ++j) 
846                 o->_k[i*(1+o->_N/2)+j] = sqrt(o->_kx[i]*o->_kx[i] + o->_kz[j]*o->_kz[j] );
847     
848     /*srand(seed);*/
849     BLI_srand(seed);
850         
851     for (i = 0 ; i  < o->_M ; ++i)
852     {
853         for (j = 0 ; j  < o->_N ; ++j)
854         {
855                 float r1 = gaussRand();
856                 float r2 = gaussRand();
857                 
858             fftw_complex r1r2;
859             init_complex(r1r2, r1, r2);
860             mul_complex_f(o->_h0[i*o->_N+j], r1r2, (float)(sqrt(Ph(o,  o->_kx[i], o->_kz[j]) / 2.0f)));
861             mul_complex_f(o->_h0_minus[i*o->_N+j], r1r2, (float)(sqrt(Ph(o, -o->_kx[i],-o->_kz[j]) / 2.0f)));
862         }
863     }
864         
865     o->_fft_in = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in");
866     o->_htilda = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_htilda");
867     
868     if (o->_do_disp_y){
869         o->_disp_y = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_disp_y");
870         o->_disp_y_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in, o->_disp_y, FFTW_ESTIMATE);
871     }
872     
873     if (o->_do_normals){
874                 o->_fft_in_nx = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_nx");
875                 o->_fft_in_nz = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_nz");
876                 
877         o->_N_x = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_N_x");
878         /*o->_N_y = (float*) fftwf_malloc(o->_M * o->_N * sizeof(float)); (MEM01)*/
879         o->_N_z = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_N_z");
880                 
881         o->_N_x_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in_nx, o->_N_x, FFTW_ESTIMATE);
882         o->_N_z_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in_nz, o->_N_z, FFTW_ESTIMATE);
883     }
884     
885     if (o->_do_chop){
886                 o->_fft_in_x = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_x");
887                 o->_fft_in_z = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_z");
888                 
889         o->_disp_x = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_disp_x");
890         o->_disp_z = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_disp_z");
891                 
892         o->_disp_x_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in_x, o->_disp_x, FFTW_ESTIMATE);
893         o->_disp_z_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in_z, o->_disp_z, FFTW_ESTIMATE);
894     }
895     if (o->_do_jacobian){
896                 o->_fft_in_jxx = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_jxx");
897                 o->_fft_in_jzz = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_jzz");
898                 o->_fft_in_jxz = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_jxz");
899                 
900         o->_Jxx = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_Jxx");
901         o->_Jzz = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_Jzz");
902         o->_Jxz = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_Jxz");
903                 
904         o->_Jxx_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in_jxx, o->_Jxx, FFTW_ESTIMATE);
905         o->_Jzz_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in_jzz, o->_Jzz, FFTW_ESTIMATE);
906         o->_Jxz_plan = fftw_plan_dft_c2r_2d(o->_M,o->_N, o->_fft_in_jxz, o->_Jxz, FFTW_ESTIMATE);
907     }
908         
909         BLI_rw_mutex_unlock(&o->oceanmutex);
910         
911     set_height_normalize_factor(o);
912         
913 }  
914
915 void BKE_free_ocean_data(struct Ocean *oc) 
916 {
917         if(!oc) return;
918         
919         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_WRITE);
920         
921         if (oc->_do_disp_y)
922         {
923                 fftw_destroy_plan(oc->_disp_y_plan);
924                 MEM_freeN(oc->_disp_y);
925         }
926         
927         if (oc->_do_normals)
928         {
929                 MEM_freeN(oc->_fft_in_nx);
930                 MEM_freeN(oc->_fft_in_nz);
931                 fftw_destroy_plan(oc->_N_x_plan);
932                 fftw_destroy_plan(oc->_N_z_plan);
933                 MEM_freeN(oc->_N_x);
934                 /*fftwf_free(oc->_N_y); (MEM01)*/
935                 MEM_freeN(oc->_N_z);
936         }
937         
938         if (oc->_do_chop)
939         {
940                 MEM_freeN(oc->_fft_in_x);
941                 MEM_freeN(oc->_fft_in_z);
942                 fftw_destroy_plan(oc->_disp_x_plan);
943                 fftw_destroy_plan(oc->_disp_z_plan);
944                 MEM_freeN(oc->_disp_x);
945                 MEM_freeN(oc->_disp_z);
946         }
947         
948         if (oc->_do_jacobian)
949         {
950                 MEM_freeN(oc->_fft_in_jxx);
951                 MEM_freeN(oc->_fft_in_jzz);
952                 MEM_freeN(oc->_fft_in_jxz);
953                 fftw_destroy_plan(oc->_Jxx_plan);
954                 fftw_destroy_plan(oc->_Jzz_plan);
955                 fftw_destroy_plan(oc->_Jxz_plan);
956                 MEM_freeN(oc->_Jxx);
957                 MEM_freeN(oc->_Jzz);
958                 MEM_freeN(oc->_Jxz);
959         }
960         
961         if (oc->_fft_in)
962                 MEM_freeN(oc->_fft_in);
963         
964         /* check that ocean data has been initialised */
965         if (oc->_htilda) {
966                 MEM_freeN(oc->_htilda);
967                 MEM_freeN(oc->_k);
968                 MEM_freeN(oc->_h0);
969                 MEM_freeN(oc->_h0_minus);
970                 MEM_freeN(oc->_kx);
971                 MEM_freeN(oc->_kz);
972         }
973         
974         BLI_rw_mutex_unlock(&oc->oceanmutex);
975 }
976
977 void BKE_free_ocean(struct Ocean *oc) 
978 {
979         if(!oc) return;
980         
981         BKE_free_ocean_data(oc);
982         BLI_rw_mutex_end(&oc->oceanmutex);
983         
984         MEM_freeN(oc);
985 }
986
987 #undef GRAVITY
988
989
990 /* ********* Baking/Caching ********* */
991
992
993 #define CACHE_TYPE_DISPLACE     1
994 #define CACHE_TYPE_FOAM         2
995 #define CACHE_TYPE_NORMAL       3
996
997 static void cache_filename(char *string, char *path, int frame, int type)
998 {
999         char *cachepath=NULL;
1000         
1001         switch(type) {
1002                 case CACHE_TYPE_FOAM:
1003                         cachepath = BLI_strdupcat(path, "foam_");
1004                         break;
1005                 case CACHE_TYPE_NORMAL:
1006                         cachepath = BLI_strdupcat(path, "normal_");
1007                         break;
1008                 case CACHE_TYPE_DISPLACE:
1009                 default:
1010                         cachepath = BLI_strdupcat(path, "disp_");
1011                         break;
1012         }               
1013         
1014         BKE_makepicstring(string, cachepath, frame, R_OPENEXR, 1, TRUE);
1015         
1016         MEM_freeN(cachepath);
1017 }
1018
1019 void BKE_free_ocean_cache(struct OceanCache *och)
1020 {
1021         int i, f=0;
1022         
1023         if (!och) return;
1024         
1025         if (och->ibufs_disp) {
1026                 for (i=och->start, f=0; i<=och->end; i++, f++)
1027                 {
1028                         if (och->ibufs_disp[f]) {
1029                                 IMB_freeImBuf(och->ibufs_disp[f]);
1030                         }
1031                 }
1032                 MEM_freeN(och->ibufs_disp);
1033         }
1034         
1035         if (och->ibufs_foam) {
1036                 for (i=och->start, f=0; i<=och->end; i++, f++)
1037                 {
1038                         if (och->ibufs_foam[f]) {
1039                                 IMB_freeImBuf(och->ibufs_foam[f]);
1040                         }
1041                 }
1042                 MEM_freeN(och->ibufs_foam);
1043         }
1044         
1045         if (och->ibufs_norm) {
1046                 for (i=och->start, f=0; i<=och->end; i++, f++)
1047                 {
1048                         if (och->ibufs_norm[f]) {
1049                                 IMB_freeImBuf(och->ibufs_norm[f]);
1050                         }
1051                 }
1052                 MEM_freeN(och->ibufs_norm);
1053         }
1054         
1055         if (och->time)
1056                 MEM_freeN(och->time);
1057         MEM_freeN(och);
1058 }
1059
1060 void BKE_ocean_cache_eval_uv(struct OceanCache *och, struct OceanResult *ocr, int f, float u, float v)
1061 {
1062         int res_x = och->resolution_x;
1063         int res_y = och->resolution_y;
1064         float result[4];
1065         
1066         u = fmod(u, 1.0);
1067         v = fmod(v, 1.0);
1068         
1069         if (u < 0) u += 1.0f;
1070     if (v < 0) v += 1.0f;
1071         
1072         if (och->ibufs_disp[f]) {
1073                 ibuf_sample(och->ibufs_disp[f], u, v, (1.0/(float)res_x), (1.0/(float)res_y), result);
1074                 ocr->disp[0] = result[0];
1075                 ocr->disp[1] = result[1];
1076                 ocr->disp[2] = result[2];
1077         }
1078         
1079         if (och->ibufs_foam[f]) {
1080                 ibuf_sample(och->ibufs_foam[f], u, v, (1.0/(float)res_x), (1.0/(float)res_y), result);
1081                 ocr->foam = result[0];
1082         }
1083         
1084         if (och->ibufs_norm[f]) {
1085                 ibuf_sample(och->ibufs_norm[f], u, v, (1.0/(float)res_x), (1.0/(float)res_y), result);
1086                 ocr->normal[0] = result[0];
1087                 ocr->normal[1] = result[1];
1088                 ocr->normal[2] = result[2];
1089         }
1090 }
1091
1092 void BKE_ocean_cache_eval_ij(struct OceanCache *och, struct OceanResult *ocr, int f, int i, int j)
1093 {       
1094         int res_x = och->resolution_x;
1095         int res_y = och->resolution_y;
1096         
1097         i = abs(i) % res_x;
1098     j = abs(j) % res_y;
1099
1100         if (och->ibufs_disp[f]) {
1101                 ocr->disp[0] = och->ibufs_disp[f]->rect_float[4*(res_x*j + i) + 0];
1102                 ocr->disp[1] = och->ibufs_disp[f]->rect_float[4*(res_x*j + i) + 1];
1103                 ocr->disp[2] = och->ibufs_disp[f]->rect_float[4*(res_x*j + i) + 2];
1104         }
1105         
1106         if (och->ibufs_foam[f]) {
1107                 ocr->foam = och->ibufs_foam[f]->rect_float[4*(res_x*j + i) + 0];
1108         }
1109         
1110         if (och->ibufs_norm[f]) {
1111                 ocr->normal[0] = och->ibufs_norm[f]->rect_float[4*(res_x*j + i) + 0];
1112                 ocr->normal[1] = och->ibufs_norm[f]->rect_float[4*(res_x*j + i) + 1];
1113                 ocr->normal[2] = och->ibufs_norm[f]->rect_float[4*(res_x*j + i) + 2];
1114         }
1115 }
1116
1117 struct OceanCache *BKE_init_ocean_cache(char *bakepath, int start, int end, float wave_scale, 
1118                                                   float chop_amount, float foam_coverage, float foam_fade, int resolution)
1119 {
1120         OceanCache *och = MEM_callocN(sizeof(OceanCache), "ocean cache data");
1121         
1122         och->bakepath = bakepath;
1123         och->start = start;
1124         och->end = end;
1125         och->duration = (end - start) + 1;
1126         och->wave_scale = wave_scale;
1127         och->chop_amount = chop_amount;
1128         och->foam_coverage = foam_coverage;
1129         och->foam_fade = foam_fade;
1130         och->resolution_x = resolution*resolution;
1131         och->resolution_y = resolution*resolution;
1132
1133         och->ibufs_disp = MEM_callocN(sizeof(ImBuf *)*och->duration, "displacement imbuf pointer array");
1134         och->ibufs_foam = MEM_callocN(sizeof(ImBuf *)*och->duration, "foam imbuf pointer array");
1135         och->ibufs_norm = MEM_callocN(sizeof(ImBuf *)*och->duration, "normal imbuf pointer array");
1136         
1137         och->time = NULL;
1138         
1139         return och;
1140 }
1141
1142 void BKE_simulate_ocean_cache(struct OceanCache *och, int frame)
1143 {
1144         char string[FILE_MAX];
1145         int f = frame;
1146         
1147         /* ibufs array is zero based, but filenames are based on frame numbers */
1148         /* still need to clamp frame numbers to valid range of images on disk though */
1149         CLAMP(frame, och->start, och->end);
1150         f = frame - och->start; // shift to 0 based
1151         
1152         /* if image is already loaded in mem, return */
1153         if (och->ibufs_disp[f] != NULL ) return;
1154         
1155         
1156         cache_filename(string, och->bakepath, frame, CACHE_TYPE_DISPLACE);
1157         och->ibufs_disp[f] = IMB_loadiffname(string, 0);
1158         //if (och->ibufs_disp[f] == NULL) printf("error loading %s \n", string);
1159         //else printf("loaded cache %s \n", string);
1160         
1161         cache_filename(string, och->bakepath, frame, CACHE_TYPE_FOAM);
1162         och->ibufs_foam[f] = IMB_loadiffname(string, 0);
1163         //if (och->ibufs_foam[f] == NULL) printf("error loading %s \n", string);
1164         //else printf("loaded cache %s \n", string);
1165         
1166         cache_filename(string, och->bakepath, frame, CACHE_TYPE_NORMAL);
1167         och->ibufs_norm[f] = IMB_loadiffname(string, 0);
1168         //if (och->ibufs_norm[f] == NULL) printf("error loading %s \n", string);
1169         //else printf("loaded cache %s \n", string);
1170 }
1171
1172
1173 void BKE_bake_ocean(struct Ocean *o, struct OceanCache *och, void (*update_cb)(void *, float progress, int *cancel), void *update_cb_data)
1174 {
1175         int f, i=0, x, y, cancel=0;
1176         float progress;
1177         OceanResult ocr;
1178         ImBuf *ibuf_foam, *ibuf_disp, *ibuf_normal;
1179         float *prev_foam;
1180         int res_x = och->resolution_x;
1181         int res_y = och->resolution_y;
1182         char string[FILE_MAX];
1183         
1184         if (!o) return;
1185         
1186         prev_foam = MEM_callocN(res_x*res_y*sizeof(float), "previous frame foam bake data");
1187         
1188         BLI_srand(0);
1189         
1190         for (f=och->start, i=0; f<=och->end; f++, i++) {
1191                 
1192                 /* create a new imbuf to store image for this frame */
1193                 ibuf_foam = IMB_allocImBuf(res_x, res_y, 32, IB_rectfloat);
1194                 ibuf_disp = IMB_allocImBuf(res_x, res_y, 32, IB_rectfloat);
1195                 ibuf_normal = IMB_allocImBuf(res_x, res_y, 32, IB_rectfloat);
1196                 
1197                 ibuf_disp->profile = ibuf_foam->profile = ibuf_normal->profile = IB_PROFILE_LINEAR_RGB;
1198                 
1199                 BKE_simulate_ocean(o, och->time[i], och->wave_scale, och->chop_amount);
1200                 
1201                 /* add new foam */
1202                 for (y=0; y < res_y; y++) {
1203                         for (x=0; x < res_x; x++) {
1204                                 float r, pr=0.0, foam_result;
1205                                 float neg_disp, neg_eplus;
1206                                 
1207                                 BKE_ocean_eval_ij(o, &ocr, x, y);
1208                                                                 
1209                                 normalize_v3(ocr.normal);
1210                                 
1211                                 /* foam */
1212                                 ocr.foam = BKE_ocean_jminus_to_foam(ocr.Jminus, och->foam_coverage);
1213                                 
1214                                 /* accumulate previous value for this cell */
1215                                 if (i>0)
1216                                         pr = prev_foam[res_x*y + x];
1217
1218                                 r = BLI_frand();        // randomly reduce foam
1219                                 
1220                                 //pr = pr * och->foam_fade;             // overall fade
1221                                 
1222                                 // remember ocean coord sys is Y up!
1223                                 // break up the foam where height (Y) is low (wave valley), 
1224                                 // and X and Z displacement is greatest
1225                                 
1226                                 /*
1227                                  vec[0] = ocr.disp[0];
1228                                 vec[1] = ocr.disp[2];
1229                                 hor_stretch = len_v2(vec);
1230                                 CLAMP(hor_stretch, 0.0, 1.0);
1231                                 */
1232                                 
1233                                 neg_disp = ocr.disp[1]<0.0?1.0+ocr.disp[1]:1.0;
1234                                 neg_disp = neg_disp<0.0?0.0:neg_disp;
1235                                 
1236                                 neg_eplus = ocr.Eplus[2]<0.0?1.0+ocr.Eplus[2]:1.0;
1237                                 neg_eplus = neg_eplus<0.0?0.0:neg_eplus;
1238
1239                                 //if (ocr.disp[1] < 0.0 || r > och->foam_fade)
1240                                 //      pr *= och->foam_fade;
1241                                 
1242                                 
1243                                 //pr = pr * (1.0 - hor_stretch) * ocr.disp[1];
1244                                 //pr = pr * neg_disp * neg_eplus;
1245                                 
1246                                 if (pr < 1.0) pr *=pr;
1247                                 
1248                                 pr *= och->foam_fade * (0.75+neg_eplus*0.25);
1249                                 
1250                                 
1251                                 foam_result = pr + ocr.foam;
1252                                 
1253                                 prev_foam[res_x*y + x] = foam_result;
1254                                                                 
1255                                 /* add to the image */
1256                                 ibuf_disp->rect_float[4*(res_x*y + x) + 0] = ocr.disp[0];
1257                                 ibuf_disp->rect_float[4*(res_x*y + x) + 1] = ocr.disp[1];
1258                                 ibuf_disp->rect_float[4*(res_x*y + x) + 2] = ocr.disp[2];
1259                                 ibuf_disp->rect_float[4*(res_x*y + x) + 3] = 1.0;
1260                                 
1261                                 if (o->_do_jacobian) {
1262                                         ibuf_foam->rect_float[4*(res_x*y + x) + 0] = foam_result;
1263                                         ibuf_foam->rect_float[4*(res_x*y + x) + 1] = foam_result;
1264                                         ibuf_foam->rect_float[4*(res_x*y + x) + 2] = foam_result;
1265                                         ibuf_foam->rect_float[4*(res_x*y + x) + 3] = 1.0;
1266                                 }
1267                                 
1268                                 if (o->_do_normals) {
1269                                         ibuf_normal->rect_float[4*(res_x*y + x) + 0] = ocr.normal[0];
1270                                         ibuf_normal->rect_float[4*(res_x*y + x) + 1] = ocr.normal[1];
1271                                         ibuf_normal->rect_float[4*(res_x*y + x) + 2] = ocr.normal[2];
1272                                         ibuf_normal->rect_float[4*(res_x*y + x) + 3] = 1.0;
1273                                 }
1274
1275                         }
1276                 }
1277                 
1278                 /* write the images */
1279                 cache_filename(string, och->bakepath, f, CACHE_TYPE_DISPLACE);
1280                 if(0 == BKE_write_ibuf(ibuf_disp, string, R_OPENEXR, R_OPENEXR_HALF, 2))  // 2 == ZIP exr codec
1281                         printf("Cannot save Displacement File Output to %s\n", string);
1282                 
1283                 if (o->_do_jacobian) {
1284                         cache_filename(string, och->bakepath, f, CACHE_TYPE_FOAM);
1285                         if(0 == BKE_write_ibuf(ibuf_foam, string, R_OPENEXR, R_OPENEXR_HALF, 2))  // 2 == ZIP exr codec
1286                                 printf("Cannot save Foam File Output to %s\n", string);
1287                 }
1288                 
1289                 if (o->_do_normals) {
1290                         cache_filename(string, och->bakepath, f, CACHE_TYPE_NORMAL);
1291                         if(0 == BKE_write_ibuf(ibuf_normal, string, R_OPENEXR, R_OPENEXR_HALF, 2))  // 2 == ZIP exr codec
1292                                 printf("Cannot save Normal File Output to %s\n", string);
1293                 }
1294                 
1295                 IMB_freeImBuf(ibuf_disp);
1296                 IMB_freeImBuf(ibuf_foam);
1297                 IMB_freeImBuf(ibuf_normal);
1298                 
1299                 progress = (f - och->start) / (float)och->duration;
1300                 
1301                 update_cb(update_cb_data, progress, &cancel);
1302                 
1303                 if (cancel) {
1304                         MEM_freeN(prev_foam);
1305                         return;
1306                 }
1307         }
1308         
1309         MEM_freeN(prev_foam);
1310         och->baked = 1;
1311 }
1312
1313 #else // WITH_OCEANSIM
1314
1315 /* stub */
1316 typedef struct Ocean {
1317 } Ocean;
1318
1319
1320 float BKE_ocean_jminus_to_foam(float UNUSED(jminus), float UNUSED(coverage)) {
1321         return 0.0f;
1322 }
1323
1324 void BKE_ocean_eval_uv(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(u),float UNUSED(v))
1325 {
1326 }
1327
1328 // use catmullrom interpolation rather than linear
1329 void BKE_ocean_eval_uv_catrom(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(u),float UNUSED(v))
1330 {
1331 }
1332
1333 void BKE_ocean_eval_xz(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(x),float UNUSED(z))
1334 {
1335 }
1336
1337 void BKE_ocean_eval_xz_catrom(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(x),float UNUSED(z))
1338 {
1339 }
1340
1341 void BKE_ocean_eval_ij(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), int UNUSED(i),int UNUSED(j))
1342 {
1343 }
1344
1345 void BKE_simulate_ocean(struct Ocean *UNUSED(o), float UNUSED(t), float UNUSED(scale), float UNUSED(chop_amount))
1346 {
1347 }
1348
1349 struct Ocean *BKE_add_ocean(void)
1350 {
1351         Ocean *oc = MEM_callocN(sizeof(Ocean), "ocean sim data");
1352
1353         return oc;
1354 }
1355
1356 void BKE_init_ocean(struct Ocean* UNUSED(o), int UNUSED(M),int UNUSED(N), float UNUSED(Lx), float UNUSED(Lz), float UNUSED(V), float UNUSED(l), float UNUSED(A), float UNUSED(w), float UNUSED(damp), 
1357                                            float UNUSED(alignment), float UNUSED(depth), float UNUSED(time), short UNUSED(do_height_field), short UNUSED(do_chop), short UNUSED(do_normals), short UNUSED(do_jacobian), int UNUSED(seed))
1358 {
1359 }
1360
1361 void BKE_free_ocean_data(struct Ocean *UNUSED(oc)) 
1362 {
1363 }
1364
1365 void BKE_free_ocean(struct Ocean *oc) 
1366 {
1367         if(!oc) return;
1368         MEM_freeN(oc);
1369 }
1370
1371
1372 /* ********* Baking/Caching ********* */
1373
1374
1375 void BKE_free_ocean_cache(struct OceanCache *och)
1376 {
1377         if (!och) return;
1378         
1379         MEM_freeN(och);
1380 }
1381
1382 void BKE_ocean_cache_eval_uv(struct OceanCache *UNUSED(och), struct OceanResult *UNUSED(ocr), int UNUSED(f), float UNUSED(u), float UNUSED(v))
1383 {
1384 }
1385
1386 void BKE_ocean_cache_eval_ij(struct OceanCache *UNUSED(och), struct OceanResult *UNUSED(ocr), int UNUSED(f), int UNUSED(i), int UNUSED(j))
1387 {       
1388 }
1389
1390 struct OceanCache *BKE_init_ocean_cache(char *UNUSED(bakepath), int UNUSED(start), int UNUSED(end), float UNUSED(wave_scale), 
1391                                                   float UNUSED(chop_amount), float UNUSED(foam_coverage), float UNUSED(foam_fade), int UNUSED(resolution))
1392 {
1393         OceanCache *och = MEM_callocN(sizeof(OceanCache), "ocean cache data");
1394         
1395         return och;
1396 }
1397
1398 void BKE_simulate_ocean_cache(struct OceanCache *UNUSED(och), int UNUSED(frame))
1399 {
1400 }
1401
1402 void BKE_bake_ocean(struct Ocean *UNUSED(o), struct OceanCache *UNUSED(och), void (*update_cb)(void *, float progress, int *cancel), void *UNUSED(update_cb_data))
1403 {
1404         /* unused */
1405         (void)update_cb;
1406 }
1407 #endif // WITH_OCEANSIM