float m_size;
float m=0.2;
float h=0.01;
float[][] coord = new float[400][2];
int k = 0;
float[] r0 = {  2.00, -0.01};
float[] rh = {  2.00, 0.00};

float x0,x1,y0,y1,z0,z1;
float dx,dy,dz,dc;
float nowX, nowY, oldX, oldY;

void setup() {
  frameRate(30);
  size( 500, 300, P3D);
  coord[0]= rh;
  coord[1]= r0; 
  for (int i=0; i < 391; i++) {
    coord[i+2] = verlet(coord[i+1], coord[i]);
  }
  nowX=nowY=oldX=oldY=0.0;
  set_mesh( -1, 2.5, -1.5, 1.5, -2.5, 0, 200);
}

float[] force(float pos[]) {
  float x, y, L;
  float[] posforce = new float[2];
  x=pos[0];
  y=pos[1];
  L=pow( (x*x+y*y), 1.5);
  posforce[0] = -x/L;
  posforce[1] = -y/L;
  return posforce;
}

//現在と過去の座標から次の座標
float[] verlet(float now[], float past[]) {
  float[] posforce = new float[2];
  float[] newpos = new float[2];
  posforce = force(now);
  newpos[0]=2*now[0]-past[0]+h*h/m*posforce[0];
  newpos[1]=2*now[1]-past[1]+h*h/m*posforce[1];
  return newpos;
} 

float func(float x, float y) {
  return -1.0/sqrt(x*x+y*y);
  // return -1.0*x+2.0*y-3.0;
}
float x_size,y_size,z_size;
int x_min,x_max,y_min,y_max,z_min,z_max;
void set_mesh(float tx0, float tx1, float ty0, float ty1, float tz0, float tz1, int tsize) {
  x_size=tsize;
  m_size=tsize;
  x0=tx0;x1=tx1;
  y0=ty0;y1=ty1;
  z0=tz0;z1=tz1;
  y_size=(y1-y0)/(x1-x0)*x_size;
  z_size=(z1-z0)/(x1-x0)*x_size;
  dx=(x1-x0)/x_size;
  x_min=int(x0/dx);
  x_max=int(x1/dx);
  dy=(y1-y0)/y_size;
  y_min=int(y0/dy);
  y_max=int(y1/dy);
  dz=(z1-z0)/z_size;
  z_min=int(z0/dz);
  z_max=int(z1/dz);
  dc=128/m_size;
}

void plot_point(int i, int j){
  float z = func( i*dx, j*dy )/dz;
  if (z >= z_min && z <= z_max) {
    stroke(255, 127+dc*i, 127+dc*j);
    strokeWeight(2);
    point(i, j, z); // z0/dz+(z-z0)/dz
  }
}
void plot_axez(){ //default plot
  plot_axez(0.0);
}
void plot_axez(float z){  //z<>0.0
  stroke(255, 255, 255);
  strokeWeight(2);
  line(x_min,0,z/dz, x_max,0,z/dz);
  line(0,y_min,z/dz, 0,y_max,z/dz);
  line(0,0,z0/dz, 0,0,z1/dz);
}

void plot_mesh(){
  for (int i=x_min; i<=x_max; i+=4) {
    for (int j=0; j<=y_max; j+=20) {
      plot_point(i,j);
    }
    for (int j=-20; j>=y_min; j-=20) {
      plot_point(i,j);
    }
  }
  for (int j=y_min; j<=y_max; j+=4) {
    for (int i=0; i<=x_max; i+=20) {
      plot_point(i,j);
    }
    for (int i=-20; i>=x_min; i-=20) {
      plot_point(i,j);
    }
  }
}

float[] get_xyz(float x, float y, float z){
  float [] res = {x/dx,y/dy,z/dz};
  return res;
}

void draw() {
  float mvX=0.0, mvY=0.0;
  background(10,10,10);
  camera(0, 0, m_size, 
         0, 0, -m_size/2, 
         0, 1, 0);

  if (mousePressed) {
    mvX =mouseX-oldX;
    mvY =mouseY-oldY;
  }
  oldX=mouseX;
  oldY=mouseY;
  nowX += mvX;
  nowY -= mvY;
  rotateX(radians(nowY/3));
  rotateY(radians(18.0/3));
  rotateZ(radians(-nowX/3));

  // mesh drawing
  plot_mesh();
  plot_axez(-1.25);
  // center sphere
  float [] pos=new float[3];
  pos=get_xyz(0.0,0.0,-1.25);
  translate(pos[0],pos[1],pos[2]);
  sphereDetail(12, 12);
  stroke(128, 0, 0);
  fill(128, 0, 0);
  sphere(15);
  translate(-pos[0],-pos[1],-pos[2]);

  //球を描画する座標　x、yは　座標*（単位座標のピクセル数）+　原点の座標
  // orbitar 
  pos=get_xyz(coord[k][0],coord[k][1],func(coord[k][0],coord[k][1]));
  translate(pos[0],pos[1],pos[2]);

  stroke(128, 128, 255);
  fill(128, 128, 255);
  sphere(5);

  k += 10;  // for loop
  if (k >392) k = 0;
}