#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <X11/Xlib.h>


Display *dpy;
Window win;
GC gc;
double radian=(180.0/M_PI);


void pixel(int x,int y,int color)
{
    XSetForeground(dpy,gc,color);
    XDrawPoint(dpy, win, gc, x,y);
}
double sqr(double n) // x^2
{
    return n*n;
}
double doublerand() //random szam  0.0-tol 1.0-ig
{
    return (double)(rand()%10000)/10000.0;
}



struct vec2d
{
    double x,y;
};
void add_amp(vec2d *v,double phase,double a1,double a2)
{
    v->x += sin(phase)*a1;
    v->y += cos(phase)*a2;
};
void add_polarizer(vec2d *v,double phase)
{
    v->x *= cos(phase);
    v->y *= sin(phase);
};
double probalbility(vec2d *v)
{
    return (sqr(v->x) + sqr(v->y));
}


int main()
{
    dpy = XOpenDisplay(0);
    win = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0,0, 800, 550, 0,0,0);
    
    XSelectInput(dpy, win, StructureNotifyMask);
    XMapWindow(dpy, win);
    gc = XCreateGC(dpy, win, 0, 0);
    
    for(;;) { XEvent e; XNextEvent(dpy, &e); if (e.type == MapNotify) break; }



    for(int x=0;x<400;x++)    
    {
        pixel(x,500,0x008800);
        pixel(x,300,0x008800);
        pixel(x,100,0x008800);
    }
    
    

    for(int ds_x=0;ds_x<400;ds_x++)//Ds position  -+4mm    
    {
        int photon_counter=0;
        int maxphoton=400;
        
        for(int w=0;w<maxphoton;w++)//slit wide and max number of photon
        {
            double phase_a=M_PI*2*doublerand();
            double ds_distance=1250.0-420.0;//mm  125-42 cm
            double dp_distance=980.0; //98 cm
            double wavelength=702.2e-6;//mm      e-9m
            double k=2.0*M_PI/wavelength;
            vec2d amp_dsdp;
            amp_dsdp.x=0;
            amp_dsdp.y=0;



            ds_distance+=0.0;
            double hole_dst=0.2;//0.2
            double hole_wide=0.2;//200 micrometer wide   
            


//Dp side
            dp_distance+=0.0;//hatrebb?
            add_amp(&amp_dsdp,    phase_a + dp_distance*k,0.5,0.5);
            add_polarizer(&amp_dsdp,3*45.0/radian);//polarizator Dp elott


//Ds side            
            int hole_side=(rand()%1000)/500;//vagy a) vagy b) mert a masik Dp fele megy
            double ds_pos=4.0*(double)(ds_x-200)/200.0;//+-4mm     Ds position
//Ds side a)Hole
            double holex=hole_dst/2.0 + hole_wide*(double)w/maxphoton;//hole
            double t=sqrt(sqr(ds_pos - holex) + sqr(ds_distance));
            if(hole_side==0) add_amp(&amp_dsdp,    phase_a + t*k,0.5,0.5);

//Ds side b)Hole
            holex=-hole_dst/2.0 - hole_wide*(double)w/maxphoton;
            t=sqrt(sqr(ds_pos - holex) + sqr(ds_distance));            
            if(hole_side==1) add_amp(&amp_dsdp,    phase_a + t*k,0.5,-0.5);// forditott forgasi irany


            if((probalbility(&amp_dsdp))>doublerand())
                photon_counter+=1;
        }
        pixel(ds_x,500-photon_counter,0xffff00);
    }

    XFlush(dpy);
    getchar();
        
    return 0;
}