Last active
April 21, 2018 09:19
-
-
Save matejak/c53ae7adb5ff016b8d9b6b3f593d3333 to your computer and use it in GitHub Desktop.
Some calculation
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 characters
| # -*- 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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment