http://www.commonsensescience.org/spinning_ring.html
http://en.wikipedia.org/wiki/Toroidal_ring_model
#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 gyok sqrt
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);
}
float sqr(float n) { return n*n;}
int main()
{
int i,j,x;
float sugar[2];
dpy = XOpenDisplay((0));
w = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0,0, 1000, 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; }
for(x=0;x<7000;x+=100)
{
int db=0;
float ampx=0,ampy=0,ido;
float forras[]={100,50};
float cel[]={700,100};
for(j=-600;j<601;j++)
{
float felulet[]={(x+j)/10.0,400.0};//felulet pontjai x +j kornyezo pontok osszessege
sugar[0]=felulet[0]-forras[0]; sugar[1]=felulet[1]-forras[1];
float tavolsag=gyok(sugar[0]*sugar[0] + sugar[1]*sugar[1]);//forrastol a feluletig
sugar[0]=cel[0]-felulet[0]; sugar[1]=cel[1]-felulet[1];//felulettol a celig
tavolsag+=gyok(sugar[0]*sugar[0] + sugar[1]*sugar[1]);
float fazis=tavolsag*M_PI/10.0;
ampx+=cos(fazis);//forgo vektor, nem komplex szam!
ampy+=sin(fazis);
if(j==0) ido=tavolsag;//sebesseg 1
db++;
}
float m=sqr((float)db);
float amp=(ampx*ampx + ampy*ampy);//amp negyzete, nincs gyok!
int y=(int)(500.0-amp*100.0/m);;
pont(x/10,y,0xff);
pont(x/10,400,0x00ff00);
pont(x/10,1000-(int)(ido/2.0),0x0088ff);
int szin=(int)(amp*255.0/m);
vonal(x/10,400, (int)forras[0],(int)forras[1],szin);
vonal(x/10,400, (int)cel[0],(int)cel[1],szin);
}
for(x=0;x<7000;x+=100)
{
int db=0;
float ampx=0,ampy=0,ido;
float forras[]={100,50};
float cel[]={700,600};
for(j=-600;j<601;j++)
{
float felulet[]={(x+j)/10.0,400.0};
sugar[0]=felulet[0]-forras[0]; sugar[1]=felulet[1]-forras[1];
float tavolsag=gyok(sugar[0]*sugar[0] + sugar[1]*sugar[1]);
sugar[0]=cel[0]-felulet[0]; sugar[1]=cel[1]-felulet[1];
float tavolsag2=gyok(sugar[0]*sugar[0] + sugar[1]*sugar[1]);
float fazis=tavolsag*M_PI/10.0+tavolsag2*M_PI/8.0;
ampx+=cos(fazis);
ampy+=sin(fazis);
if(j==0) ido=tavolsag+tavolsag2;
db++;
}
float m=sqr((float)db);
float amp=(ampx*ampx + ampy*ampy);
int y=(int)(500.0-amp*100.0/m);;
pont(x/10,y,0xff0000);
pont(x/10,400,0x00ff00);
pont(x/10,1000-(int)(ido/2.0),0xff8800);
int szin=(int)(amp*255.0/m)<<16;
vonal(2+x/10,400, (int)forras[0],(int)forras[1],szin);
vonal(2+x/10,400, (int)cel[0],(int)cel[1],szin);
}
pont(100,50,0xffffff);
pont(700,100,0xffffff);
pont(700,600,0xffffff);
XFlush(dpy);
getchar();
return 0;
}
Szólj hozzá!
#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,meresi_szog,m=1;
num pozicio[3], sebesseg[3], normal[3], I[3], B[3], F[3];
int stat[2],szamlalo=0,sz3=0,meres,B_valtas;
void set()
{
def(pozicio,100,0,0);
szog=M_PI/2*((float)(rand()%10000))/10000.0;
def(sebesseg,0.0,cos(szog)*10.0,sin(szog)*10.0);
def(B,0,0,1);
szamlalo=0;
B_valtas=100000+((rand()%1000)*100+(rand()%1000));
meres=B_valtas+50000+(rand()%1000)*50+(rand()%1000);
}
void stern_gerlach()
{
int a;
for(a=0;a<360;a+=10)
{
meresi_szog=a;
t=sin((meresi_szog/2.0)*M_PI*2/360.0);
t*=t;
printf("%f fok %f \n",meresi_szog,t);
pont(a,500-(int)(t*200.0),0xff0000);
dt=0.02;
set();
stat[0]=0;
stat[1]=0;
int szin=0;
while(1)
{
pozicio[0] += sebesseg[0]*dt;
pozicio[1] += sebesseg[1]*dt;
pozicio[2] += sebesseg[2]*dt;
normalize(normal,pozicio);//r=100
pozicio[0] = normal[0]*100.0;
pozicio[1] = normal[1]*100.0;
pozicio[2] = normal[2]*100.0;
num h=dot(sebesseg,normal);
sebesseg[0] -= normal[0]*h;
sebesseg[1] -= normal[1]*h;
sebesseg[2] -= normal[2]*h;
normalize(sebesseg,sebesseg);
sebesseg[0] *= 10.0;
sebesseg[1] *= 10.0;
sebesseg[2] *= 10.0;
cross(F,sebesseg , B );
h=dt*0.02/m;
sebesseg[0] += F[0]*h;
sebesseg[1] += F[1]*h;
sebesseg[2] += F[2]*h;
szamlalo++;
if(szamlalo==meres)//impulzusmomentum iranya a magneses terhez kepest
{
cross(I,pozicio,sebesseg);
t=dot(B,I);
if(t<0.0) stat[0]++;//ellentetes
else stat[1]++;//megegyezo
set();
}
else
if(szamlalo==B_valtas)//magneses ter elfordul
{
szog=meresi_szog*M_PI*2/360.0;
def(B,sin(szog),0.0,cos(szog));
}
#if 0
//beloves
pont(
500+(int)pozicio[1],
500+(int)pozicio[0],szin++);
cross(I,pozicio,sebesseg);
vonal(
500+(int)pozicio[1],
500+(int)pozicio[0],
500+(int)pozicio[1]+I[1]*0.1,
500+(int)pozicio[0]+I[0]*0.1,
0xff00ff );
#endif
sz3--;
if(sz3<0)//neha statisztika kiiras
{
sz3=1000000;
if(stat[1]) printf("%d/%d/%d %d %\n",stat[0],stat[1],stat[0]+stat[1],stat[0]*100/(stat[0]+stat[1]));
}
if((stat[0]+stat[1])>100) break;
}
num st=stat[0]*100/(stat[0]+stat[1]);
pont(a,500-(int)(st*2.0),0xffff00);
XFlush(dpy);
}
}
int main()
{
dpy = XOpenDisplay((0));
w = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0,0, 1000, 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; }
stern_gerlach();
XFlush(dpy);
getchar();
return 0;
}
Szólj hozzá!
#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;
}
