# -*- 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, 2, 55, 18), 'numpy.mean', 'np.mean', ({(55, 10, 55, 17): 'K2[25:]'}, {}), '(K2[25:])', True, 'import numpy as np\n'), ((57, 0, 57, 59), 'matplotlib.pyplot.plot', 'plt.plot', (), '', True, 'import matplotlib.pyplot as plt\n'), ((59, 0, 59, 36), 'matplotlib.pyplot.xlabel', 'plt.xlabel', ({(59, 11, 59, 35): '"""Integration range $R$"""'}, {}), "('Integration range $R$')", True, 'import matplotlib.pyplot as plt\n'), ((60, 0, 60, 34), 'matplotlib.pyplot.ylabel', 'plt.ylabel', ({(60, 11, 60, 33): '"""Detection probabilty"""'}, {}), "('Detection probabilty')", True, 'import matplotlib.pyplot as plt\n'), ((61, 0, 61, 22), 'matplotlib.pyplot.legend', 'plt.legend', (), '', True, 'import matplotlib.pyplot as plt\n'), ((62, 0, 62, 29), 'matplotlib.pyplot.text', 'plt.text', ({(62, 9, 62, 12): '(2.4)', (62, 13, 62, 16): '(0.2)', (62, 17, 62, 28): '"""$k=1000$"""'}, {}), "(2.4, 0.2, '$k=1000$')", True, 'import matplotlib.pyplot as plt\n'), ((64, 0, 64, 36), 'matplotlib.pyplot.savefig', 'plt.savefig', (), '', True, 'import matplotlib.pyplot as plt\n'), ((65, 0, 65, 10), 'matplotlib.pyplot.show', 'plt.show', ({}, {}), '()', True, 'import matplotlib.pyplot as plt\n'), ((21, 7, 21, 20), 'numpy.tan', 'np.tan', ({(21, 14, 21, 19): 'alpha'}, {}), '(alpha)', True, 'import numpy as np\n'), ((23, 7, 23, 19), 'numpy.tan', 'np.tan', ({(23, 14, 23, 18): 'beta'}, {}), '(beta)', True, 'import numpy as np\n'), ((37, 9, 37, 43), 'numpy.linspace', 'np.linspace', (), '', True, 'import numpy as np\n'), ((38, 9, 38, 20), 'numpy.empty', 'np.empty', ({(38, 18, 38, 19): 'N'}, {}), '(N)', True, 'import numpy as np\n'), ((39, 9, 39, 20), 'numpy.empty', 'np.empty', ({(39, 18, 39, 19): 'N'}, {}), '(N)', True, 'import numpy as np\n'), ((18, 21, 18, 34), 'numpy.sin', 'np.sin', ({(18, 28, 18, 33): 'alpha'}, {}), '(alpha)', True, 'import numpy as np\n'), ((50, 11, 50, 35), 'numpy.sqrt', 'np.sqrt', ({(50, 19, 50, 34): '(real ** 2 + imag ** 2)'}, {}), '(real ** 2 + imag ** 2)', True, 'import numpy as np\n'), ((26, 14, 26, 38), 'numpy.sqrt', 'np.sqrt', ({(26, 22, 26, 37): '((xa - x) ** 2 + ya ** 2)'}, {}), '((xa - x) ** 2 + ya ** 2)', True, 'import numpy as np\n'), ((26, 42, 26, 66), 'numpy.sqrt', 'np.sqrt', ({(26, 50, 26, 65): '((xb - x) ** 2 + yb ** 2)'}, {}), '((xb - x) ** 2 + yb ** 2)', True, 'import numpy as np\n'), ((29, 15, 29, 31), 'numpy.sqrt', 'np.sqrt', ({(29, 23, 29, 30): '(xa ** 2 + 1)'}, {}), '(xa ** 2 + 1)', True, 'import numpy as np\n'), ((29, 45, 29, 61), 'numpy.sqrt', 'np.sqrt', ({(29, 53, 29, 60): '(xa ** 2 + 1)'}, {}), '(xa ** 2 + 1)', True, 'import numpy as np\n')]