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

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

int main (void){
  int nmax = 1357122, lmax = 360, mmax = 181;
  static float x[1357122], y[1357122], Vmag[1357122],
    u[1357122], v[1357122], w[1357122],
    density[181][360];
  float R0 = 1.0;
  int n, m, l, count;
  float x_grid, y_grid, u_grid, v_grid, w_grid, cosR0, cosy0, siny0;
  FILE *fp;

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

  fp = fopen ("datag.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 = m*180.0/(mmax-1)-90.0;
    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 = 360.0-l*360.0/lmax;
      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 ("densityg.dat", "wb"); // 出力ファイルを開く。
  fwrite (density, sizeof(float), mmax*lmax, fp);
                                     // データを書き出す。
  fclose (fp);                       // 出力ファイルを閉じる。

  return 0;
}