// planet.pde ― 太陽系惑星軌道シミュレーション（7体問題）
// ・惑星の質量は実際の約100倍にしてある（太陽の質量を約1/100にしたのと同じ）。
//   そのため惑星どうしの引力が強く、現実の太陽系では起きない接近（木星と土星のフライバイなど）が起こる。
// ・原作は岡本（Y.Okamoto）が1991年に作った N88-BASIC のプログラム planet.bas。
//   2026年に Claude（Anthropic）が Processing に書き直した。
//
// 以下は原作の見出し
// [planet.bas] [1.00]
// calculation for planetary orbits 7-planets
//       standard model
//       sun's ,mass reduced about 1/100
//       prog. by Y.Okamoto  1991.04/29-
// 画面は原作(640×400)の2倍の1280×800。表示の数値（X00, Y00, MAG, RMAG）を2倍にしてある。

// Processing 版で加えた値：1フレーム（約1/60秒）に進める計算の回数
int nstep = 1;

// Processing 版で加えた値：文字の大きさ（ドット）。
// 文字の位置も、1文字の幅を fs/2、高さを fs として原作の LOCATE から計算するので、
// これを変えると大きさと位置がまとめて変わる
float fs = 32;

// 色（PC-98の色番号 0〜7）
color[] COL = {
  color(0, 0, 0),       // 0: 黒
  color(0, 0, 255),     // 1: 青
  color(255, 0, 0),     // 2: 赤
  color(255, 0, 255),   // 3: 紫
  color(0, 255, 0),     // 4: 緑
  color(0, 255, 255),   // 5: 水色
  color(255, 255, 0),   // 6: 黄
  color(255, 255, 255), // 7: 白
  color(255, 160, 0)    // 8: 橙（Processing 版で追加。土星用）
};

// Processing 版で加えた値：惑星を描くときの最小の半径（ドット）
float rmin = 3;

int pmax = 7;
float g = 1;
// 配列は原作と同じく 1〜pmax の番号で使う
float[] m = new float[pmax+1], x = new float[pmax+1], y = new float[pmax+1];
float[] vx = new float[pmax+1], vy = new float[pmax+1];
float[] ax = new float[pmax+1], ay = new float[pmax+1];
float[] px = new float[pmax+1], py = new float[pmax+1];
float[][] r = new float[pmax+1][pmax+1];
int[] cl = new int[pmax+1];
// 惑星の目印の文字と、その横の位置（各惑星の最初の位置から計算する。Processing 版で変更）
String[] lb = {"", "s", "m", "v", "e", "m", "j", "s"};
float[] lx = new float[pmax+1];
float ly = 256;  // 目印の縦の位置（ドット）。原作の8行目にあたる。文字の大きさ fs を変えても動かない
float x00 = 640, y00 = 400, mag = 60, rmag = 12, vmag = 0.0333;
float dt = 0.01;
float t;

void setup() {
  size(1280, 800);
  noSmooth();  // 色0で描き直して消すので、線のふちをぼかさない
  textFont(createFont("SansSerif", fs));
  textAlign(LEFT, TOP);

  // data for planets
  //   mass        radius of orbit   velocity on orbits
  m[1] = 1;       x[1] = 0;      y[1] = 0;  vx[1] = 0;  vy[1] = 0;            // sun
  m[2] = 0.00001; x[2] = 0.3871; y[2] = 0;  vx[2] = 0;  vy[2] = 47.36*vmag;   // mercury
  m[3] = 0.00024; x[3] = 0.7233; y[3] = 0;  vx[3] = 0;  vy[3] = 35.02*vmag;   // venus
  m[4] = 0.0003;  x[4] = 1;      y[4] = 0;  vx[4] = 0;  vy[4] = 29.78*vmag;   // earth
  m[5] = 0.00003; x[5] = 1.524;  y[5] = 0;  vx[5] = 0;  vy[5] = 24.08*vmag;   // mars
  m[6] = 0.09548; x[6] = 5.203;  y[6] = 0;  vx[6] = 0;  vy[6] = 13.06*vmag;   // jupiter
  m[7] = 0.02859; x[7] = 9.555;  y[7] = 0;  vx[7] = 0;  vy[7] = 9.65*vmag;    // saturn

  for (int i = 1; i <= pmax; i++) {
    cl[i] = (i % 7) + 1;  // planet color
  }
  cl[7] = 8;  // 原作では土星が色1（青）になり背景と同じで見えないので、橙に変えた

  for (int i = 1; i <= pmax; i++) {
    lx[i] = mag*x[i] + x00;  // 目印の横の位置＝惑星の最初の位置
  }

  background(COL[1]);  // 原作の LINE (0,0)-(639,399),1,BF
  drawText();

  t = 0;
  rCalc();
  aCalc();

  for (int i = 1; i <= pmax; i++) {
    vx[i] = vx[i] + ax[i]*(dt/2);
    vy[i] = vy[i] + ay[i]*(dt/2);
    x[i] = x[i] + vx[i]*dt;  px[i] = x[i]*mag;
    y[i] = y[i] + vy[i]*dt;  py[i] = y[i]*mag;
    drawCircle(mag*x[i] + x00, -mag*y[i] + y00, m[i]*rmag, COL[cl[i]]);
  }
}

// 原作の1440〜1600行のくり返しを、1フレームに nstep 回ずつ行う
void draw() {
  for (int n = 0; n < nstep; n++) {
    t = t + dt;
    rCalc();
    aCalc();
    for (int i = 1; i <= pmax; i++) {
      drawCircle(px[i] + x00, -py[i] + y00, m[i]*rmag, COL[0]);
      stroke(COL[0]);
      point(px[i] + x00, -py[i] + y00);
    }
    x00 = x00 - vx[1]*dt*mag;  y00 = y00 + vy[1]*dt*mag;
    for (int i = 1; i <= pmax; i++) {
      vx[i] = vx[i] + ax[i]*dt;
      vy[i] = vy[i] + ay[i]*dt;
      x[i] = x[i] + vx[i]*dt;  px[i] = x[i]*mag;
      y[i] = y[i] + vy[i]*dt;  py[i] = y[i]*mag;
      drawCircle(px[i] + x00, -py[i] + y00, m[i]*rmag, COL[cl[i]]);
    }
  }
  // PC-98 では文字は図形と別の画面に出るので、図形に消されない。
  // Processing では同じ画面なので、毎回書き直して同じ見え方にする
  drawText();
  // 原作の ESC キーでメニューに戻る処理は、Processing では ESC キーで終了になる
}

// 原作の CIRCLE。
// 半径が rmin に満たない惑星は、見やすいように半径 rmin の塗りつぶした円で描く（Processing 版で変更）
void drawCircle(float cx, float cy, float rad, color c) {
  stroke(c);
  strokeWeight(2);
  if (rad < rmin) {
    fill(c);
    ellipse(cx, cy, 2*rmin, 2*rmin);
  } else {
    noFill();
    ellipse(cx, cy, 2*rad, 2*rad);
  }
}

void drawText() {
  fill(COL[7]);
  text("＊太陽系惑星軌道シミュレーション(7 body problem)＊", 0*fs/2, 0*fs);
  text("７惑星相互の引力を計算して軌道上の運動を示す", 0*fs/2, 2*fs);
  text("LIST上で惑星の位置,速度等変更可能", 0*fs/2, 3*fs);
  text("参考文献：新ＢＡＳＩＣによる物理／平田邦男著／共立出版", 0*fs/2, 4*fs);
  text("（ＥＳＣキーでメニューに戻る）", 0*fs/2, 5*fs);
  // 原作は LOCATE 39,8 などの桁で置いていた。ここでは各惑星の最初の位置の真上に置く
  textAlign(CENTER, TOP);
  for (int i = 1; i <= pmax; i++) {
    text(lb[i], lx[i], ly);
  }
  textAlign(LEFT, TOP);
}

void rCalc() {
  for (int i = 1; i <= pmax; i++) {
    for (int j = 1; j <= pmax; j++) {
      if (i <= j) continue;
      r[i][j] = sqrt(pow(x[i] - x[j], 2) + pow(y[i] - y[j], 2));
      r[j][i] = r[i][j];
    }
  }
}

void aCalc() {
  for (int i = 1; i <= pmax; i++) {
    float sgmx = 0, sgmy = 0;
    for (int j = 1; j <= pmax; j++) {
      if (i == j) continue;
      sgmx = sgmx - g*m[j]*(x[i] - x[j])/pow(r[i][j], 3);
      sgmy = sgmy - g*m[j]*(y[i] - y[j])/pow(r[i][j], 3);
    }
    ax[i] = sgmx;
    ay[i] = sgmy;
  }
}