#include <iostream>
#include <complex>
using namespace std;
//use units
#include <boost/typeof/std/complex.hpp>
#include <boost/units/cmath.hpp>
#include <boost/units/systems/si.hpp>
#include <boost/units/systems/si/prefixes.hpp>
using namespace boost::units;
using namespace boost::units::si;
//use units end

double integral(
	  const quantity<length> x	// x - pozycja dla ktorej licze jasnosc
	, const quantity<length> R	// R - odleglosc szczeliny od ekranu
	, const quantity<length> D	// D - szerokosc szczeliny
	, const quantity<length> lambda	// lambda - dlugosc fali
	, double steps			// steps - ilosc krokow calkowania == ilosc podzialow szer. D
)
{
    complex<double> i(0,1);
    complex<double> a(0);
    static const double pi = 3.14159265358979323846;
    quantity<wavenumber,complex<double> > k=2*pi/lambda;
    quantity<length> d=-0.5*D;   // d - pozycja na szczelinie
    quantity<length> dd=D/steps; // dd - krok calkowy
    // integrate from -d to d: exp(i*k*r)/r
    for( ; d<= 0.5*D ; d+=dd )
    {
        quantity<length> r=sqrt(pow<2>(R)+pow<2>(x-d));
        a += (dd*exp(i*k*r)/r).value();
    }
    return pow<2>(abs(a));
}

int main()
{
    quantity<length> rozmiar_ekranu         (  2048*7.4*micro*meters);
    quantity<length> odleglosc_do_ekranu    (  86.5   *centi*meters);
    quantity<length> rozmiar_szczeliny      (   0.515 *milli*meters);
    //quantity<length> dlugosc_fali           ( 543.0   *nano*meters); // 0.53mm
    quantity<length> dlugosc_fali           ( 632.8   *nano*meters); // 0.49mm

/*
    quantity<length> rozmiar_ekranu         (  27.8528*milli*meters);
    quantity<length> odleglosc_do_ekranu    (  86.5   *centi*meters);
    quantity<length> rozmiar_szczeliny      (   0.265  *milli*meters);
    //quantity<length> dlugosc_fali           ( 543.0   *nano*meters); // 0.29mm
    quantity<length> dlugosc_fali           ( 501.0   *nano*meters);
    //quantity<length> dlugosc_fali           ( 632.8   *nano*meters); // 0.27mm
 */

    quantity<length> s(rozmiar_ekranu*0.5);
    quantity<length> ds(0.01*milli*meter);
    for(quantity<length> x=-s ; x<=s ; x+=ds)
    {
        cout << x << " " << integral(x,odleglosc_do_ekranu,rozmiar_szczeliny,dlugosc_fali,1000) << "\n";
    }

    return false;
}