//g++ x.cpp -lX11


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


#include <map>
using namespace std;


#define igaz true
#define hamis false
#define visszater return



Display *dpy;
Window win;
GC gc;

void pont(int x,int y,int szin)
{
    XSetForeground(dpy,gc,szin);
    XDrawPoint(dpy, win, gc, x,y);
}
void vonal(int x1,int y1,int x2,int y2,int szin)
{
    XSetForeground(dpy,gc,szin);
    XDrawLine(dpy, win, gc, x1,y1, x2,y2);
}


#define gyok sqrt


//#define float1 double
#define float1 float


struct float2
{
float1 x,y;

float2() {x=y=0;};
float2(float1 x2,float1 y2) {x=x2;y=y2;}
};

float1 skalar(float2 v1,float2 v2) { return (v1.x*v2.x + v1.y*v2.y );}

float1 hossz(float2 v1) {return gyok(skalar(v1,v1));}

float2 normalt(float2 v1)
{
    float1 q=hossz(v1);
    v1.x/=q;
    v1.y/=q;
    visszater v1;
}
float1 frand(float1 n)
{
    return n*(float1)(rand()%20000)/20000.0;
}









struct adat
{
    float2 v1;
    int tipus;
};
map<float1,adat> pontok;



void metszet(float2 &v1,float2 &v2,float2 &u1,float2 &u2,int i)
{
    float2 ir1,ir2, nor,tmp,metszespont;


    ir1.x=v2.x-v1.x;
    ir1.y=v2.y-v1.y;
    float1 h1=hossz(ir1);    
    ir1.x /= h1;
    ir1.y /= h1;
    
    ir2.x=u2.x-u1.x;
    ir2.y=u2.y-u1.y;
    float1 h2=hossz(ir2);    
    ir2.x /= h2;
    ir2.y /= h2;
    nor.x= ir2.y;
    nor.y=-ir2.x;    



    tmp.x=v1.x-u1.x;//u1
    tmp.y=v1.y-u1.y;
    float1 tav= (0.0- skalar(tmp,nor)) / skalar(ir1,nor);
    if(tav>=0.0)
    if(tav<=h1)
    {
        metszespont.x = v1.x + ir1.x*tav;
        metszespont.y = v1.y + ir1.y*tav;
        
        tmp.x = metszespont.x - u1.x;
        tmp.y = metszespont.y - u1.y;
        float1 tav2=skalar(tmp,ir2);
        
        if(tav2>=0.0)
        if(tav2<=h2)
        {
            adat ad;
            ad.v1=metszespont;
            ad.tipus=i;
            
            pontok[tav]=ad;
        }
    }
}
float2 sarkok[2][5];
int prizman_belul[600*600];
float1 amplitudo[600*600] ,sebesseg[600*600];

bool belul(float2 v1)
{
    for(int i=0;i<2;i++)
    {
        pontok.clear();
        metszet(v1,sarkok[i][4],sarkok[i][0],sarkok[i][1],i);
        metszet(v1,sarkok[i][4],sarkok[i][1],sarkok[i][2],i);
        metszet(v1,sarkok[i][4],sarkok[i][3],sarkok[i][2],i);
        metszet(v1,sarkok[i][4],sarkok[i][0],sarkok[i][3],i);
        if(pontok.size()==0)
                        visszater igaz;
    }
    visszater hamis;
}



int main()
{
 
dpy = XOpenDisplay(0);
win = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0,0, 600, 600, 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; }



    sarkok[0][0]=float2(135+10,95);
    sarkok[0][1]=float2(200-10,95);
    sarkok[0][2]=float2(200,145);
    sarkok[0][3]=float2(135,145);

    sarkok[0][4].x=(sarkok[0][0].x+sarkok[0][1].x+sarkok[0][2].x+sarkok[0][3].x)/4.0;
    sarkok[0][4].y=(sarkok[0][0].y+sarkok[0][1].y+sarkok[0][2].y+sarkok[0][3].y)/4.0;
    
    
    sarkok[1][0]=float2( 70,150);
    sarkok[1][1]=float2(200,150);
    sarkok[1][2]=float2(200-10,200);
    sarkok[1][3]=float2( 70+10,200);

    sarkok[1][4].x=(sarkok[1][0].x+sarkok[1][1].x+sarkok[1][2].x+sarkok[1][3].x)/4.0;
    sarkok[1][4].y=(sarkok[1][0].y+sarkok[1][1].y+sarkok[1][2].y+sarkok[1][3].y)/4.0;
    


 

int u=256,u2=8;// 2^u2 = u
for(int y=0;y<u;y+=1)
for(int x=0;x<u;x+=1)
{
    if(belul(float2(x,y)))
    {
        prizman_belul[y*u+x]=1;
    }
    else
        prizman_belul[y*u+x]=0;
        
    pont(x,y,prizman_belul[y*u+x]*250);
}

float1 time=0;
float1 D=0.2,dt=0.5;//0.3

while(1)
{

    for(int y=200;y<222;y+=1)
    for(int x=10;x<30;x+=1)
    {
        int index=(y<<u2)+x+y-200;

        amplitudo[index]=sin(time*0.4);
    }
    for(int y=0;y<u;y+=1)
    for(int x=0;x<u;x+=1)
    {
        int index=(y<<u2)+x;
    
        if(amplitudo[index]> 1.0) amplitudo[index]=1.0;
        if(amplitudo[index]<-1.0) amplitudo[index]=-1.0;
        amplitudo[index] += sebesseg[index]*dt;
    }
    for(int y=1;y<u-1;y+=1)
    for(int x=1;x<u-1;x+=1)
    {
        int index=(y<<u2)+x;
        int index1=(((y+1)&(u-1))<<u2)+((x));
        int index2=(((y-1)&(u-1))<<u2)+((x));
        int index3=(((y))<<u2)   +((x+1)&(u-1));
        int index4=(((y))<<u2)   +((x-1)&(u-1));
    
        float1 alap=amplitudo[index1]+amplitudo[index2]+amplitudo[index3]+amplitudo[index4];
        alap*=(1.0/4.0);
        if(prizman_belul[index])
            sebesseg[index] -= (amplitudo[index]-alap)*D*dt;
        else
            sebesseg[index] -= (amplitudo[index]-alap)*D*dt*1.5;
    }
    for(int y=0;y<u;y+=1)
    for(int x=0;x<u;x+=1)
    {
        int index=(y<<u2)+x; // <<u2 == *u

#if 1
        int szin=127+amplitudo[index]*120.0;
#else
        int szin=amplitudo[index]*250.0;
#endif
        if(szin<0) szin=0;
        if(prizman_belul[index]) szin<<=16;
        pont(x,y,szin);
    }

    time+=dt;
    XFlush(dpy);
}
    XFlush(dpy);
getchar();

return 0;
}