nemmtomm

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#include <X11/Xlib.h>
#include <assert.h>
#include <unistd.h>


Display *dpy;
Window w;
GC gc;





#define num double

void def(num * v1,num x,num y,num z) {v1[0]=x;v1[1]=y;v1[2]=z;}

void cross(num * v1,num * v2,num * v3)
{
    v1[0] = (v2[1]*v3[2])-(v3[1]*v2[2]);
    v1[1] = (v2[2]*v3[0])-(v3[2]*v2[0]);
    v1[2] = (v2[0]*v3[1])-(v3[0]*v2[1]);
};
num dot(num * v1,num *v2) {    num t = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];return t;}

void normalize(num *v1,num *v2)
{
    num t = sqrt(dot(v2,v2));
    v1[0] = v2[0]/t;
    v1[1] = v2[1]/t;
    v1[2] = v2[2]/t;
}




void pont(int x,int y,int szin)
{
    XSetForeground(dpy,gc,szin);
    XDrawPoint(dpy, w, gc, x,y);
}
void vonal(int x1,int y1,int x2,int y2,int szin)
{
    XSetForeground(dpy,gc,szin);
    XDrawLine(dpy, w, gc, x1,y1,x2,y2);
}
num sqr(num n)
{
    return n*n;
}




num t,dt,szog,m=1;
num pozicio[10][3], sebesseg[10][3], normal[3], I[3], B[3], F[3];
int sz=0;



void spin_wave()
{
    int a;
   

    dt=0.005;
    for(a=0;a<10;a++)
    {
        def(pozicio[a],100,0,0);
        def(sebesseg[a],0.0,10.0,1.0);
    }
    def(B,0,0,1);
    def(sebesseg[0],0.0,-10.0,1.0);
    def(sebesseg[2],0.0,-10.0,1.0);


    int szin=0;



    while(1)
    {
        for(a=0;a<10;a++)
        {
        pozicio[a][0] += sebesseg[a][0]*dt;
        pozicio[a][1] += sebesseg[a][1]*dt;
        pozicio[a][2] += sebesseg[a][2]*dt;

        normalize(normal,pozicio[a]);//r=100
        pozicio[a][0] = normal[0]*50.0;
        pozicio[a][1] = normal[1]*50.0;
        pozicio[a][2] = normal[2]*50.0;

        num h=dot(sebesseg[a],normal);
        sebesseg[a][0] -= normal[0]*h;
        sebesseg[a][1] -= normal[1]*h;
        sebesseg[a][2] -= normal[2]*h;
       
  
        cross(F,sebesseg[a],B);
        h=dt*0.02/m;
        sebesseg[a][0] += F[0]*h;
        sebesseg[a][1] += F[1]*h;
        sebesseg[a][2] += F[2]*h;

        normalize(sebesseg[a],sebesseg[a]);
        sebesseg[a][0] *= 10.0;
        sebesseg[a][1] *= 10.0;
        sebesseg[a][2] *= 10.0;


        sz++;
        if(sz>10)
        {
            sz=0;
            int x,y;
            cross(I,pozicio[a],sebesseg[a]);

            t=dot(B,I);
           
            if(t<0.0) szin=0xff0000;
            else szin=0x0000ff;
           
            for(y=0;y<50;y++)
            for(x=0;x<50;x++) pont(x+a*100,y,szin);
           
            pont(a*100+pozicio[a][0],200+pozicio[a][1],szin);
            pont(a*100+pozicio[a][0],300+pozicio[a][2],szin);
        }
     }   

        for(a=0;a<10;a++)
        {
        num h;
       
        if(a>0)
        {
            cross(I,pozicio[a-1],sebesseg[a-1]);
            cross(F,sebesseg[a],I);
            h=dt*0.0002/m;
            sebesseg[a][0] += F[0]*h;
            sebesseg[a][1] += F[1]*h;
            sebesseg[a][2] += F[2]*h;
            sebesseg[a-1][0] -= F[0]*h;
            sebesseg[a-1][1] -= F[1]*h;
            sebesseg[a-1][2] -= F[2]*h;
        }
        if(a<9)
        {
            cross(I,pozicio[a+1],sebesseg[a+1]);
            cross(F,sebesseg[a],I);
            h=dt*0.0002/m;
            sebesseg[a][0] += F[0]*h;
            sebesseg[a][1] += F[1]*h;
            sebesseg[a][2] += F[2]*h;
            sebesseg[a+1][0] -= F[0]*h;
            sebesseg[a+1][1] -= F[1]*h;
            sebesseg[a+1][2] -= F[2]*h;
        }
        normalize(sebesseg[a],sebesseg[a]);
        sebesseg[a][0] *= 10.0;
        sebesseg[a][1] *= 10.0;
        sebesseg[a][2] *= 10.0;

       
        }
    }
    XFlush(dpy);
}


int main()
{

    dpy  =  XOpenDisplay((0));
    w  =  XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0,0, 1200, 800, 0,0,0);
  
    XSelectInput(dpy, w, StructureNotifyMask);
    XMapWindow(dpy, w);
  
    gc  =  XCreateGC(dpy, w, 0, (0));
    XSetForeground(dpy,gc,0);
  
    for(;;) {    XEvent e;    XNextEvent(dpy, &e);    if (e.type  ==  MapNotify)break;    }
  
  
  spin_wave();
 
 
    XFlush(dpy);
    getchar();

    return 0;

}