#! usr/bin/python3

# this files plots the analytical solution for acoustic wave generation by 2 spinning vortices.
# the domane normaly is 2D square can be seperated into inner region and outer region.
# the inner region is described by incompressible fluid motion. 
# The analytical solution is given below.
# The outer region is the acoustic field.
# The acoustic pressure analytical solution is also given below.
# For more information see:
# Manuel Hasert, Joerg Bernsdorf and Sabine Roller, "Towards aeroacoustic sound generation by flow through porous media"
# Dipl.–Math. techn. Roland Fortenbach, "Mehrskalenmodellierung von aeroakustischen Quellen in schwach kompressiblen Stroemungen"


# Raphael Haupt
# 20.Sept.2016

import sys
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hankel2
from subprocess import call

# generating var.py file to import variables from seeder and musubi lua files.
# print.lua print out the parameters of musubi and seeder lua files.
call ("lua print.lua > var.py", shell=True)

from var import rho0 as rho0
from var import circulation as gamma
from var import r0 as r0
from var import Cs as Cs
from var import Length as length
from var import Level as level
from var import Ma as Ma

r= np.linspace(0,(2*length), (2**level))  # radius
x= np.linspace(-20,20, (2**level))        # x-values
y=0                                       # y position
theta=0                                   # angle
t=0                                       # time
p0=0                                      # background pressure
ommega=gamma/(4*math.pi*r0**2)            # angular velocity 
k=2*ommega/Cs                             
z=x+1j*y
b=r0*np.exp(1j*ommega*t)
ommegaz = gamma/(math.pi*1j)*z/(z**2-b**2)

# velocity analytical solution
u=np.real(ommegaz)                        # vel u
v=-1*np.imag(ommegaz)                     # vel v

# accustic pressure amplitude analytical solution
p_amp = np.real((-1j*Ma**4*(rho0*gamma**4)/( 64 * math.pi**3 * r0**4 * Cs**2 ))*hankel2(k,r)*np.exp(2*1j*(ommega*t-theta))) 

# fluid pressure analytical solution
PHIt = np.real( (-gamma*ommega/math.pi)*(b**2/(z**2-b**2)))
p = p0 - Ma**2*rho0*(PHIt + u**2 + v**2 / 2)

plt.figure(1)
plt.subplot(211)
plt.plot(r,p_amp,label='p_amp')

plt.subplot(212)
plt.plot(x,u, x,v,x,p,)

plt.show()


