# set terminal pngcairo  transparent enhanced font "arial,10" fontscale 1.0 size 500, 350 
# set output 'bivariat.8.png'
set key inside right bottom vertical Right noreverse enhanced autotitles nobox
set samples 500, 500
set title "Finite summation of 10, 100, 1000 fourier coefficients" 
set xrange [ -10.0000 : 10.0000 ] noreverse nowriteback
set yrange [ -0.400000 : 1.20000 ] noreverse nowriteback
integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta))
int1a(x,d) = (x<=d*.1) ? 0 : (int1a(x-d,d)+(f(x-d)+4*f(x-d*.5)+f(x))*d/6.)
int1b(x,d) = (x>=-d*.1) ? 0 : (int1b(x+d,d)+(f(x+d)+4*f(x+d*.5)+f(x))*d/6.)
f(x)=sin(x-1)-.75*sin(2*x-1)+(x**2)/8-5
integral2_f(x,y) = (x<y)?int2(x,y,(y-x)/ceil((y-x)/delta)):                         -int2(y,x,(x-y)/ceil((x-y)/delta))
int2(x,y,d) = (x>y-d*.5) ? 0 : (int2(x+d,y,d) + (f(x)+4*f(x+d*.5)+f(x+d))*d/6.)
ack(m,n) = (m == 0) ? n + 1 : (n == 0) ? ack(m-1,1) : ack(m-1,ack(m,n-1))
min(x,y) = (x < y) ? x : y
max(x,y) = (x > y) ? x : y
gcd(x,y) = gcd1(rnd(max(x,y)),rnd(min(x,y)))
rnd(x) = int(x+0.5)
gcd1(x,y) = (y == 0) ? x : gcd1(y, x - x/y * y)
fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)
sum10(x)   = 1./2 + sum [k=1:10]   fourier(k, x)
sum100(x)  = 1./2 + sum [k=1:100]  fourier(k, x)
sum1000(x) = 1./2 + sum [k=1:1000] fourier(k, x)
delta = 0.2
GPFUN_integral_f = "integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta))"
GPFUN_int1a = "int1a(x,d) = (x<=d*.1) ? 0 : (int1a(x-d,d)+(f(x-d)+4*f(x-d*.5)+f(x))*d/6.)"
GPFUN_int1b = "int1b(x,d) = (x>=-d*.1) ? 0 : (int1b(x+d,d)+(f(x+d)+4*f(x+d*.5)+f(x))*d/6.)"
GPFUN_integral2_f = "integral2_f(x,y) = (x<y)?int2(x,y,(y-x)/ceil((y-x)/delta)):                         -int2(y,x,(x-y)/ceil((x-y)/delta))"
GPFUN_int2 = "int2(x,y,d) = (x>y-d*.5) ? 0 : (int2(x+d,y,d) + (f(x)+4*f(x+d*.5)+f(x+d))*d/6.)"
GPFUN_f = "f(x)=sin(x-1)-.75*sin(2*x-1)+(x**2)/8-5"
GPFUN_ack = "ack(m,n) = (m == 0) ? n + 1 : (n == 0) ? ack(m-1,1) : ack(m-1,ack(m,n-1))"
GPFUN_min = "min(x,y) = (x < y) ? x : y"
GPFUN_max = "max(x,y) = (x > y) ? x : y"
GPFUN_gcd = "gcd(x,y) = gcd1(rnd(max(x,y)),rnd(min(x,y)))"
GPFUN_gcd1 = "gcd1(x,y) = (y == 0) ? x : gcd1(y, x - x/y * y)"
GPFUN_rnd = "rnd(x) = int(x+0.5)"
GPFUN_fourier = "fourier(k, x) = sin(3./2*k)/k * 2./3*cos(k*x)"
k = 1000
GPFUN_sum10 = "sum10(x)   = 1./2 + sum [k=1:10]   fourier(k, x)"
GPFUN_sum100 = "sum100(x)  = 1./2 + sum [k=1:100]  fourier(k, x)"
GPFUN_sum1000 = "sum1000(x) = 1./2 + sum [k=1:1000] fourier(k, x)"
plot     sum10(x)   title "1./2 + sum [k=1:10]   sin(3./2*k)/k * 2./3*cos(k*x)",     sum100(x)  title "1./2 + sum [k=1:100]  sin(3./2*k)/k * 2./3*cos(k*x)",     sum1000(x) title "1./2 + sum [k=1:1000] sin(3./2*k)/k * 2./3*cos(k*x)"