rendered paste body#include <SDL/SDL.h>#include <GL/gl.h>#include <omp.h>#include <unistd.h>#include <stdlib.h>#include <stdio.h>#include <math.h>typedef float real;typedef unsigned int uint;static real epsilon = 5.96e-08;//static real epsilon = 1.11e-16; // doublereal randr(){ return 2.0 * rand() / RAND_MAX - 1.0;}struct simulation{ uint n; real *x, *y, *z; real *vx, *vy, *vz; real *nx, *ny, *nz; real *m;};struct simulation* sim_create(uint n, real radius, real mass){ struct simulation* sim; uint i; sim = malloc(sizeof(struct simulation)); sim->n = n; sim->x = malloc(n * sizeof(real)); sim->y = malloc(n * sizeof(real)); sim->z = malloc(n * sizeof(real)); sim->vx = malloc(n * sizeof(real)); sim->vy = malloc(n * sizeof(real)); sim->vz = malloc(n * sizeof(real)); sim->nx = malloc(n * sizeof(real)); sim->ny = malloc(n * sizeof(real)); sim->nz = malloc(n * sizeof(real)); sim->m = malloc(n * sizeof(real)); for (i=0; i<n; i++) { sim->m[i] = mass; sim->x[i] = randr()*radius; sim->y[i] = randr()*radius; sim->z[i] = randr()*radius; sim->vx[i] = 0.0; sim->vy[i] = 0.0; sim->vz[i] = 0.0; } return sim;}void sim_destroy(struct simulation *sim){ free(sim->x); free(sim->y); free(sim->z); free(sim->vx); free(sim->vy); free(sim->vz); free(sim->nx); free(sim->ny); free(sim->nz); free(sim->m); free(sim);}void sim_simulate(struct simulation *sim, real dt){ uint i;#pragma omp parallel for for (i=0; i<sim->n; i++) { real ax, ay, az; ax = 0.0; ay = 0.0; az = 0.0; uint j; for (j=0; j<sim->n; j++) { real dx, dy, dz; dx = sim->x[j] - sim->x[i]; dy = sim->y[j] - sim->y[i]; dz = sim->z[j] - sim->z[i]; real inv, f; inv = 1.0 / sqrt(dx*dx + dy*dy + dz*dz + epsilon*epsilon); inv = inv*inv*inv; f = sim->m[j]*inv; ax += f*dx; ay += f*dy; az += f*dz; } sim->nx[i] = sim->x[i] + dt*sim->vx[i] + 0.5*dt*dt*ax; sim->ny[i] = sim->y[i] + dt*sim->vy[i] + 0.5*dt*dt*ay; sim->nz[i] = sim->z[i] + dt*sim->vz[i] + 0.5*dt*dt*az; sim->vx[i] += dt*ax; sim->vy[i] += dt*ay; sim->vz[i] += dt*az; } for (i=0; i<sim->n; i++) { sim->x[i] = sim->nx[i]; sim->y[i] = sim->ny[i]; sim->z[i] = sim->nz[i]; }}void sim_draw(struct simulation *sim){ uint i; glColor3f(1.0, 1.0, 1.0); glBegin(GL_POINTS); for (i=0; i<sim->n; i++) { glVertex3f(sim->x[i], sim->y[i], sim->z[i]); } glEnd();}int main(int argc, char **argv){ int n = 0; int frames = 500; int continuous = 0; real radius = 100.0; real mass = 1.0; if (argc < 2) { fprintf(stderr, "%s - num [-r radius] [-m mass] [-f frames] [-c]\n", argv[0]); return 1; } n = atoi(argv[1]); char c; while ((c = getopt(argc, argv, "r:m:f:c")) != -1) { switch(c) { case 'r': radius = atof(optarg); break; case 'm': mass = atof(optarg); break; case 'f': frames = atoi(optarg); break; case 'c': continuous = 1; break; } } struct simulation *sim = sim_create(n, radius, mass); if (SDL_Init(SDL_INIT_VIDEO) < 0) { fprintf(stderr, "Could not init video: %s\n", SDL_GetError()); return 2; } if (SDL_SetVideoMode(800, 600, 32, SDL_OPENGL) == 0) { fprintf(stderr, "Could not set video mode: %s\n", SDL_GetError()); return 2; } glClearColor(0, 0, 0, 0); glViewport(0, 0, 800, 600); glMatrixMode(GL_PROJECTION); glLoadIdentity(); glOrtho(-radius, radius, -radius, radius, -radius*10, radius*10); glMatrixMode(GL_MODELVIEW); glLoadIdentity(); glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT); int done = 0; while (frames-- || continuous) { glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT); sim_simulate(sim, 0.1); sim_draw(sim); SDL_PumpEvents(); if (SDL_GetKeyState(NULL)[SDLK_ESCAPE]) break; if (SDL_GetKeyState(NULL)[SDLK_SPACE]) { sim_destroy(sim); sim = sim_create(n, radius, mass); } SDL_GL_SwapBuffers(); } sim_destroy(sim); return 0;}