/*
   Diffraction.

   Kode til at udregne intensitetsmønstret for double slit eksperimentet. Hvor hver
   spalter også har en endelig bredde.

   I teorien skelner man mellem  Frauenhofer grænsen (langt fra spalte til skærmen),
   mens i Fresnell grænsen er skærmen tæt på spalten. Koden kan simulerer begge.

   Amplituden for en vej med længde d fra hul til skærm for en lyskilde med bølgelænge
   lambda er   A = exp(I*2.0*M_PI*phi)/d    hvor phi = d/lambda er fasen.
   eksponential funktionen bølgens komplekse amplitude, mens jeg dividerer med d fordi
   amplituden falder med afstanden fra lyskilden for en sphærisk bølge. Dette er en meget
   lille effekt.

   Jeg approximerer en endelig bred spalte som N punktkilder fordelt langs spalten
   som er en tilnærmelse af integralet over spaltens bredde som anvendes i Frauenhofer
   teorien.

   matlab:
   data=importdata("doubleslit.dat"); plot(data(:,1),data(:,2),'.');
*/

#include <iostream>
#include <fstream>
#include <cmath>
#include <complex>

using namespace std;

int main()
{
   ofstream fo("doubleslit.dat");

   const double m=1;     // m
   const double mm=1e-3; // m
   const double nm=1e-9; // m

   double lambda=633*nm;  // Wavelength of laser
   double L=0.05*m;          // Distance slit to screen
   double h=1*mm;         // Distance between slits
   double w=0.1*mm;       // width of slits
   double N=100;          // Points summed per slit.

   // Range to plot
   double ymin=-100*mm;
   double ymax= 100*mm;
   complex<double> I(0,1);           // Imaginary number i

   for (double y=ymin ; y<ymax ; y+=0.1*mm )
     {
       complex<double> a(0,0);
       double xh=0;        // x position of hole
       double yh;          // y position of hole
       double yc;          // centre position of hole
       double d,phi;       // distance, phase

// Calculate amplitudes from N point sources in top hole
       yc=h/2;
       for (int i=0; i<N; i++)
         {
           double yh=yc+w*(i+0.5-N/2)/N;                 // (i+0.5-N/2)/N  number from -1/2 ... +1/2
           d=sqrt( (L-xh)*(L-xh) + (y-yh)*(y-yh) );
           phi=d/lambda;
           a+=exp(I*2.0*M_PI*phi)/(N*d);                 // 1/N due to integral approximated as sum.
         }

// Calculate amplitudes from N point sources in bottom hole
       yc=-h/2;
       for (int i=0; i<N; i++)
         {
           double yh=yc+w*(i+0.5-N/2)/N;
           d=sqrt( (L-xh)*(L-xh) + (y-yh)*(y-yh) );
           phi=d/lambda;
           a+=exp(I*2.0*M_PI*phi)/(N*d);
         }

       double Intensity=norm(a);                        // norm(z) = x^2 + y^2
       fo << y << "\t" << Intensity << "\n";
   }

   fo.close();
}