A polarizer nem mindig jó ebben a formában.
Javítva
//#define _COMPLEXNUMBER
enum{
prgst_noQWP,
prgst_QWPpresent,
prgst_ERASERpresent,
prgst_ERASERpresentClassic,
prgst_backwardcausality
};
// int prgstate=prgst_noQWP;
// int prgstate=prgst_QWPpresent;
// int prgstate=prgst_ERASERpresent;
// int prgstate=prgst_ERASERpresentClassic;
int prgstate=prgst_backwardcausality;
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <X11/Xlib.h>
Display *dpy;
Window win;
GC gc;
double iradian=(M_PI/180.0);
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);
}
void initX11()
{
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; }
}
double sqr(double n) // x^2
{
return n*n;
}
double doublerand() //rand 0.0-1.0
{
return (double)(rand()%10000)/10000.0;
}
#ifdef _COMPLEXNUMBER
#define vec2d complex
struct complex
{
double real,img;
complex() {real=img=0;};
complex(double r,double i) {real=r;img=i;};
void set(double r,double i) {real=r;img=i;};
void set(complex c) {real=c.real;img=c.img;};
void conjugate() {img*=-1.0;}
complex operator +(complex c)
{
complex e;
e.real = this->real + c.real;
e.img = this->img + c.img;
return e;
}
complex operator -(complex c)
{
complex e;
e.real = this->real - c.real;
e.img = this->img - c.img;
return e;
}
complex operator *(complex c)
{
complex e;
e.real = this->real*c.real - this->img*c.img;
e.img = this->img*c.real + this->real*c.img;
return e;
}
complex operator /(complex c)
{
complex e;
double denominator=c.real*c.real + c.img*c.img;
e.real = (this->real * c.real + this->img * c.img)/denominator;
e.img = (this->img * c.real - this->real * c.img)/denominator;
return e;
}
complex operator *(double c)
{
complex e;
e.real = this->real*c;
e.img = this->img*c;
return e;
}
};
void add_amp(complex *v,double phase,double ax)
{
*v = *v + complex(sin(phase)*ax,cos(phase)*ax);
};
void add_polarizer2(complex *v,double phase)
{
complex pol;
add_amp(&pol,phase,1.0);
complex amp=(*v)*pol; amp.img=0; //v'=pol><pol|v>
*v = pol*amp;
};
double probability(complex c1,complex c2)
{
c1.conjugate();// <c1|
complex P=c1*c2; //<c1|c2>
return P.real;
}
void add_quarterwaveplatel2c(vec2d *v,double phase,double dist,double qwp_angle,double ax)
{
complex axis_fast,axis_slow,input_wave;
add_amp(&axis_fast,qwp_angle,1.0);
add_amp(&axis_slow,qwp_angle+90.0*iradian,1.0);
add_amp(&input_wave,phase,1.0);
complex amp_fast=axis_fast*input_wave; amp_fast.img=0;
complex amp_slow=axis_slow*input_wave; amp_slow.img=0;
amp_fast=amp_fast*complex(cos(dist)*ax,0.0);
amp_slow=amp_slow*complex(sin(dist)*ax,0.0);//-+90 phase shift
*v=*v+axis_fast*amp_fast + axis_slow*amp_slow;
};
#else //###########################################################################################
struct vec2d
{
double x,y;
vec2d() {x=0;y=0;}
vec2d operator *(double c)
{
vec2d u;
u.x=x*c;
u.y=y*c;
return u;
}
};
double dot(vec2d *v1,vec2d *v2)
{
return v1->x*v2->x + v1->y*v2->y;
}
void add_amp(vec2d *v,double phase,double ax)
{
v->x += sin(phase)*ax;
v->y += cos(phase)*ax;
};
void add_polarizer2(vec2d *v,double phase)
{
vec2d pol;
add_amp(&pol,phase,1.0);
double amp=dot(v,&pol); //v'=pol><pol|v>
v->x = pol.x*amp;
v->y = pol.y*amp;
};
double probability(vec2d v1,vec2d v2)
{
return dot(&v1,&v2);
}
//linear to circular
void add_quarterwaveplatel2c(vec2d *v,double phase,double dist,double qwp_angle,double ax)
{
vec2d axis_fast,axis_slow,input_wave;
add_amp(&axis_fast,qwp_angle,1.0);
add_amp(&axis_slow,qwp_angle+90.0*iradian,1.0);
add_amp(&input_wave,phase,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)*ax;
}
#endif //###########################################################################################
int main()
{
initX11();
for(int i=100;i<500;i+=20) line(0,i,400,i,0x005500);
line(0,500,400,500,0x00aa00);
line(0,300,400,300,0x00aa00);
line(0,100,400,100,0x00aa00);
double eraser_angle=(-45.0*1.0)*iradian;//+-45
for(int ds_x=0;ds_x<200;ds_x++)//Ds position -+4mm
{
double P=0,Pds=0,Pdp=0;
double step=0.01;
for(double p=0;p<M_PI*2.0;p+=step)
{
vec2d amp_dp,amp_ds;
double photon_pol_a, photon_pol_b;
// if prgstate!=prgst_ERASERpresentClassic AND prgstate!=prgst_backwardcausality
photon_pol_a=0;// QM for BBO
photon_pol_b=0+M_PI/2;
if(rand()%100<50) { photon_pol_a=M_PI/2;photon_pol_b=0; }
if(prgstate==prgst_ERASERpresentClassic)
{
photon_pol_a=p;//classic
photon_pol_b=p+M_PI/2;
}
if(prgstate==prgst_backwardcausality)
{
photon_pol_a=eraser_angle;//source modified by Dp polarizer (eraser)
photon_pol_b=eraser_angle+M_PI/2;
if(rand()%100<50) { photon_pol_a=eraser_angle+M_PI/2;photon_pol_b=eraser_angle; }
}
double dphase1=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;//wave number
double hole_dist=0.2;//0.2
double hole_wide=0.2;//0.2mm = 200 micrometer wide
double ds_pos=4.0*(double)(ds_x-100)/100.0;//+-4mm Ds position
int slit_wide=10;
for(int w=0;w<slit_wide;w++)
{
double hole1x=hole_dist/2.0 + hole_wide*(double)w/slit_wide;//hole A
double hole2x=-hole_dist/2.0 - hole_wide*(double)w/slit_wide;//hole B
double dist1=sqrt(sqr(ds_pos - hole1x) + sqr(ds_distance));
double dist2=sqrt(sqr(ds_pos - hole2x) + sqr(ds_distance));
if(prgstate==prgst_noQWP)
{
add_amp(&_dp,photon_pol_a ,1.0);
add_amp(&_ds,photon_pol_b +dist1*k +dphase1,0.5);
add_amp(&_ds,photon_pol_b +dist2*k +dphase1,0.5);
}
else
//if(prgstate==prgst_QWPpresent || prgstate==prgst_ERASERpresent ||prgstate==prgst_backwardcausality || prgst_ERASERpresentClassic)
{
add_amp(&_dp,photon_pol_a ,1.0);
add_quarterwaveplatel2c(&_ds,photon_pol_b, dist1*k +dphase1,-45.0*iradian,0.5);
add_quarterwaveplatel2c(&_ds,photon_pol_b, dist2*k +dphase1, 45.0*iradian,0.5);
}
}
double normalization_factor=1.0/(double)slit_wide;
amp_ds=amp_ds*normalization_factor;
amp_dp=amp_dp*normalization_factor;
if(prgstate==prgst_ERASERpresent || prgstate==prgst_backwardcausality ||prgst_ERASERpresentClassic)
add_polarizer2(&_dp,eraser_angle);// polarizer(eraser) before Dp
P+=probability(amp_dp,amp_dp)*probability(amp_ds,amp_ds)*step;
Pds+=probability(amp_ds,amp_ds)*step;
Pdp+=probability(amp_dp,amp_dp)*step;
}
P/=(M_PI*2.0); pixel(ds_x*2,501-P*400.0,0xffff00);//520
Pds/=(M_PI*2.0); pixel(ds_x*2,500-Pds*400.0,0xff0000);
Pdp/=(M_PI*2.0); pixel(ds_x*2,502-Pdp*400.0,0x00ff);
}
XFlush(dpy);
getchar();
return 0;
}