Some calculation

vypocet.py

# -*- 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í.
    """
    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):
    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)



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)