R = 40.
H = 60
x_w = 30.
y_w = 90.
Rw = 20
seed = 1

def f(x,y):
    np.random.seed(seed)
    return - 3.*np.exp( -((x-x_w)**2. + (y-y_w)**2.)/(Rw**2.) ) + \
           - 3.*np.exp( -((x-30)**2. + (y-40)**2.)/(20**2.) ) + \
           - 1.5*np.exp( -((x-np.random.uniform(-60,60))**2. + (y-np.random.uniform(40,70))**2.)/(60.**2.) ) #+ 4.*(y/H)**0.0 
    

"""
Explanation: ```python
R = 40.
H = 60
x_w = 0.
y_w = 60.
f = lambda x,y: - 4.np.exp( -((x-x_w)2. + (y-y_w)2.)/(30.2.) ) + \
                - 4.np.exp( -((x+70.)2. + (y-60.)2.)/(60.2.) ) #+ 15.*(y/H)0.05 
```
End of explanation
"""


x_MC = np.random.uniform(-R,R,9e6)
y_MC = np.random.uniform(-R,R,9e6) + H

r = np.sqrt(x_MC**2. + (y_MC-H)**2.)
ind_in = np.where(r <= R)

u_eq_MC = np.mean(f(x_MC[ind_in],y_MC[ind_in]))
u_eq_MC

x_mesh,y_mesh = np.meshgrid(np.arange(-50.,100.,0.1),np.arange(0.,120.,0.1))

F = f(x_mesh,y_mesh)

plt.contourf(x_mesh,y_mesh,F,20)
cb = plt.colorbar()
circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)
plt.plot(x_MC[ind_in],y_MC[ind_in],'.k',alpha=0.005)
x_plot = np.arange(-50.,50.1,0.1)

"""
Explanation: Monte-Carlo integration
End of explanation
"""


N_x = 5
N_y = 5
root_x_1D, weight_x = np.polynomial.legendre.leggauss(N_x)
root_y_1D, weight_y = np.polynomial.legendre.leggauss(N_y)

root_x,root_y = np.meshgrid(root_x_1D,root_y_1D)
root_x,root_y = root_x.flatten(), root_y.flatten()

print root_x
print root_y
print
print root_x_1D
print root_y_1D
print weight_x
print weight_y

y1 = lambda x: np.sqrt(R**2. - x**2.) + H
y2 = lambda x: -np.sqrt(R**2. - x**2.) + H

xi_1D = R*root_x_1D
xi = R*root_x
yij = 0.5*(y1(xi) - y2(xi))*root_y + 0.5*(y1(xi) + y2(xi))
fij = f(xi,yij)

fij.reshape(N_x,N_y)

np.dot(fij.reshape(N_x,N_y),weight_y)

G_xi = 0.5*(y1(xi_1D) - y2(xi_1D)) * np.dot(fij.reshape(N_x,N_y),weight_y)

I = np.dot(R*G_xi,weight_x)/(np.pi*R**2.)
print I
print 100*(I-u_eq_MC)/u_eq_MC, '%'

plt.contourf(x_mesh,y_mesh,F,20)
plt.colorbar()

circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)

x_plot = np.arange(-50.,50.2,0.1)
plt.plot(x_plot,y1(x_plot))
plt.plot(x_plot,y2(x_plot))
plt.plot(xi,yij,'ok')



"""
Explanation: Cartessian Gausss-Legendre quadrature
End of explanation
"""


N_x = 5
N_y = 5
root_x_1D, weight_x = np.polynomial.legendre.leggauss(N_x)
root_y_1D, weight_y = np.polynomial.legendre.leggauss(N_y)

root_x,root_y = np.meshgrid(root_x_1D,root_y_1D)
root_x,root_y = root_x.flatten(), root_y.flatten()

x1 = lambda y: np.sqrt(R**2. - (y-H)**2.)
x2 = lambda y: -np.sqrt(R**2. - (y-H)**2.)

yi_1D = R*root_y_1D+H
yi = R*root_y+H
xij = 0.5*(x1(yi) - x2(yi))*root_x + 0.5*(x1(yi) + x2(yi))
fij = f(xij,yi)

fij.reshape(N_x,N_y)

G_yi = 0.5*(x1(yi_1D) - x2(yi_1D)) * np.dot(weight_x,fij.reshape(N_x,N_y))

I = np.dot(R*G_yi,weight_y)/(np.pi*R**2.)
print I
print 100*(I-u_eq_MC)/u_eq_MC, '%'

plt.contourf(x_mesh,y_mesh,F,20)
plt.colorbar()

circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)

y_plot = np.arange(-50.,50.2,0.1)+H
plt.plot(x1(y_plot),y_plot)
plt.plot(x2(y_plot),y_plot)
plt.plot(xij,yi,'ok')


"""
Explanation: integrating before in y
End of explanation
"""


f_r = lambda r,th: f(-r*np.sin(th),r*np.cos(th)+H)

def gaussN(R, func, varargin, NGr=4, NGth=4):
    """Calculate numerically the gauss integration.
    [1] eq. 38

    Inputs
    ----------
    R (float): Wind turbine radius [m]
    func (function): Wind speed function
    varargin: Other arguments for the function besides [r,te]
    NG (int): Number of Ga

    Outputs
    ----------
    Ua (float):
    """
    A = np.pi*R**2
    #coefficients
    if  (NGr==4)&(NGth==4): #for speed give the full values
        rt = np.array([[ -0.339981043584856, -0.861136311594053,
            0.339981043584856, 0.861136311594053]])
        te = rt.T
        w = np.array([[0.652145154862546, 0.347854845137454,
            0.652145154862546, 0.347854845137454]])
        wt=w
    else:
        rt,w = np.polynomial.legendre.leggauss(NGr)
        rt = np.array([rt])
        #te = rt.T
        w = np.array([w])
        te,wt = np.polynomial.legendre.leggauss(NGr)
        te = np.array([te]).T
        wt = np.array([wt])
        
    return np.sum((np.pi/4.0)*(R**2./A)*w*wt.T*func(R*(rt+1.0)/2.0,
        np.pi*(te+1.0),*varargin)*(rt+1.0))

N_r = 4
N_th = 4
I = gaussN(R,f_r,[],N_r,N_th)
print I
print 100*(I-u_eq_MC)/u_eq_MC, '%'

plt.contourf(x_mesh,y_mesh,F,20)
plt.colorbar()

circ = plt.Circle([0.,H],R,ec='k',fill=False)
ax = plt.gca()
ax.add_patch(circ)


rt,w = np.polynomial.legendre.leggauss(N_r)
rt = np.array([rt])
te,wt = np.polynomial.legendre.leggauss(N_th)
te = np.array([te])
re = R*(rt+1.0)/2.0
te = np.pi*(te.T+1.0)
plt.plot(-re*np.sin(te),re*np.cos(te)+H,'ok')



"""
Explanation: Radial Gauss-Leg
End of explanation
"""


