#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn
from scipy.integrate import quad
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def icart(x, y, delta_l, l_m):
    lower_bound = 1 - delta_l/(2*l_m)
    upper_bound = 1 + delta_l/(2*l_m)
    prefactor = 4E4/(np.pi**2 * (x**2 + y**2))
    return prefactor * quad(lambda labda: (jn(1, np.pi*np.sqrt(x**2 + y**2)/labda))**2,
                            lower_bound, upper_bound)[0]

def cartesian_norm_airy(xvals, yvals, delta_l, l_m):
    lower_bound = 1 - delta_l/(2*l_m)
    upper_bound = 1 + delta_l/(2*l_m)

    X, Y = np.meshgrid(xvals, yvals)
    prefactor = 4E4/(np.pi**2 * (X**2+Y**2))
    returnarray = np.zeros_like(X)
    for rowidx, rowval in enumerate(X):
        for colidx, colval in enumerate(rowval):
            returnarray[rowidx][colidx] = icart(X[rowidx][colidx], Y[rowidx][colidx], delta_l, l_m)
    return X, Y, returnarray


x = np.linspace(0,5,100)-2.5+0.0001
y = x
z = cartesian_norm_airy(x, y, 0.3, 45)