In [23]:
import numpy as np
from math import *
import pandas as pd
from pylab import *

In [24]:
points = float(254)
modes = float(1)
seed = 5

In [25]:
def f(a, b, k1, k2, x, y):
    return a*cos(k1*x+k2*y) + b*sin(k1*x+k2*y)
def fx(a, b, k1, k2, x, y):
    return k1*(- a*sin(k1*x+k2*y) + b*cos(k1*x+k2&y) )
def fy(a, b, k1, k2, x, y):
    return k2*(- a*sin(k1*x+k2*y) + b*cos(k1*x+k2&y) )
def fxx(a, b, k1, k2, x, y):
    return k1*k1*(-a*cos(k1*x+k2*y) - b*sin(k1*x+k2&y))
def fyy(a, b, k1, k2, x, y):
    return k2*k2*(-a*cos(k1*x+k2*y) - b*sin(k1*x+k2&y))
def fxy(a, b, k1, k2, x, y):
    return k1*k2*(-a*cos(k1*x+k2*y) - b*sin(k1*x+k2&y))
def fxxx(a, b, k1, k2, x, y):
    return k1*k1*k1*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))
def fyyy(a, b, k1, k2, x, y):
    return k2*k2*k2*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))
def fxyy(a, b, k1, k2, x, y):
    return k1*k2*k2*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))
def fyxx(a, b, k1, k2, x, y):
    return k1*k1*k2*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))

Пример для трёх гармоник


In [26]:
mu, sigma = 0.0, 1.0 

a_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))
b_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))

In [27]:
def P(k1, k2, K0):
    if (k1 > K0) | (k2 > K0):
        return 0
    else:
        return 1

In [28]:
def func (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + ((fabs(k1-int(modes)))**2 - k2**2)*(a_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))

    return y
#((float(k1 - int(modes)))**2-(float(k2))**2)**P(k1, k2, int(modes))*

In [29]:
def func_x (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + (k1 - int(modes))*((fabs(k1-int(modes)))**2 - k2**2)*(-a_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))

    return y
#((float(k1 - int(modes)))**2-(float(k2))**2)**P(k1, k2, int(modes))*

In [30]:
def func_y (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + k2*((fabs(k1-int(modes)))**2 - k2**2)*(-a_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))

    return y
#((float(k1 - int(modes)))**2-(float(k2))**2)**P(k1, k2, int(modes))*

In [31]:
file = open("f1.dat", 'w')
for j1 in xrange(0, int(points)):
    for j2 in xrange (0, int(points)):
        x = j1/(points)
        y = j2/(points)
        file.write('{}    {}    {}\n'.format(y, x, func(j1, j2, a_coef, b_coef, int(modes), points)))
file.close()

In [32]:
file = open("derivatives1.dat", 'w')
for j1 in xrange(0, int(points)):
    for j2 in xrange (0, int(points)):
        x = j1/(points)
        y = j2/(points)
        file.write('{}    {}    {}    {}\n'.format(y, x, func_x(j1, j2, a_coef, b_coef, int(modes), points), func_y(j1, j2, a_coef, b_coef, int(modes), points)))
file.close()

In [33]:
def func (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + (k1-int(modes))*2.0*k2*(a_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))
    
    return y
#2*((float(k1 - int(modes)))*(float(k2)))*

In [34]:
def func_x (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + 2*(k1 - int(modes))*(k1-int(modes))*k2*(-a_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))

    return y
#((float(k1 - int(modes)))**2-(float(k2))**2)**P(k1, k2, int(modes))*

In [35]:
def func_y (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + 2*k2*(k1-int(modes))*k2*(-a_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))

    return y
#((float(k1 - int(modes)))**2-(float(k2))**2)**P(k1, k2, int(modes))*

In [36]:
mu, sigma = 0.0, 1.0 

a_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))
b_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))

In [37]:
file = open("f2.dat", 'w')
for j1 in xrange(0, int(points)):
    for j2 in xrange (0, int(points)):
        x = j1/(points)
        y = j2/(points)
        file.write('{}    {}    {}\n'.format(y, x, func(j1, j2, a_coef, b_coef, int(modes), points)))
file.close()

In [38]:
file = open("derivatives2.dat", 'w')
for j1 in xrange(0, int(points)):
    for j2 in xrange (0, int(points)):
        x = j1/(points)
        y = j2/(points)
        file.write('{}    {}    {}    {}\n'.format(y, x, func_x(j1, j2, a_coef, b_coef, int(modes), points), func_y(j1, j2, a_coef, b_coef, int(modes), points)))
file.close()

In [39]:
print a_coef


[[ 1.10471934 -1.21705879  0.45360344 ..., -1.76914744 -1.4691392
   1.15120101]
 [-0.23643842 -0.64194739 -0.74948277 ..., -0.30400607  1.11294846
  -0.79222186]
 [-0.31663419  0.28994381  1.71435204 ..., -0.54695454  0.69111755
  -2.02858315]
 ..., 
 [ 1.44574228  0.64694699 -1.20176958 ...,  1.72642793 -0.83883564
   2.70010344]
 [-0.66717089  0.50017932 -1.11163304 ...,  0.99294532 -0.88541582
   0.40333067]
 [ 0.82501874 -0.72485442  1.114924   ...,  0.99211124 -1.14004993
  -0.05524027]]

In [ ]:


In [ ]: