API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 1st, 2020
139
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 9.69 KB | None | 0 0
raw download clone embed print report
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdlib.h>
  4.  
  5. #include <GL\glut.h>
  6. #include <math.h>
  7. #include <time.h>
  8. #include "struktury.h"
  9.  
  10. #define YS 1.0f
  11. #define PI 3.1415
  12.  
  13.  
  14. double rot=0.0f;
  15. int solver=2;
  16.  
  17. //--------------------------------------------------------------
  18. //
  19. // do rozwiazania uklad rownan
  20. // {dr/dt=v
  21. // {d2r/dt2=F/m
  22. //
  23. // Solvery:
  24. // solveEuler
  25. // solveMidpoint
  26. // solveRK4
  27. //
  28. // podobnie jak w pendulum
  29. //--------------------------------------------------------------
  30.  
  31. void derivatives(Wektor *in, Wektor *out, Punkt *p){
  32. // dr/dt =v
  33. // d2r/dt2=F/m
  34. out[0]=in[1];
  35. out[1]=p->f*(1.0/p->m);
  36. }
  37.  
  38.  
  39.  
  40. void solveEuler(Punkt *p, float dt){
  41.  
  42.  
  43. Wektor stan_poczatkowy[2];
  44. Wektor stan_aktualny[2];
  45.  
  46.  
  47. stan_poczatkowy[0]=p->r;
  48. stan_poczatkowy[1]=p->v;
  49.  
  50. derivatives(stan_poczatkowy, stan_aktualny, p);
  51.  
  52.  
  53. p->r=p->r+stan_aktualny[0]*dt;
  54. p->v=p->v+stan_aktualny[1]*dt;
  55.  
  56. }
  57.  
  58. void solveMidPoint(Punkt *p, float dt){
  59. //r,v,f ->Wektor
  60. Wektor k1[2],k2[2];
  61. Wektor y_k[2];
  62. Wektor y_h[2];
  63. Wektor yout[2];
  64.  
  65. y_h[0]=p->r;
  66. y_h[1]=p->v;
  67. derivatives(y_h,yout,p);
  68.  
  69. k1[1]=yout[1];
  70. k1[0]=yout[0];
  71. y_k[0]=y_h[0]+k1[0]*0.5*dt;
  72. y_k[1]=y_h[1]+k1[1]*0.5*dt;
  73. derivatives(y_k,yout,p);
  74. k2[1]=yout[1];
  75. k2[0]=yout[0];
  76. p->r=p->r+k2[0]*dt;
  77. p->v=p->v+k2[1]*dt;
  78. }
  79.  
  80. void solveRK4(Punkt *p, float dt){
  81. //r,v,f ->Wektor
  82. Wektor y_in[2];
  83. Wektor y_temp[2];
  84. Wektor k[4][2];
  85. y_in[0] = p->r;
  86. y_in[1] = p->v;
  87. derivatives(y_in, k[0], p);
  88. y_temp[0] = y_in[0] + k[0][0] * 0.5*dt;
  89. y_temp[1] = y_in[1] + k[0][1] * 0.5*dt;
  90. derivatives(y_temp, k[1], p);
  91. y_temp[0] = y_in[0] + k[1][0] * 0.5*dt;
  92. y_temp[1] = y_in[1] + k[1][1] * 0.5*dt;
  93. derivatives(y_temp, k[2], p);
  94. y_temp[0] = y_in[0] + k[2][0] * dt;
  95. y_temp[1] = y_in[1] + k[2][1] * dt;
  96. derivatives(y_temp, k[3], p);
  97. p->r = p->r + (k[0][0] + 2* k[1][0] + 2* k[2][0] + k[3][0] )* dt * (1.0/6);
  98. p->v = p->v + (k[0][1] + 2 * k[1][1] + 2 * k[2][1] + k[3][1]) * dt * (1.0/6);
  99. }
  100.  
  101.  
  102. //--------------------------------------------------------------
  103. // RESZTA KODU
  104. //--------------------------------------------------------------
  105.  
  106.  
  107.  
  108.  
  109. struct SferaN{
  110. Wektor r1;
  111. float r;
  112. float t;
  113. float color[3];
  114. SferaN *right;
  115. };
  116.  
  117.  
  118.  
  119. Wektor g(0,-1.0,0);
  120. Punkt *root=NULL;
  121. SferaN *sroot=NULL;
  122. int loop=0;
  123. SferaN *SphereAlloc(float R, Wektor r1, float t, float c[3])
  124. {
  125. SferaN *tmp;
  126. if (!(tmp=(SferaN*)malloc(sizeof(SferaN))))
  127. return NULL;
  128.  
  129. tmp->right=NULL;
  130. tmp->r=R;
  131. tmp->r1=r1;
  132. tmp->t=t;
  133. tmp->color[0]=c[0];
  134. tmp->color[1]=c[1];
  135. tmp->color[2]=c[2];
  136.  
  137. return tmp;
  138. }
  139.  
  140. void SphereTest(SferaN *s, Punkt *p)
  141. {
  142. float d;
  143. Wektor n;
  144. float z;
  145. d=(s->r1-p->r).len();
  146. if (d-s->r<0)
  147. {
  148. n=(s->r1-p->r);
  149. n.norm();
  150. z=d-s->r;
  151. p->r=p->r+n*z;
  152.  
  153. Wektor vs,vn;
  154. vn=n*(n*p->v);
  155. vs=p->v-vn;
  156. p->v=(vs-vn*s->t);
  157. }
  158.  
  159.  
  160. }
  161.  
  162. void AddSphere(SferaN *ro, float R, Wektor r1, float t, float c[3])
  163. {
  164. SferaN *tmp;
  165. for (tmp=ro; tmp->right!=NULL; tmp=tmp->right);
  166.  
  167. tmp->right=SphereAlloc(R,r1,t,c);
  168. }
  169.  
  170.  
  171. Punkt *PointAlloc(float m, int flaga, Wektor r, Wektor v)
  172. {
  173. Punkt *tmp;
  174. if (!(tmp=(Punkt*)malloc(sizeof(Punkt))))
  175. return NULL;
  176.  
  177. tmp->m=m;
  178. tmp->flag=flaga;
  179. tmp->r=r;
  180. tmp->v=v;
  181. tmp->right=NULL;
  182. return tmp;
  183. }
  184.  
  185. void AddPoint(Punkt *ro, float m, int flag, Wektor r, Wektor v)
  186. {
  187. Punkt *tmp;
  188. for (tmp=ro; tmp->right!=NULL; tmp=tmp->right);
  189.  
  190. tmp->right=PointAlloc(m,flag,r,v);
  191. }
  192.  
  193. Wektor W_Euler(Wektor f, float h)
  194. {
  195. return (f*h);
  196. }
  197.  
  198. void calcDeriv(Punkt *p){
  199.  
  200. }
  201.  
  202. void Solver(Punkt *ro, float dt)
  203. {
  204. //0 euler, 1 midpoint, 2 RK4
  205. Punkt *tmp;
  206. for (tmp=ro;tmp!=NULL;tmp=tmp->right)
  207. {
  208. switch (solver){
  209. case 0:
  210. solveEuler(tmp,dt);
  211. break;
  212. case 1:
  213. solveMidPoint(tmp,dt);
  214. break;
  215. case 2:
  216. solveRK4(tmp,dt);
  217. break;
  218. default:
  219. solveEuler(tmp,dt);
  220. }
  221. /*
  222. tmp->dv=W_Euler(tmp->f*(1/tmp->m),dt);
  223. tmp->v=tmp->v+tmp->dv;
  224.  
  225. tmp->dr=tmp->v*dt;
  226. tmp->r=tmp->r+tmp->dr;
  227. */
  228.  
  229. /*
  230. Wektor a=tmp->f*(1.0/tmp->m);
  231. //--------------
  232.  
  233. Wektor k1=W_Euler(a,dt);
  234. Wektor k2=W_Euler(a+0.5*k1,dt);
  235. Wektor k3=W_Euler(a+0.5*k2,dt);
  236. Wektor k4=W_Euler(a+k3,dt);
  237. Wektor kk=(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4);
  238. tmp->v=tmp->v+kk;
  239. tmp->r=tmp->r+tmp->v*dt;
  240. */
  241.  
  242. }
  243. }
  244.  
  245. void Sily(Punkt *ro)
  246. {
  247. Punkt *tmp;
  248.  
  249. for (tmp=ro;tmp!=NULL; tmp=tmp->right)
  250. {
  251. tmp->f=g*tmp->m;
  252. }
  253. }
  254.  
  255.  
  256.  
  257.  
  258. void AnimateScene(void)
  259. {
  260. glClearColor(0.0f,0.0f,0.0f,0.0f);
  261.  
  262. //if (loop==0)
  263. glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  264.  
  265. Punkt *tmp;
  266. //Sily(root);
  267. Solver(root,0.008);
  268. Sily(root);
  269. // glRotatef(0.1,0.0f,1.0f,0.0f);
  270.  
  271. rot+=0.01;
  272. //if (rot>36) rot=0.0f;
  273. glPointSize(5);
  274. glDisable(GL_LIGHTING);
  275. glDisable(GL_LIGHT0);
  276. for(tmp=root;tmp!=NULL;tmp=tmp->right)
  277. {
  278. //glBegin(GL_POINTS);
  279.  
  280. //glColor3f(rand()/(float)RAND_MAX,rand()/(float)RAND_MAX,rand()/(float)RAND_MAX);
  281. glColor3f(1,1,1);
  282. // glVertex3f(tmp->r.x/1.0,tmp->r.y/1.0,tmp->r.z/1.0);
  283.  
  284. glPushMatrix();
  285. glTranslatef(tmp->r.x/1.0,tmp->r.y/1.0,tmp->r.z/1.0);
  286. glutSolidSphere(0.011,5,5);
  287. glPopMatrix();
  288. SferaN *stmp;
  289.  
  290. for(stmp=sroot;stmp!=NULL;stmp=stmp->right)
  291. {
  292.  
  293. SphereTest(stmp,tmp);
  294. }
  295. //glVertex2f(0,0);
  296. //printf("%f %f\n",tmp->r.x, tmp->r.y);
  297. // glEnd();
  298.  
  299. glBegin(GL_LINES);
  300.  
  301. glColor3f(1,0,0);
  302. double x1,x2,y1,y2,z1,z2;
  303. x1=tmp->r.x;
  304. y1=tmp->r.y;
  305. x2=tmp->v.x;
  306. y2=tmp->v.y;
  307. z1=tmp->r.z;
  308. z2=tmp->v.z;
  309.  
  310.  
  311. glVertex3f(x1,y1,z1);
  312. glVertex3f(x1+(x2)*0.05,y1+(y2)*0.05,z1+(z2)*0.05);
  313. glEnd();
  314.  
  315. }
  316. glEnable(GL_LIGHTING);
  317. glEnable(GL_LIGHT0);
  318. SferaN *stmp;
  319. for(stmp=sroot;stmp!=NULL;stmp=stmp->right)
  320. {
  321. glPushMatrix();
  322. glTranslatef(stmp->r1.x,stmp->r1.y,stmp->r1.z);
  323. glColor3fv(stmp->color);
  324.  
  325. glutSolidSphere(stmp->r,50,50);
  326. glPopMatrix();
  327. }
  328.  
  329.  
  330. // printf("****\n");
  331.  
  332. glFlush();
  333. glutSwapBuffers();
  334. loop++;
  335.  
  336.  
  337. if (loop>2000)
  338. {
  339. double vx,vy,vz;
  340. double zz=0.01;
  341.  
  342. for(tmp=root;tmp!=NULL;tmp=tmp->right)
  343. {
  344. vx=0.5-rand()/(float)RAND_MAX;
  345. vy=1.0-rand()/(float)RAND_MAX;
  346. vz=0.5-rand()/(float)RAND_MAX;
  347. vz=0.0f;
  348. vy*=YS;
  349. vy+=zz;
  350. zz+=0.001;
  351. tmp->r.x=0;
  352. tmp->r.y=0;
  353. tmp->v=Wektor(vx,vy,vz);
  354. }
  355. loop=0;
  356. }
  357. }
  358.  
  359. void InitGraphics()
  360. {
  361. GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 };
  362. GLfloat mat_shininess[] = { 90.0 };
  363. GLfloat light_position[] = {1.0, 1.0, 1.0, 0.0};
  364.  
  365. glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
  366. glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
  367. glLightfv(GL_LIGHT0, GL_POSITION, light_position);
  368. glEnable(GL_LIGHTING);
  369. glEnable(GL_LIGHT0);
  370. glDepthFunc(GL_LEQUAL);
  371. glEnable(GL_DEPTH_TEST);
  372.  
  373.  
  374. glShadeModel(GL_SMOOTH);
  375.  
  376.  
  377.  
  378.  
  379. glEnable(GL_COLOR_MATERIAL);
  380.  
  381.  
  382. }
  383.  
  384.  
  385. void myReshape(GLsizei w, GLsizei h)
  386. {
  387. glViewport(0, 0, w, h);
  388. glMatrixMode(GL_PROJECTION);
  389. glLoadIdentity();
  390.  
  391. double p1=1.0f;
  392.  
  393. if(w <= h)
  394. glOrtho(-p1,p1,-p1*(GLfloat)h/(GLfloat)w,p1*(GLfloat)h/(GLfloat)w,-10.0,10.0);
  395. else
  396. glOrtho(-p1*(GLfloat)w/(GLfloat)h,p1*(GLfloat)w/(GLfloat)h,-p1,p1,-10.0,10.0);
  397.  
  398. glMatrixMode(GL_MODELVIEW);
  399. glLoadIdentity();
  400. }
  401. void keyboard(unsigned char key, int x, int y){
  402. switch (key){
  403. case 'e':
  404. solver=0;
  405. printf("Solver --> Euler\n");
  406. break;
  407. case 'm':
  408. solver=1;
  409. printf("Solver --> MidPoint\n");
  410. break;
  411. case 'r':
  412. solver=2;
  413. printf("Solver --> RK4\n");
  414. break;
  415.  
  416. }
  417. }
  418. void idle(void)
  419. {
  420. glutPostRedisplay();
  421. }
  422.  
  423. int main(int argc, char *argv[])
  424. {
  425. double vx,vy,vz;
  426. int i;
  427. srand(time(NULL));
  428.  
  429. double zz=0.01;
  430.  
  431. for (i=1; i<620; i++)
  432. {
  433. vx=0.5-rand()/(float)RAND_MAX;
  434. vy=1.0-rand()/(float)RAND_MAX;
  435. vz=0.5-rand()/(float)RAND_MAX;
  436. vz=0.0f;
  437.  
  438. vy*=YS;
  439. vy+=zz;
  440. zz+=0.001;
  441.  
  442. if (!root)
  443. root=PointAlloc(1,0,Wektor(0,0,0),Wektor(vx,vy,vz));
  444. else
  445. AddPoint(root,1,0,Wektor(0,0,0),Wektor(vx,vy,vz));
  446. }
  447. zz=-0.2;
  448.  
  449. for (i=1; i<140; i++)
  450. {
  451.  
  452. AddPoint(root,10,0,Wektor(zz,1,0),Wektor(0,0,0));
  453. zz+=0.01;
  454. }
  455.  
  456. float c1[] = {0,0,1};
  457. float c2[] = {0,1,0};
  458. float c3[] = {1,1,0};
  459. float c4[] = {1,0,1};
  460. float c5[] = {0.6,0.5,1};
  461. float c6[] = {0.6,0.5,0.3};
  462.  
  463. sroot=SphereAlloc(0.5,Wektor(0,-0.5,0),0.9,c1);
  464. AddSphere(sroot,0.1,Wektor(-0.5,0.5,0),0.8,c2);
  465. AddSphere(sroot,0.1,Wektor(0.5,0.5,0),0.3,c3);
  466. AddSphere(sroot,0.1,Wektor(1,0.0),0.3,c4);
  467. AddSphere(sroot,0.15,Wektor(0,0.7,0),0.9,c5);
  468. AddSphere(sroot,0.2,Wektor(-1,-0.2,0),1.0,c6);
  469.  
  470. Sily(root);
  471.  
  472. glutInit(&argc, argv);
  473. glutInitWindowSize(600,600);
  474. glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
  475. glutCreateWindow("GLUT example");
  476.  
  477.  
  478.  
  479. InitGraphics();
  480.  
  481. glutDisplayFunc(AnimateScene);
  482. glutIdleFunc(idle);
  483. glutKeyboardFunc(keyboard);
  484. glutReshapeFunc(myReshape);
  485. glutMainLoop();
  486.  
  487. return 0;
  488. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!