#include <math.h>
#include <X11/Xlib.h>
typedef long double skalar;
struct vektor
{
skalar x,y,z;
vektor() {x=0;y=0;z=0;};
vektor(skalar x_,skalar y_,skalar z_) {x=x_;y=y_;z=z_;};
vektor operator + (vektor v) {vektor v_;v_.x=x+v.x;v_.y=y+v.y;v_.z=z+v.z;return v_;};
vektor operator - (vektor v) {vektor v_;v_.x=x-v.x;v_.y=y-v.y;v_.z=z-v.z;return v_;};
vektor operator * (skalar s) {vektor v_;v_.x=x*s;v_.y=y*s;v_.z=z*s;return v_;};
vektor operator / (skalar s) {vektor v_;v_.x=x/s;v_.y=y/s;v_.z=z/s;return v_;};
};
skalar skalar_szorzat(vektor v1,vektor v2) { return (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);}
skalar hossz(vektor v1) {return sqrtl(skalar_szorzat(v1,v1));}
vektor normalt(vektor v1) { return (v1/hossz(v1));}
vektor vektor_szorzat(vektor va,vektor vb)
{
vektor v_;
v_.x=(va.y*vb.z)-(vb.y*va.z);
v_.y=(va.z*vb.x)-(vb.z*va.x);
v_.z=(va.x*vb.y)-(vb.x*va.y);
return v_;
};
Display *dpy;
Window w;
GC gc;
void pont(int x,int y,int color)
{
XSetForeground(dpy,gc,color);
XDrawPoint(dpy, w, gc, x,y);
}
vektor a(skalar u,skalar v,skalar R)
{
vektor v_;
v_.x=R*cos(u)*sin(v);
v_.y=R*sin(u)*sin(v);
v_.z=R*cos(v);
return v_;
}
int main()
{
dpy = XOpenDisplay(0);
w = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0,0, 800, 600, 0,0,0);
XSelectInput(dpy, w, StructureNotifyMask);
XMapWindow(dpy, w);
gc = XCreateGC(dpy, w, 0, 0);
for(;;) { XEvent e; XNextEvent(dpy, &e); if (e.type == MapNotify) break; }
int i,ii=10000;
skalar radian=M_PI/180;
skalar u=22*radian;
skalar v=40*radian;
skalar du=0.1*radian;
skalar dv=0.05*radian;
skalar z=du*0.01;
skalar R=200;
vektor p1=a(u,v,R);
vektor v1=a(u+du,v+dv,R)-a(u,v,R);
for(i=0;i<ii;i++)
{
vektor n1=normalt(p1);
v1=v1-n1*skalar_szorzat(n1,v1);
p1=p1+v1;
pont(400+(int)p1.x,300+(int)p1.y,0xffff00);
}
p1=a(u,v,R);
v1=a(u+du,v+dv,R)-a(u,v,R);
for(i=0;i<ii;i++)
{
vektor au = a(u+z,v ,R)-a(u,v,R);
vektor av = a(u ,v+z,R)-a(u,v,R);
vektor au_du = a(u+z+du ,v ,R) - a(u+du ,v ,R);
vektor au_dv = a(u+z ,v +dv,R) - a(u ,v+dv,R);
vektor av_du = a(u +du ,v+z ,R) - a(u+du ,v ,R);
vektor av_dv = a(u ,v+z+dv,R) - a(u ,v+dv,R);
skalar E = skalar_szorzat(au,au);
skalar F = skalar_szorzat(au,av);
skalar G = skalar_szorzat(av,av);
skalar u_uu = skalar_szorzat((au_du-au)/du,au);
skalar u_uv = skalar_szorzat((au_dv-au)/dv,au);
skalar u_vu = skalar_szorzat((av_du-av)/du,au);
skalar u_vv = skalar_szorzat((av_dv-av)/dv,au);
skalar v_uu = skalar_szorzat((au_du-au)/du,av);
skalar v_uv = skalar_szorzat((au_dv-au)/dv,av);
skalar v_vu = skalar_szorzat((av_du-av)/du,av);
skalar v_vv = skalar_szorzat((av_dv-av)/dv,av);
skalar g = 1/(E*G-F*F);
// vektor normal = vektor_szorzat(au,av); skalar g=1/skalar_szorzat(normal,normal);
u_uu*=E;
u_uv*=G;
u_vu*=G;
u_vv*=F;
v_uu*=E;
v_uv*=G;
v_vu*=G;
v_vv*=F;
skalar T0_00 = (u_uu + u_uu - u_uu)*g;
skalar T0_01 = (v_uu + u_uv - v_uu)*g;
skalar T0_10 = (u_vu + u_vu - u_uv)*g;
skalar T0_11 = (v_vu + u_vv - v_uv)*g;
skalar T1_00 = (u_uv + v_uu - u_vu)*g;
skalar T1_01 = (v_uv + v_uv - v_vu)*g;
skalar T1_10 = (u_vv + v_vu - u_vv)*g;
skalar T1_11 = (v_vv + v_vv - v_vv)*g;
skalar ddu= -T0_00*du*du - T0_01*du*dv - T0_10*dv*du - T0_11*dv*dv;
skalar ddv= -T1_00*du*du - T1_01*du*dv - T1_10*dv*du - T1_11*dv*dv;
u += du;
v += dv;
du += ddu;
dv += ddv;
p1=a(u,v,R);
pont(400+(int)p1.x,300+(int)p1.y,0xff0000);
}
XFlush(dpy);
getchar();
return 0;
}