diff --git a/examples_ptx/nbody/plummer.h b/examples_ptx/nbody/plummer.h
new file mode 100644
index 00000000..bb26f017
--- /dev/null
+++ b/examples_ptx/nbody/plummer.h
@@ -0,0 +1,138 @@
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include "vector3.h"
+
+struct Plummer{
+	std::vector<double> mass;
+	std::vector<dvec3> pos, vel;
+	Plummer(
+			unsigned long n, 
+			unsigned int  seed = 19810614, 
+			const char   *filename = "plummer.dat") 
+		: mass(n), pos(n), vel(n)
+	{
+		{
+			std::ifstream ifs(filename);
+			if(!ifs.fail()){
+				unsigned long ntmp, stmp;
+				ifs.read((char *)&ntmp, sizeof(unsigned long));
+				ifs.read((char *)&stmp, sizeof(unsigned long));
+				if(n == ntmp  && seed == stmp){
+					ifs.read((char *)&mass[0], n*sizeof(double));
+					ifs.read((char *)& pos[0], n*sizeof(dvec3));
+					ifs.read((char *)& vel[0], n*sizeof(dvec3));
+					if(!ifs.fail()){
+						fprintf(stdout, "plummer : read from %s\n", filename);
+					}
+					return;
+				}
+			}
+		}
+		srand(seed);
+		unsigned long i = 0;
+		while(i < n){
+			double X1 = my_rand();
+			double X2 = my_rand();
+			double X3 = my_rand();
+			double R = 1.0/sqrt( (pow(X1,-2.0/3.0) - 1.0) );
+			if(R < 100.0) {
+				double Z = (1.0 - 2.0*X2)*R;
+				double X = sqrt(R*R - Z*Z) * cos(2.0*M_PI*X3);
+				double Y = sqrt(R*R - Z*Z) * sin(2.0*M_PI*X3);
+
+				double Ve = sqrt(2.0)*pow( (1.0 + R*R), -0.25 );
+
+				double X4 = 0.0; 
+				double X5 = 0.0;
+
+				while( 0.1*X5 >= X4*X4*pow( (1.0-X4*X4), 3.5) ) {
+					X4 = my_rand(); X5 = my_rand(); 
+				} 
+
+				double V = Ve*X4;
+
+				double X6 = my_rand();
+				double X7 = my_rand();
+
+				double Vz = (1.0 - 2.0*X6)*V;
+				double Vx = sqrt(V*V - Vz*Vz) * cos(2.0*M_PI*X7);
+				double Vy = sqrt(V*V - Vz*Vz) * sin(2.0*M_PI*X7);
+
+				double conv = 3.0*M_PI/16.0;
+				X *= conv; Y *= conv; Z *= conv;    
+				Vx /= sqrt(conv); Vy /= sqrt(conv); Vz /= sqrt(conv);
+
+				double M = 1.0;
+				mass[i] = M/n;
+
+				pos[i][0] = X;
+				pos[i][1] = Y;
+				pos[i][2] = Z;
+
+				vel[i][0] = Vx;
+				vel[i][1] = Vy;
+				vel[i][2] = Vz;
+
+				/*
+				tmp_i = ldiv(i, 256);
+				if(tmp_i.rem == 0) printf("i = %d \n", i);
+				*/
+
+				ldiv_t tmp_i = ldiv(i, n/64);
+
+				if(tmp_i.rem == 0) {
+					printf(".");
+					fflush(stdout); 
+				}
+				i++; 
+			}		
+		} // while (i<n)
+		double mcm = 0.0;
+
+		double xcm[3], vcm[3];
+		for(int k=0;k<3;k++) {
+			xcm[k] = 0.0; vcm[k] = 0.0;
+		} /* k */
+
+		for(i=0; i<n; i++) {
+			mcm += mass[i];
+			for(int k=0;k<3;k++) {
+				xcm[k] += mass[i] * pos[i][k]; 
+				vcm[k] += mass[i] * vel[i][k]; 
+			} /* k */ 
+		}  /* i */
+		for(int k=0;k<3;k++) {
+			xcm[k] /= mcm; vcm[k] /= mcm;
+		} /* k */
+
+		for(i=0; i<n; i++) {
+			for(int k=0;k<3;k++) {
+				pos[i][k] -= xcm[k]; 
+				vel[i][k] -= vcm[k]; 
+			} /* k */ 
+		} /* i */ 
+		printf("\n");
+		{
+			std::ofstream ofs(filename);
+			if(!ofs.fail()){
+				unsigned long ntmp = n;
+				unsigned long stmp = seed;
+				ofs.write((char *)&ntmp, sizeof(unsigned long));
+				ofs.write((char *)&stmp, sizeof(unsigned long));
+				ofs.write((char *)&mass[0], n*sizeof(double));
+				ofs.write((char *)& pos[0], n*sizeof(dvec3));
+				ofs.write((char *)& vel[0], n*sizeof(dvec3));
+				if(!ofs.fail()){
+					fprintf(stdout, "plummer : wrote to %s\n", filename);
+				}
+			}
+		}
+	}
+	static double my_rand(void) {
+		return rand()/(1. + RAND_MAX);
+	}
+};
diff --git a/examples_ptx/nbody/vector3.h b/examples_ptx/nbody/vector3.h
new file mode 100644
index 00000000..7d772023
--- /dev/null
+++ b/examples_ptx/nbody/vector3.h
@@ -0,0 +1,135 @@
+#pragma once
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+template <class REAL> struct vector3{
+  public:
+    REAL x, y, z;
+
+    vector3(){
+      x = y = z = REAL(0);
+    }
+    vector3(const REAL &r){
+      x = y = z = r;
+    }
+    vector3(const REAL &_x, const REAL &_y, const REAL &_z){
+      x = _x;  y = _y;  z = _z;
+    }
+    vector3(const REAL *p){
+      x = p[0]; y = p[1]; z = p[2];
+    }
+    ~vector3(){}
+
+    REAL &operator [](int i){
+      return (&x)[i];
+    }
+    const REAL &operator [](int i) const{
+      return (&x)[i];
+    }
+    template <class real> 
+      operator vector3<real> () const {
+        return vector3<real> (real(x), real(y), real(z));
+      }
+    operator REAL *(){
+      return &x;
+    }
+    REAL (*toPointer())[3]{
+      return (REAL (*)[3])&x;
+    }
+    typedef REAL (*pArrayOfReal3)[3];
+    operator pArrayOfReal3(){
+      return toPointer();
+    }
+
+    void outv(std::ostream &ofs = std::cout) const{
+      ofs << "(" << x << ", " << y << ", " << z << ")" << std::endl;
+    }
+    bool are_numbers () const{
+      // returns false if *this has (a) NaN member(s)
+      return (norm2() >= REAL(0));
+    }
+
+    REAL norm2() const{
+      return (*this)*(*this);
+    }
+    REAL abs() const{
+      return std::sqrt(norm2());
+    }
+
+    friend std::ostream &operator << (std::ostream &ofs, const vector3<REAL> &v){
+      ofs << v.x << " " << v.y << " " << v.z;
+      return ofs;
+    }
+    friend std::istream &operator >> (std::istream &ifs, vector3<REAL> &v){
+      ifs >> v.x >> v.y >> v.z;
+      return ifs;
+    }
+    const vector3<REAL> operator + (const vector3<REAL> &v) const{
+      return vector3<REAL> (x+v.x, y+v.y, z+v.z);
+    }
+    const inline vector3<REAL> operator - (const vector3<REAL> &v) const{
+      return vector3<REAL> (x-v.x, y-v.y, z-v.z);
+    }
+    const vector3<REAL> operator * (const REAL &s) const{
+      return vector3<REAL> (x*s, y*s, z*s);
+    }
+    friend const vector3<REAL> operator * (const REAL &s, const vector3<REAL> &v){
+      return v*s;
+    }
+    // dot product
+    const inline REAL operator * (const vector3<REAL> &v) const{
+      return (x*v.x + y*v.y + z*v.z);
+    }
+    // vector product
+    const vector3<REAL> operator % (const vector3<REAL> &v) const{
+      return vector3<REAL> (
+          y*v.z - z*v.y, 
+          z*v.x - x*v.z, 
+          x*v.y - y*v.x);
+    }
+    const vector3<REAL> operator / (const REAL &s) const{
+      REAL r = REAL(1)/s;
+      return (*this)*r;
+    }
+    const vector3<REAL> operator = (const vector3<REAL> &v){
+      x = v.x; y=v.y; z=v.z;
+      return *this;
+    }
+
+    const vector3<REAL> operator - (){
+      return vector3<REAL> (-x, -y, -z);
+    }
+    const vector3<REAL> &operator += (const vector3<REAL> &v){
+      *this = *this + v;
+      return *this;
+    }
+    const vector3<REAL> &operator -= (const vector3<REAL> &v){
+      *this = *this - v;
+      return *this;
+    }
+    const vector3<REAL> &operator *= (const REAL &s){
+      *this = *this * s;
+      return *this;
+    }
+    const vector3<REAL> &operator /= (const REAL &s){
+      *this = *this / s;
+      return *this;
+    }
+
+    friend const vector3<REAL> maxeach (const vector3<REAL> &a, const vector3<REAL> &b){
+      return vector3<REAL> (std::max(a.x, b.x), std::max(a.y, b.y), std::max(a.z, b.z));
+    }
+    friend const vector3<REAL> mineach (const vector3<REAL> &a, const vector3<REAL> &b){
+      return vector3<REAL> (std::min(a.x, b.x), std::min(a.y, b.y), std::min(a.z, b.z));
+    }
+    const vector3<REAL> abseach(){
+      return vector3<REAL> (std::fabs(x), std::fabs(y), std::fabs(z));
+    }
+};
+
+typedef vector3<double> dvec3;
+typedef vector3<float>  fvec3;
+
+