#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=550;
        int maxwide=20;
        
        for(int p=0;p<maxphoton;p++)// max number of photon
        {
            double phase_a=M_PI*2*doublerand();
            double ds_distance=1250.0-420.0;//mm  125-42 cm
            double ds_distance2=420.0+980.0;//mm  visszafele a restol BBO-ig
            double dp_distance=980.0; //98 cm
            double wavelength=702.2e-6;//mm      e-9m
            double k=2.0*M_PI/wavelength;
            vec2d amp_dp,amp_ds;
            amp_dp.x=0;
            amp_dp.y=0;
            amp_ds.x=0;
            amp_ds.y=0;



            ds_distance+=0.0;
            double hole_dist=0.2;//0.2
            double hole_wide=0.2;//200 micrometer wide   
            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

            for(int w=0;w<maxwide;w++)//slit wide
            {
            double hole1x=hole_dist/2.0 + hole_wide*(double)w/maxwide;//hole
            double dist1=sqrt(sqr(ds_pos - hole1x) + sqr(ds_distance));
            double dist1b=sqrt(sqr(ds_pos - hole1x) + sqr(ds_distance2));//backward
            double hole2x=-hole_dist/2.0 - hole_wide*(double)w/maxwide;
            double dist2=sqrt(sqr(ds_pos - hole2x) + sqr(ds_distance));            
            double dist2b=sqrt(sqr(ds_pos - hole2x) + sqr(ds_distance2));            


            double rot_direction=1.0;// +-1
//Dp side
            dp_distance+=0.0;//hatrebb?
            //ha Ds normaliranyu, ez forditott
            if(hole_side==0) add_amp(&amp_dp,    phase_a + dp_distance*k + dist1b*k + M_PI,0.5, 0.5*rot_direction);
            if(hole_side==1) add_amp(&amp_dp,    phase_a + dp_distance*k + dist2b*k + M_PI,0.5, 0.5);//normal            
//Ds side            
            if(hole_side==0) add_amp(&amp_ds,    phase_a + dist1*k ,0.5, 0.5);// normal
            if(hole_side==1) add_amp(&amp_ds,    phase_a + dist2*k ,0.5, 0.5*rot_direction);// forditott forgasi irany
            }
//            add_polarizer(&amp_dp,(0*45)/radian);//0 1 2*45 polarizator Dp elott
            amp_ds.x+=amp_dp.x;
            amp_ds.y+=amp_dp.y;

            amp_ds.x/=maxwide;//normalizalas
            amp_ds.y/=maxwide;
        
            if((probalbility(&amp_ds))>doublerand())//mindig csak egy foton van
                photon_counter+=1;
        }
        pixel(ds_x,500-photon_counter,0xffff00);
    }

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