#include <stdio.h>
#include <math.h>

#define V0    20.0   // 初速
#define THETA 45.0*M_PI/180.0
                     // 打ち上げ角度 [rad]
#define d_t    0.001 // サンプリング時間

void main(void)
{
  const double m = 0.1,   // 質量
               g = 9.8,   // 重力加速度
               k = 0.01;  // 慣性抵抗係数
  double x[2],  // 位置ベクトル
    v[2],       // 速度ベクトル
    a[2],       // 加速度ベクトル
    t = 0.0;    // 現在時刻
  double speed; // 速さ
  FILE *fp;     // ファイルポインタ


  // 位置の初期化 (初期位置は原点)
  x[0] = 0.0;
  x[1] = 0.0;

  // 速度の初期化
  v[0] = V0 * cos(THETA);
  v[1] = V0 * sin(THETA);

  // 出力ファイルの初期化
  fp = fopen("out.csv", "w");
  if (fp == NULL) {
    fprintf(stderr, "File cannot open.\n");
    return;
  }

  // 初期状態をファイルへ出力
  fprintf(fp, "%f,%f,%f\n", t, x[0], x[1]);

  // 時刻ごとのループ
  for (t = 0.0; ; t += d_t) {

    // 速さを計算
    speed = sqrt(v[0]*v[0]+v[1]*v[1]);

    // 加速度ベクトル a[] の計算
    a[0] = - k / m * speed * v[0];
    a[1] = - g - k / m * speed * v[1];

    // 位置ベクトル x[] の更新
    x[0] += v[0] * d_t;
    x[1] += v[1] * d_t;

    // 速度ベクトル v[] の更新
    v[0] += a[0] * d_t;
    v[1] += a[1] * d_t;

    // ファイルへの出力
    fprintf(fp, "%f,%f,%f\n", t, x[0], x[1]);

    // 落下していたらループを抜ける
    if (x[1] < 0.0) break;
  }

  // ファイルのクローズ
  fclose(fp);

  // 飛行距離の画面表示
  printf("Flight distance: %f.\n", x[0]);
}