N = 100                               # Number of points in each direction
x_start, x_end = -5, 7            # x-direction boundaries
y_start, y_end = -7, 7            # y-direction boundaries
x = numpy.linspace(x_start, x_end, N)    # computes a 1D-array for x
y = numpy.linspace(y_start, y_end, N)    # computes a 1D-array for y
X, Y = numpy.meshgrid(x, y)              # generates a mesh grid

"""
Explanation: Line segment source
Generate meshgrid for pressure and velocity field
End of explanation
"""


#coordinate transformation
origin_new_x=0            # A point in Fig 1
origin_new_y=0            

X_local=X-origin_new_x
Y_local=Y-origin_new_y

"""
Explanation: Line segment source integration
Line source solution for pressure and velocity (Datta-Gupta, 2007)
\begin{equation}
P(x,y)=B{{Q}{w}}=-\frac{70.60\mu }{h\sqrt{{{k}{x}}{{k}{y}}}}\ln \left{ {{(x-{{x}{w}})}^{2}}+\frac{{{k}{x}}}{{{k}{y}}}{{(y-{{y}{w}})}^{2}} \right}{{Q}{w}}+{{P}_{avg}}
\end{equation}
\begin{equation}
\frac{\partial P}{\partial x}=u=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}{x}}}{{{k}{y}}}}\sum\limits_{k=1}^{{{N}{w}}}{{{Q}{k}}}\frac{x-{{x}{k}}}{{{\left( x-{{x}{k}} \right)}^{2}}+\frac{{{k}{x}}}{{{k}{y}}}{{(y-{{y}_{k}})}^{2}}}
\end{equation}
\begin{equation}
\frac{\partial P}{\partial y}=v=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}{x}}}{{{k}{y}}}}\sum\limits_{k=1}^{{{N}{w}}}{{{Q}{k}}}\frac{y-{{y}{k}}}{{{\left( x-{{x}{k}} \right)}^{2}}+\frac{{{k}{x}}}{{{k}{y}}}{{(y-{{y}_{k}})}^{2}}}
\end{equation}
Local coordinate transformation
<img src="./resources/PanelLocal.png" width="300">
<center>Figure 1. Nomenclature of the boundary element in the local coordinates</center>
Based on the nomenclatures in Fig. 1, a boundary element could be determined as its start point (A), end point (B), length (lj) and inclination (αj). The coordinate transformation from global system to local system can be given as follows:
\begin{matrix}
   {x}'=x-{{x}{A}} & {y}'=y-{{y}{A}}  \
\end{matrix}
Where (xA,yA) is the location of start point (point A in Fig. A-1) of a BE in global coordinate system. Trigonometric sine and cosine function for a BE could be determined from start point (A) and end point (B) of the boundary elements, which can be expressed as follows:
\begin{matrix}
   \sin {{\alpha }{j}}=\frac{{{y}{B}}-{{y}{A}}}{{{l}{j}}} & \cos {{\alpha }{j}}=\frac{{{x}{B}}-{{x}{A}}}{{{l}{j}}}  \
\end{matrix}
In local coordinate system, line source are placed along the boundary element at a point (x’j,y’j). The location can be easily given as follows:
\begin{equation}
\begin{matrix}
   \begin{matrix}
   {{{{x}'}}{j}}=t\cos {{\alpha }{j}} & {{{{y}'}}{j}}=t\sin {{\alpha }{j}}  \
\end{matrix} & 0\le t  \
\end{matrix}\le {{l}_{j}}
\end{equation}
End of explanation
"""


#More General method
Q = 2000    #strength of the source-sheet,stb/d
h=25.26        #thickness of local gridblock,ft
phi=0.2        #porosity 
kx=200         #pemerability in x direction,md
ky=200         #pemerability in y direction,md
kr=kx/ky       #pemerability ratio
miu=1          #viscosity,cp

#Line segment source parameters
L1=1                      #Line segment source length
alpha=30*numpy.pi/180     #Line segment angle from x direction, angle alpha in Fig.1

# computes the pressure-velcity field generated by the source-sheet
p_sheet1 = numpy.empty((N, N), dtype=float)
u_sheet1 = numpy.empty((N, N), dtype=float)
v_sheet1 = numpy.empty((N, N), dtype=float)
for i in range(N):
    for j in range(N):  
        a=math.cos(alpha)**2+kr*math.sin(alpha)**2
        b=X_local[i,j]*math.cos(alpha)+kr*math.sin(alpha)*Y_local[i,j]
        c=Y_local[i,j]*math.cos(alpha)-X_local[i,j]*math.sin(alpha)
        dp=70.6*Q*miu/h/math.sqrt(kx*ky)
        dv=-0.4468*Q/h/phi*math.sqrt(kx/ky)
        
        p_sheet1[i,j] = dp/a*(
                             (
                              b*math.log(X_local[i,j]**2-2*b*L1+a*L1**2+kr*Y_local[i,j]**2)
                             -L1*a*math.log((X_local[i,j]-L1*math.cos(alpha))**2+kr*(Y_local[i,j]-L1*math.sin(alpha))**2)
                             +2*math.sqrt(kr)*c*math.atan((b-a*L1)/math.sqrt(kr)/c)
                             )
                             -
                             (
                               b*math.log(X_local[i,j]**2+kr*Y_local[i,j]**2)
                               +2*math.sqrt(kr)*c*math.atan((b)/math.sqrt(kr)/c)
                             )         
                            )
        u_sheet1[i,j] = dv/a*(
                             ( 
                            math.cos(alpha)*math.log(X_local[i,j]**2-2*b*L1+a*L1**2+kr*Y_local[i,j]**2)+ 2*math.sqrt(kr)*math.sin(alpha)*math.atan((a*L1-b)/math.sqrt(kr)/c) 
                             )
                            -
                             (
                            math.cos(alpha)*math.log(X_local[i,j]**2+kr*Y_local[i,j]**2)+2*math.sqrt(kr)*math.sin(alpha)*math.atan((-b)/math.sqrt(kr)/c)
                             )    
                            )
        v_sheet1[i,j] = dv/a*(
                             ( 
                            math.sin(alpha)*math.log(X_local[i,j]**2-2*b*L1+a*L1**2+kr*Y_local[i,j]**2)+ 2*math.sqrt(1/kr)*math.cos(alpha)*math.atan((b-a*L1)/math.sqrt(kr)/c) 
                             )
                            -
                             (
                            math.sin(alpha)*math.log(X_local[i,j]**2+kr*Y_local[i,j]**2)+2*math.sqrt(1/kr)*math.cos(alpha)*math.atan((b)/math.sqrt(kr)/c)
                             )    
                            )

# superposition of the source-sheet to the uniform flow
p1 = p_sheet1+ p_sheet1.min()*-1 #make it looks better
#p1 = p_sheet1 #make it looks better
u1 = u_sheet1
v1 = v_sheet1

"""
Explanation: Line segment source integration
Integrating the line source solution (Eqs. 3-4) from A to B over a length l, the velocity can be given as follows:
\begin{equation}
P({x}',{y}')=-\frac{70.60\mu }{h\sqrt{{{k}{x}}{{k}{y}}}}{{Q}{w}}\int{t=0}^{t={{l}{j}}}{\ln \left{ {{({x}'-t\cos {{\alpha }{j}})}^{2}}+\frac{{{k}{x}}}{{{k}{y}}}{{({y}'-t\sin {{\alpha }_{j}})}^{2}} \right}dt}
\end{equation}
\begin{equation}
u\left( {x}',{y}' \right)=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}{x}}}{{{k}{y}}}}{{Q}{j}}\int{t=0}^{t={{l}{j}}}{\frac{{x}'-t\cos {{\alpha }{j}}}{{{\left( {x}'-t\cos {{\alpha }{j}} \right)}^{2}}+\frac{{{k}{x}}}{{{k}{y}}}{{({y}'-t\sin {{\alpha }{j}})}^{2}}}dt}
\end{equation}
\begin{equation}
v\left( {x}',{y}' \right)=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}{x}}}{{{k}{y}}}}{{Q}{j}}\int{t=0}^{t={{l}{j}}}{\frac{{y}'-t\sin {{\alpha }{j}}}{{{\left( {x}'-t\cos {{\alpha }{j}} \right)}^{2}}+\frac{{{k}{x}}}{{{k}{y}}}{{({y}'-t\sin {{\alpha }{j}})}^{2}}}dt}
\end{equation}
The integrals can be rewritten as follows:
\begin{equation}
\left. P({x}',{y}')={{Q}{w}}\frac{{{d}{p}}}{a}\left[ b\ln ({{{{x}'}}^{2}}-2bt+a{{t}^{2}}+{{k}{r}}{{{{y}'}}^{2}})-ta\ln \left{ {{({x}'-t\cos {{\alpha }{j}})}^{2}}+{{k}{r}}{{({y}'-t\sin {{\alpha }{t}})}^{2}} \right}+2\sqrt{{{k}{r}}}c\arctan \left( \frac{b-at}{\sqrt{{{k}{r}}}c} \right) \right] \right|{t=0}^{t={{l}{j}}}
\end{equation}
\begin{equation}
\left. u\left( {x}',{y}' \right)={{Q}{j}}\frac{{{d}{v}}}{a}\left[ \cos {{\alpha }{j}}\ln ({{{{x}'}}^{2}}-2bt+a{{t}^{2}}+{{k}{r}}{{{{y}'}}^{2}})+2\sqrt{{{k}{r}}}\sin {{\alpha }{j}}\arctan (\frac{at-b}{\sqrt{{{k}{r}}}c}) \right] \right|{t=0}^{t={{l}_{j}}}
\end{equation}
\begin{equation}
v\left( {x}',{y}' \right)=\left. {{Q}{j}}\frac{{{d}{v}}}{a}\left[ \sin {{\alpha }{j}}\ln ({{{{x}'}}^{2}}-2tb+a{{t}^{2}}+{{k}{r}}{{{{y}'}}^{2}})+2\sqrt{\frac{1}{{{k}{r}}}}\cos {{\alpha }{j}}\arctan (\frac{b-at}{\sqrt{{{k}{r}}}c}) \right] \right|{t=0}^{t={{l}_{j}}}
\end{equation}
where,
\begin{equation}
a={{\cos }^{2}}{{\alpha }{j}}+{{k}{r}}{{\sin }^{2}}{{\alpha }{j}}\end{equation}
\begin{equation}
b={x}'\cos {{\alpha }{j}}+{y}'{{k}{r}}\sin {{\alpha }{j}}\end{equation}
\begin{equation}
c={y}'\cos {{\alpha }{j}}-{x}'\sin {{\alpha }{j}}\end{equation}
\begin{equation}
{{d}{v}}=-\frac{0.4468}{h\phi }\sqrt{{{k}{r}}}\begin{matrix}
   {} & {{d}{p}}=\frac{70.60\mu }{h\sqrt{{{k}{x}}{{k}{y}}}}  \
\end{matrix}\end{equation}
\begin{equation}
{{k}{r}}={{k}{x}}/{{k}{y}}\end{equation}
Code implementation
End of explanation
"""


# plots the streamlines
%matplotlib inline

size = 8
pyplot.figure(figsize=(size, size-1))
pyplot.title('Streamline field')
pyplot.grid(True)
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.plot((0+origin_new_x,L1*math.cos(alpha)+origin_new_x),(0+origin_new_y,L1*math.sin(alpha)+origin_new_y), color='#CD2305', linewidth=4) #red source line

#stream =pyplot.streamplot(X, Y, u, v, color=numpy.sqrt(u**2+v**2),cmap='jet',density=2, linewidth=1, arrowsize=1, arrowstyle='->') #streamline
stream =pyplot.streamplot(X, Y, u1, v1, color=numpy.sqrt(u1**2+v1**2),cmap='jet',density=2, linewidth=1, arrowsize=1, arrowstyle='->') #streamline

cbar=pyplot.colorbar(orientation='vertical')
equipotential=pyplot.contour(X, Y, p1, extend='both')

size = 8
pyplot.figure(figsize=(size, size-1))
pyplot.title('Velocity field')
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)

Vtotal= numpy.sqrt(u1**2+v1**2)

pyplot.contour(X, Y, Vtotal, 10, linewidths=0.5, colors='k')
pyplot.contourf(X, Y, Vtotal, 10, cmap='rainbow')
pyplot.colorbar()  # draw colorbar

"""
Explanation: Plot the results
End of explanation
"""


