Making a function periodic This might not be the best place to ask this question, but here it goes... I'm creating a game and need 3D sea waves. Since it's for mobiles, there's no time to generate entire screen worth of waves on every frame, so I figured I'll create a pattern and copy it over and over again. I created a function $h(x,y,time)$ which returns $z$ coordinate of point $T(x,y)$ where $x\in[0,m\rangle\cap\mathbb{N}$ and $y\in[0,n\rangle\cap\mathbb{N}$ based on modified simplex noise algorithm.
I had difficulty visualising algorithm for 3D waves, so I started with a simple graph plot, where y is the height of the wave. By running simplified function $h$ as $h(x)$, I get something like this:

Obviously, this is no good. You can clearly see that the pattern doesn't smoothly fade into next instance. Therefore, I calculated $y$ coordinate as
$$y=h(x)\cdot(1-\frac{x}{m})+h(x-m)\cdot\frac{x}{m}$$
This results in something like this:

It was good enough for me, so I attempted to expand it to 3D space. My math and logic seem non-existant from this point on, as the best I could come up with is this:
if (y-x < 0) then $A(x-y,0),B(m,n-x+y)$
else $A(0,y-x),B(m+x-y,n)$
$$f=\frac{\sqrt{(A_x-x)^2+(A_y-y)^2}}{\sqrt{(A_x-B_x)^2+(A_y-B_y)^2}}$$
$$z=h(x,y,0)\cdot(1-f)+h(x-m,y-n,0)\cdot f$$
This fails terribly. Here's the result of that expression - blue square is a $z=0$ plane. You can clearly see that the pattern doesn't blend smoothly.

I know that the expression should fail, but am unable to come up with a working expression. How can I make the function repeat itself smoothly like in the 2D example on both x and y axes?
Thanks in advance, Mirac7

EDIT: I've come up with the expression which should calculate height correctly.
$$z=\big(h(x,y,0)\cdot(1-\frac{x}{m})+h(x-m,y,0)\cdot\frac{x}{m}\big)\cdot(1-\frac{y}{n})+\big(h(x,y-n,0)\cdot(1-\frac{x}{m})+h(x-m,y-n,0)\cdot\frac{x}{m}\big)\cdot\frac{y}{n}$$
The result of this expression is this:

It's still not quite there. When I increase the $m$ and $n$, thereby gaining more precision from every point, I get smoother transition. I need to set $m,n$ to at least 128 points to get somewhat smooth transition, which is way too many to be calculated per frame.
But another problem occurred. It is barely noticeable in the above rendering, but here's what heightmap of some tiny waves with $m,n=128$ looks like:

Can someone please help me figure out why these edges are occurring? Thanks, Mirac7.
 A: If your $h(x,y)$ should be periodic with period $m$ in $x$-direction and periodic with period $n$ in $y$-direction choose a finite linear combination of terms of the form
$$(x,y)\mapsto \cos{2j\pi x\over m}\cos{2k\pi x\over m}$$
with $j$, $k\in{\mathbb N}$, and similar terms with $\cos(\>)\cos(\>)$ replaced by $\cos(\>)\sin(\>)$, etc.
Experiment with Mathematica or similar to obtain a wave plot of your liking.
A: Easiest is to do it in sequential manner - one axis at a time, then update array, and repeat for second axis. However, you can rather easily also derive explicit analytic formula (bi-liner polynomial) from that iteration. 
Here is simple python/numpy script (sorry, I put code instead for math formulas, but this is programming question anyway)
#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt

# ========= functions

def blob(p,w,Xs,Ys):
    # helps generate nice test-function
    dXs = Xs-p[0]
    dYs = Ys-p[1]
    return 1/( 1 + (dXs**2+dYs**2)/(w*w) )

def matchPBC(f,Xs,Ys):
   '''
      match periodic boundary condions iteratevily (one axis at a tim) 
   '''
    f=f.copy()
    ### X-boundary
    Dxs = f[:,0] - f[:,-1];    # difference for each pixel at boundary
    f += Dxs[:,None]*(Xs-0.5)  # linear interpolation for each x-line[!
    ### Y-boundary
    Dys = f[0,:] - f[-1,:]; 
    f += Dys[None,:]*(Ys-0.5) 
    return f

# ========= main

### perpare support arryas
xs = np.linspace(0.0,1.0,100,endpoint=False)
ys = np.linspace(0.0,1.0,100,endpoint=False)
Xs,Ys = np.meshgrid(xs,ys)

#### test function
#f = np.sin(Xs*6) + np.sin(Xs*11) + np.cos(Xs*7-Ys*9) + np.sin(Xs*13+Ys*14) #+ np.sin(Xs*33+Ys*24)
#f = np.cos(Xs*3.5)  + np.sin(Ys*8)
#f = (Xs*5 + Ys*16)/21.0
f = blob([0.0,0.9],0.2,Xs,Ys) + blob([0.5,0.0],0.2,Xs,Ys) 

#### do it !!!
fpbc = matchPBC(f,Xs,Ys)

#### plot results !!!
plt.subplot(1,2,1); plt.imshow( np.tile(f   ,(2,2)) )
plt.subplot(1,2,2); plt.imshow( np.tile(fpbc,(2,2)) )
plt.show()

Results looks like this
Left: input non-periaodic function; Right: periodized resulting function ]1
