GLUT ->  OpenGL Utility Toolkit
OpenGLのお手軽なCラッパ。

***
インスコ@debian

$su
%apt-get install libghc6-glut-dev

// GLUTで重力多体問題

#include <GL/glut.h>
#include <cmath>
#include <vector>
#include <map>
#include <iostream>
#include <stdio.h>

using namespace std;

int T = 1000;
double H = 0.01;
double H2=H*H;
double HH = H/2;
double H2H = H2/2;
double G = 8.0;
int rate = 10;
int gcount = 0;
int N=0;
const double eps = 0.0001;

///////////////////////////////
double xa, xb, ya, yb;
double xc, yc, wx, wy;


void putPoint(double x, double y){
  glVertex2d((x-xc)*wx, (y-yc)*wy);
}

struct particle{
  vector<double> q, p, x, xp;
  double m;
  particle(){
    q = vector<double>(2,0);
    p = vector<double>(2,0);
    x = vector<double>(2,0);
    xp = vector<double>(2,0);
    m = 0;
  }

  inline void setQ(double a, double b){ q[0] = a; q[1] = b;  }
  inline void setP(double a, double b){ p[0] = a; p[1] = b;  }
  inline void setA(double a, double b){ x[0] = a; x[1] = b;  }
  inline void setM(double M){ m = M;  }

  inline void stepmove(){
    if (gcount==0) cout << " " << q[0] << " " << q[1];
    q[0] += H*p[0]+H2H*x[0]; q[1] += H*p[1]+H2H*x[1];
    p[0] += HH*(x[0]+xp[0]); p[1] += HH*(x[1]+xp[1]);

    xp[0] = x[0]; xp[1] = x[1];
    x[0] = x[1] = 0; // 後処理

  }
};


vector<particle> obj;

vector<double> force(particle* a, particle* b){ //力計算
  vector<double> ans(2,0);
  vector<double> r(2,0);
  r[0] = a->q[0]-b->q[0];
  r[1] = a->q[1]-b->q[1];

  double r3 = pow(r[0]*r[0]+r[1]*r[1]+eps, -1.5);
  
  ans[0] = -G*r3*r[0];
  ans[1] = -G*r3*r[1];
  return ans;
}



void putPoints(){
  glClear(GL_COLOR_BUFFER_BIT);
  glBegin(GL_POINTS);
  for(int i=0; i<N; ++i){
    putPoint(obj[i].q[0], obj[i].q[1]);
  }
  glEnd();
  glFlush();  
}

void stepforward(){ //全体の時間発展処理
  vector<double> a(2);
  for(int i=0; i<N; ++i){
    for(int j=i+1; j<N; ++j){
      a = force(&obj[i], &obj[j]);
      obj[i].x[0] += a[0]*obj[j].m;
      obj[i].x[1] += a[1]*obj[j].m;

      obj[j].x[0] -= a[0]*obj[i].m;
      obj[j].x[1] -= a[1]*obj[i].m;
    }
  }

  for(int i=0; i<N; ++i){
    obj[i].stepmove();
  }

  if(gcount==0){
    putPoints();
    cout << endl;
  }
}


void display(void){
  putPoints();
}


void init(void){
  glClearColor(0., 0., 0., 1.);
}

void mouse(int button, int state, int x, int y){
  int turn = 0;
  for(int t=0; t<T; ++t){
    stepforward();
    (++gcount) %= rate;
  }
}

int main(int argc, char *argv[]){
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_RGBA);
  glutCreateWindow(argv[0]);
  glutDisplayFunc(display);
  glutMouseFunc(mouse);
  init();

  cin >> xa >> ya >> xb >> yb;
  wx = 2/(xa-xb);
  wy = 2/(ya-yb);
  xc = (xa+xb)/2;
  yc = (ya+yb)/2;

  cin >> T >> H >> rate >> G;
  H2 = H*H;
  HH = H/2;
  H2H = H2/2;
  cin >> N;

  string str;
  for(int i=0; i<N; ++i){
    double Mass, x, y, vx, vy;
    cin >> Mass >> x >> y >> vx >> vy;
    particle T;
    T.setQ(x, y); T.setP(vx, vy); T.setM(Mass);
    obj.push_back(T);
    cerr << obj[i].m << " " << obj[i].q[0] << " " << obj[i].q[1]
     << " " << obj[i].p[0] << " " << obj[i].p[1] << endl;
  }

  glutMainLoop();

}

//GLUTでもりもり正規分布

#include <GL/glut.h>
#include <cmath>
#include <vector>
#include <iostream>
#include <boost/random.hpp>

using namespace std;
using namespace boost;

int N = 1024; // columns of histogram
unsigned long long FullLength = 1 << 20;
double LL = -0.9, RL = 0.9;
double UL = 0.9, DL = -0.9;



mt19937 gen(199999999);
normal_distribution<> dst( .0, .1 );
variate_generator<mt19937&, normal_distribution<>
> normi( gen, dst );


struct bar{
  unsigned long long H;

  bar(unsigned long long height = 0){
    H = height;
  }

  unsigned long long operator ++(){
    ++H;
    return H;
  }

  unsigned long long operator ++(int){
    ++H;
    return H-1;
  }
  
  double Height(){
    return (H+.0)/FullLength;
  }
};

vector<bar> histo;


void display(void){
  glClear(GL_COLOR_BUFFER_BIT);
  glBegin(GL_LINES);

  for(int i=0; i<N; ++i){
    glVertex2d(LL+(RL-LL)/N*i, DL);
    glVertex2d(LL+(RL-LL)/N*i, histo[i].Height()*(UL-DL)+DL);
    
  }
  glEnd();
  glFlush();  
}


void init(void){
  glClearColor(0., 0., 0., 1.);
}

void mouse(int button, int state, int x, int y){
  for(;;){
    for(int i=0; i<100000; ++i){
      int i = N*(normi()+.5);
      if(0<=i && i < N) ++(histo[i]);
    }
    display();
  }
}


int main(int argc, char *argv[]){

  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_RGBA);
  glutCreateWindow(argv[0]);
  glutDisplayFunc(display);
  glutMouseFunc(mouse);
  init();

  histo = vector<bar>(N, bar());

  glutMainLoop();
}