#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);
}
void line(int x1,int y1,int x2,int y2,int color)
{
    XSetForeground(dpy,gc,color);
    XDrawLine(dpy, win, gc, x1,y1,x2,y2);
}
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;
    
    vec2d() {x=0;y=0;}
};
void add_amp(vec2d *v,double phase,double ax,double ay)
{
    v->x += sin(phase)*ax;
    v->y += cos(phase)*ay;
};
void add_polarizer(vec2d *v,double phase)
{
    v->x *= sin(phase);
    v->y *= cos(phase);
};
double dot(vec2d *v1,vec2d *v2)
{
    return v1->x*v2->x + v1->y*v2->y;
}
void add_quarterwaveplate(vec2d *v,double phase,double dist,double qwp,double ax,double ay)
{
    vec2d axis_fast,axis_slow,input_wave;
    
    add_amp(&axis_fast,qwp,1.0,1.0);
    add_amp(&axis_slow,qwp+90.0/radian,1.0,1.0);
    add_amp(&input_wave,phase,1.0,1.0);

    double amp_fast=dot(&axis_fast,&input_wave);
    double amp_slow=dot(&axis_slow,&input_wave);
    
    amp_fast=cos(dist)*amp_fast;
    amp_slow=sin(dist)*amp_slow;//-+90 phase shift
    
    v->x += (axis_fast.x*amp_fast + axis_slow.x*amp_slow)*ax;
    v->y += (axis_fast.y*amp_fast + axis_slow.y*amp_slow)*ay;
};
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;
        int maxwide=10;

        int qwplate_on=1;
        int eraser_on=1;
        double eraser_alpha=45.0/radian;
        
        for(int p=0;p<maxphoton;p++)// max number of photon
        {
#if 0
            double phase_a=M_PI*2*doublerand();
            double phase_b=phase_a+M_PI/2;
#else            
            double phase_a=0; //polarized lightsource
            if(doublerand()>0.5) phase_a=M_PI/2;
            double phase_b=phase_a+M_PI/2;
            if(doublerand()>0.5) phase_b=phase_a-M_PI/2;
#endif


            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_dp,amp_ds;



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

            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 hole2x=-hole_dist/2.0 - hole_wide*(double)w/maxwide;
                double dist1=sqrt(sqr(ds_pos - hole1x) + sqr(ds_distance));
                double dist2=sqrt(sqr(ds_pos - hole2x) + sqr(ds_distance));            


                add_amp(&amp_dp,phase_b ,0.5,0.5);
                add_amp(&amp_dp,phase_b ,0.5,0.5);

                if(qwplate_on)
                {
                    if(eraser_on) phase_a=eraser_alpha;
                    add_quarterwaveplate(&amp_ds,phase_a , dist1*k,-45.0/radian,0.5, 0.5);
                    add_quarterwaveplate(&amp_ds,phase_a , dist2*k, 45.0/radian,0.5, 0.5);
                }
                else
                {
                    add_amp(&amp_ds,phase_a +dist1*k ,0.5,0.5);
                    add_amp(&amp_ds,phase_a +dist2*k ,0.5,0.5);
                }
            }
            if(eraser_on) add_polarizer(&amp_dp,eraser_alpha);// polarizer before Dp

            
            amp_dp.x/=maxwide;//normalized
            amp_dp.y/=maxwide;
            amp_ds.x/=maxwide;
            amp_ds.y/=maxwide;


            if((probalbility(&amp_dp))>doublerand())
            if((probalbility(&amp_ds))>doublerand())
                photon_counter+=1;
        }

        pixel(ds_x,500-photon_counter,0xffff00);
    }

/*
        double dist1=0.0;
        for(int p=0;p<1000;p++)//quarterwaveplate
        {
            vec2d amp_dp,amp_ds;
            add_quarterwaveplate(&amp_ds,0.0/radian , dist1*M_PI/100.0,45.0/radian,0.5, 0.5);

            line(200,200,200+200.0*amp_ds.x,200+200.0*amp_ds.y,0x00ff00);
            dist1+=5.0;
            
            XFlush(dpy);
            getchar();
        }
*/
    XFlush(dpy);
    getchar();
        
    return 0;
}