// #######################################################
//   星数密度分布図の作成 (解析プログラム)
// #######################################################

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

int main (void){
  int nmax = 1181214, lmax = 301, mmax = 241;
  static float x[1181214], y[1181214], Vmag[1181214],
    u[1181214], v[1181214], w[1181214],
    density[241][301];
  float R0 = 0.05, x1 = 75.0, dx = -0.05, y1 = 20.0, dy = 0.05;
  int n, m, l, count;
  float x_grid, y_grid, u_grid, v_grid, w_grid, cosR0, cosy0, siny0;
  FILE *fp;

// #######################################################
//   [1] データの読みこみ
// #######################################################

  fp = fopen ("gaia.txt", "r"); // 入力ファイルを開く。
  for (n=1; n<=nmax; n++){
    fscanf (fp, "%f %f %f", &x[n-1], &y[n-1], &Vmag[n-1]);
                                // データを読みこむ。
  }
  fclose (fp);                  // 入力ファイルを閉じる。

// #######################################################
//   [2] 解析
// #######################################################

  cosR0 = cos(R0*M_PI/180.0);
  for (n=1; n<=nmax; n++){
    u[n-1] = cos(y[n-1]*M_PI/180.0)*cos(x[n-1]*M_PI/180.);
    v[n-1] = cos(y[n-1]*M_PI/180.0)*sin(x[n-1]*M_PI/180.);
    w[n-1] = sin(y[n-1]*M_PI/180.0);
  }
  for (m=0; m<=mmax-1; m++){
    y_grid = y1+m*dy;
    cosy0 = cos(y_grid*M_PI/180.0);
    siny0 = sin(y_grid*M_PI/180.0);
    for (l=0; l<=lmax-1; l++){
      x_grid = x1+l*dx;
      u_grid = cosy0*cos(x_grid*M_PI/180.0);
      v_grid = cosy0*sin(x_grid*M_PI/180.0);
      w_grid = siny0;
      count = 0;
      for (n=1; n<=nmax; n++){
	if (u_grid*u[n-1] + v_grid*v[n-1] + w_grid*w[n-1] >= cosR0){
	  count++;
	}
      }
      density[m][l] = count / (M_PI*R0*R0);
    }
    printf (" %3d / %3d\n", m, mmax-1);
  }

// #######################################################
//   [3] 結果の書き出し
// #######################################################

  fp = fopen ("density2.dat", "wb"); // 出力ファイルを開く。
  fwrite (density, sizeof(float), mmax*lmax, fp);
                                    // データを書き出す。
  fclose (fp);                      // 出力ファイルを閉じる。

  return 0;
}