Last active
April 21, 2018 09:19
-
-
Save matejak/c53ae7adb5ff016b8d9b6b3f593d3333 to your computer and use it in GitHub Desktop.
Revisions
-
matejak revised this gist
Apr 21, 2018 . No changes.There are no files selected for viewing
-
Apr 21, 2018 . 1 changed file with 73 additions and 103 deletions.There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -1,67 +1,15 @@ # -*- coding: utf-8 -*- """ Created on Mon Mar 12 15:58:17 2018 @author: Kata """ import numpy as np import pylab as pyl K_CEIL = 50 L_CEIL = 50 def cart2pol(x, y): """ Převede kartézské souřadnice na polární. """ @@ -82,81 +30,103 @@ def fakt(x): def A(k,l): return XX(k, l, 1) def B(k,l): return XX(k, l, -1) def C(k,l): return XX(k, l, 0) def D(k,l): return XX(k, l, -2) def X0(k, l, m, c): minisum = 2*l+k-4*m+c return (np.pi*N)**minisum / fakt(minisum) # One sum element of A - one sum element of B def A_B0(k, l, m): return X0(k, l, m, 1) - X0(k, l, m, -1) def C_D0(k, l, m): return X0(k, l, m, 0) - X0(k, l, m, -2) def sum_function(k, l, fun, minimal_c): m_max = np.floor((2*l+k+ minimal_c)/4)+1 vals = [fun(k, l, m) for m in range(int(m_max))] result = sum(vals) # napr. B ma konstantu c=(-1) (=minimal_c), A ma c=1, (-1) + 2 = 1. vals2 = [X0(k, l, m, minimal_c + 2) for m in range(int(m_max), int(np.floor((2*l+k+ minimal_c + 2)/4)+1))] return result + sum(vals2) # A: c = 1 # B: c = -1 # C: c = 0 # D: c = -2 def XX(k,l,c): ms=np.arange(0, np.floor((2*l+k+c)/4)+1) suma=0 for m in ms: suma += X0(k, l, m, c) return(suma) def YY(l, rho, X1, X2): ks=np.arange(0, K_CEIL) suma=0 for k in ks: suma += (-1)**k * (np.pi*N)**(2*l+k) * fakt(2*l+k+1) / fakt(k) / fakt(2*(2*l+1)+k) * (rho/rho_0)**(2*l+k+1) * (X1(k, l) - X2(k, l)) return(suma) def YY2(l, rho, DF, min_c): ks=np.arange(0, K_CEIL) suma=0 for k in ks: suma += (-1)**k * (np.pi*N)**(2*l+k) * fakt(2*l+k+1) / fakt(k) / fakt(2*(2*l+1)+k) * (rho/rho_0)**(2*l+k+1) * sum_function(k, l, DF, min_c) return(suma) def F(l,rho): return YY2(l, rho, A_B0, -1) def Z(fi, rho, fun): ls=np.arange(0, L_CEIL) suma=0 for l in ls: minisum = 2 * l + 1 temp = ( (np.sin(2 * minisum * fi)) / minisum) * fun(l,rho) suma = suma + temp return(suma) def EE(fi,rho): return Z(fi, rho, FF) def E(fi,rho): return Z(fi, rho, F) def FF(l,rho): return YY(l, rho, C, D) def G(fi,rho): temp = ( E(fi,rho)**2) + EE(fi,rho)**2 return(temp) def I(fi,rho): temp = 16*(N**2)*G(fi,rho) return(temp) @@ -201,4 +171,4 @@ def I(fi,rho): pixel_size=x/res #mm np.save('Image', Image) np.save('pixel_size', pixel_size) -
Apr 16, 2018 . 1 changed file with 0 additions and 0 deletions.There are no files selected for viewing
File renamed without changes. -
Apr 16, 2018 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,204 @@ # -*- coding: utf-8 -*- """ Created on Mon Mar 12 15:58:17 2018 @author: Kata """ import numpy as np """ m_max=60 ms = np.linspace(0,m_max,m_max+1) for poz, m in enumerate(ms): print(m) ns = np.linspace(0,m_max,m_max+1) for n in ns: print(n) pozi=0 for k in range(0,50): print(pozi, k) pozi=pozi+1 def C(x): temp = x**2 return (temp) hodnota_C = C(2) def A(x, a=2): temp = ((a+x)/x)*C(x) return (temp) hodnota_A = A(2,1) def D(x): temp = 2*x return (temp) hodnota_D = D(2) def B(x, c=1): d=2 temp = ((c+x)/d)*D(x) return (temp) hodnota_B = B(2,2) K=1 def F(x): temp = K*((A(x))**2+(B(x))**2) return(temp) hodnota_F = F(6) """ def cart2pol(x, y): """ Převede kartézské souřadnice na polární. """ rho = np.sqrt(x**2 + y**2) phi = np.arctan2(y, x) return(rho, phi) def pol2cart(rho, phi): """ Převede polární souřadnice na kartézské. """ x = rho * np.cos(phi) y = rho * np.sin(phi) return(x, y) def fakt(x): f=np.math.factorial(x) return(f) def A(k,l): ms=np.arange(0,np.floor((2*l+k+1)/4)+1) suma=0 for m in ms: minisum = 2*l+k-4*m+1 temp = (np.pi*N)**minisum / fakt(minisum) suma=suma+temp return(suma) def B(k,l): ms=np.arange(0,np.floor((2*l+k-1)/4)+1) suma=0 for m in ms: temp = ( (np.pi*N)**(2*l+k-4*m-1) ) / ( fakt(2*l+k-4*m-1) ) suma=suma+temp return(suma) def C(k,l): ms=np.arange(0,np.floor((2*l+k)/4)+1) suma=0 for m in ms: temp = ( (np.pi*N)**(2*l+k-4*m) ) / ( fakt(2*l+k-4*m) ) suma=suma+temp return(suma) def D(k,l): ms=np.arange(0,np.floor((2*l+k-2)/4)+1) suma=0 for m in ms: temp = ( (np.pi*N)**(2*l+k-4*m-2) ) / ( fakt(2*l+k-4*m-2) ) suma=suma+temp return(suma) def E(fi,rho): ls=np.arange(0,100) suma=0 for l in ls: temp = ( (np.sin(2*(2*l+1)*fi)) / (2*l+1) )*F(l,rho) suma = suma + temp return(suma) def F_old(l,rho): ks=np.arange(0,100) suma=0 for k in ks: temp = ( ( ( (-1)**k )*( (np.pi*N)**(2*l+k) )*( fakt(2*l+k+1) ))/ ( fakt(k) * fakt(2*(2*l+1)+k) ) ) * ((rho/rho_0)**(2*l+k+1)) * (A(k,l)-B(k,l)) suma = suma + temp return(suma) def F(l,rho): ks=np.arange(0,100) suma=0 for k in ks: suma += (-1)**k * (np.pi*N)**(2*l+k) * fakt(2*l+k+1) / fakt(k) / fakt(2*(2*l+1)+k) * (rho/rho_0)**(2*l+k+1) * (A(k, l) - B(k, l)) return(suma) def EE(fi,rho): ls=np.arange(0,100) suma=0 for l in ls: temp = ( (np.sin(2*(2*l+1)*fi)) / (2*l+1) )*FF(l,rho) suma = suma + temp return(suma) def FF(l,rho): ks=np.arange(0,100) suma=0 for k in ks: temp = ( ( ( (-1)**k )*( (np.pi*N)**(2*l+k) )*( fakt(2*l+k+1) ))/ ( fakt(k) * fakt(2*(2*l+1)+k) ) ) * ((rho/rho_0)**(2*l+k+1)) * (C(k,l)-D(k,l)) suma = suma + temp return(suma) def G(fi,rho): temp = ( E(fi,rho)**2) + EE(fi,rho)**2 return(temp) def I(fi,rho): temp = 16*(N**2)*G(fi,rho) return(temp) if __name__ == "__main__": ##################### Optické parametry ################################### Lambda = 400 # vlnová délkalaseru [nm] rho_0 = 3 # poloměr masky [mm] z = 200 # poloha stínítka za maskou [mm] x=y= 4 # rozměr stínítka [mm] res = 10 # rozlišení stínítka v osách x a y [počet bodů] ###################### Převody veličin #################################### rho_0 = rho_0*10**-3 Lambda = Lambda*10**-9 z = z*10**-3 x = x*10**-3 y = y*10**-3 xxs =np.linspace(-x/2, 0, res) yys =np.linspace(-y/2, 0, res) Image=np.zeros((res, res), dtype=float) ########################## Výpočty ######################################## N = (rho_0**2)/(Lambda*z)# Fresnelovo číslo N = 20 counter=0 for X, xx in enumerate(xxs): # výpočet intenzity for Y, yy in enumerate(yys): s=0 print (counter) counter=counter+1 rho, phi =cart2pol(xx,yy) s=I(phi,rho) Image[X,Y]=s #ERR=ERR.flatten() pixel_size=x/res #mm np.save('Image', Image) np.save('pixel_size', pixel_size)