# -*- coding: utf-8 -*-
"""
Created on Fri Mar 15 16:51:16 2019

@author: Silvan
"""

import numpy as np
import scipy
import matplotlib.pyplot as plt

k=1000

n1=2.0
n2=1.0

alpha=np.pi/6.0
beta=np.arcsin(n2/n1*np.sin(alpha))

ya=1.0
xa=-ya*np.tan(alpha)
yb=-1.0
xb=-yb*np.tan(beta)

def s(x):
    return n1*np.sqrt((xa-x)**2+ya**2)+n2*np.sqrt((xb-x)**2+yb**2)

def kernel(xa,xb):
    return 1.0/np.sqrt(xa**2+1)**(3/2.0)+1.0/np.sqrt(xa**2+1)**(3/2.0)



def K(R):
    L=1000     #Maximum Number of subdivisions for integral calculations
    eps=0.01
    N=50
    x,dx=np.linspace(0.01,R,N,retstep=True)
    real=np.empty(N)
    imag=np.empty(N)
    real[0]=scipy.integrate.quad(lambda x: np.cos(k*s(x)),-x[0],x[0],epsrel=eps,limit=L)[0]
    imag[0]=scipy.integrate.quad(lambda x: np.sin(k*s(x)),-x[0],x[0],epsrel=eps,limit=L)[0]
    for i in range(1,N):
        r1=scipy.integrate.quad(lambda x: np.cos(k*s(x)),-x[i]-dx,-x[i],epsrel=eps,limit=L)[0]
        r2=scipy.integrate.quad(lambda x: np.cos(k*s(x)),x[i],x[i]+dx,epsrel=eps,limit=L)[0]
        real[i]=real[i-1]+r1+r2
        i1=scipy.integrate.quad(lambda x: np.sin(k*s(x)),-x[i]-dx,-x[i],epsrel=eps,limit=L)[0]
        i2=scipy.integrate.quad(lambda x: np.sin(k*s(x)),x[i],x[i]+dx,epsrel=eps,limit=L)[0]
        imag[i]=imag[i-1]+i1+i2
    
    return np.sqrt(real**2+imag**2),x,real,imag



K2,x,r,i=K(3)
M=np.mean(K2[25:])

plt.plot(x,K2/M,label=r'$|\int_{-R}^{R}e^{i k s(x)}dx|^2$')
#plt.errorbar(x,K2/M,0.1*K2/M)
plt.xlabel(r'Integration range $R$')
plt.ylabel('Detection probabilty')
plt.legend(loc='best')
plt.text(2.4,0.2,r'$k=1000$')
#plt.text(1.1,0.5,r'$|\int_{-R}^{R}e^{i k s(x)}dx|^2$',fontsize=20)
plt.savefig('refraction_v3',dpi=200)
plt.show()

#N=20
#
#dx=np.linspace(0,10,N)
#
#P=np.ones(N)
#
#for i in range(N):
#    print(i+1)
#    P[i]=trans_amp(dx[i])
#
#
#plt.figure(1)
#plt.plot(dx,P/np.mean(P[20:]))
#plt.text(4.0,0.5,r'$|\int_{-\Delta x}^{\Delta x} e^{ik s(x)}dx$|',fontsize=20)
#plt.ylabel('Transition Amplitude')
#plt.xlabel(r'Integration Interval $ \Delta x$')
##plt.axis([0,10,0,1.1])
#plt.legend(loc='best')
##plt.savefig('refraction',dpi=200)
#plt.show()

#x=np.linspace(-5,5,100)
#
#plt.figure(2)
#plt.plot(x,s(x))
#plt.show()
#
#d=np.linspace(0,5,100)
#xa=-d/2
#xb=d/2
#plt.figure(3)
#plt.plot(d,kernel(xa,xb)**2)
#plt.show()
[(55, 'numpy.mean', 'np.mean', (['K2[25:]'], {}), True, 'import numpy as np\n'), (57, 'matplotlib.pyplot.plot', 'plt.plot', (['x', '(K2 / M)'], {'label': '"""$|\\\\int_{-R}^{R}e^{i k s(x)}dx|^2$"""'}), True, 'import matplotlib.pyplot as plt\n'), (59, 'matplotlib.pyplot.xlabel', 'plt.xlabel', (['"""Integration range $R$"""'], {}), True, 'import matplotlib.pyplot as plt\n'), (60, 'matplotlib.pyplot.ylabel', 'plt.ylabel', (['"""Detection probabilty"""'], {}), True, 'import matplotlib.pyplot as plt\n'), (61, 'matplotlib.pyplot.legend', 'plt.legend', ([], {'loc': '"""best"""'}), True, 'import matplotlib.pyplot as plt\n'), (62, 'matplotlib.pyplot.text', 'plt.text', (['(2.4)', '(0.2)', '"""$k=1000$"""'], {}), True, 'import matplotlib.pyplot as plt\n'), (64, 'matplotlib.pyplot.savefig', 'plt.savefig', (['"""refraction_v3"""'], {'dpi': '(200)'}), True, 'import matplotlib.pyplot as plt\n'), (65, 'matplotlib.pyplot.show', 'plt.show', ([], {}), True, 'import matplotlib.pyplot as plt\n'), (21, 'numpy.tan', 'np.tan', (['alpha'], {}), True, 'import numpy as np\n'), (23, 'numpy.tan', 'np.tan', (['beta'], {}), True, 'import numpy as np\n'), (37, 'numpy.linspace', 'np.linspace', (['(0.01)', 'R', 'N'], {'retstep': '(True)'}), True, 'import numpy as np\n'), (38, 'numpy.empty', 'np.empty', (['N'], {}), True, 'import numpy as np\n'), (39, 'numpy.empty', 'np.empty', (['N'], {}), True, 'import numpy as np\n'), (18, 'numpy.sin', 'np.sin', (['alpha'], {}), True, 'import numpy as np\n'), (50, 'numpy.sqrt', 'np.sqrt', (['(real ** 2 + imag ** 2)'], {}), True, 'import numpy as np\n'), (26, 'numpy.sqrt', 'np.sqrt', (['((xa - x) ** 2 + ya ** 2)'], {}), True, 'import numpy as np\n'), (26, 'numpy.sqrt', 'np.sqrt', (['((xb - x) ** 2 + yb ** 2)'], {}), True, 'import numpy as np\n'), (29, 'numpy.sqrt', 'np.sqrt', (['(xa ** 2 + 1)'], {}), True, 'import numpy as np\n'), (29, 'numpy.sqrt', 'np.sqrt', (['(xa ** 2 + 1)'], {}), True, 'import numpy as np\n')]