/*  "<?>"の部分を自分で埋めてプログラムを完成してください。  */

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

int main(void)
{

/*  定数piとreを宣言する。
    定数piは円周率を、reは地球半径を表す。  */
  float pi=3.1415927, re=6368.*1000.;

/*  定数mmaxとlmaxを宣言する。
    定数mmaxは緯度方向の格子数を、lmaxは経度方向の格子数を表す。  */
  int mmax=73, lmax=144;

/*  配列u,v,xiを宣言する。  */
  float u[mmax][lmax], v[mmax][lmax], xi[mmax][lmax];

/*  変数を宣言する。  */
  int m, l, mm, mp, lm, lp;
  float dla, dph, ph, phm, php, um, up, vm, vp, usum, uave;

/*  ファイルポインタfpを宣言する。  */
  FILE *fp;

/*  データを読みこむ。  */

/*  東西風データのファイルを開く。
    モードは"rb"(バイナリファイルの読みこみ)を指定する。  */
  fp = fopen( "U500_201104.dat", "rb" );

/*  ファイルポインタfpで指定されたファイルの所定の位置に移動する。
    第3引数にSEEK_SETを指定し、ファイルの先頭を位置の基準とする。
    データファイルには4月1日0時UTCから6時間ごとのデータが記録されていて、
    ここでは、93番目のレコード(24日0時UTCのデータ)を読みこむので、
    まず、ファイルの先頭から
    (93-1)*mmax*lmax*sizeof(float)バイトだけ前方に移動する。  */
  fseek( fp, (93-1)*mmax*lmax*sizeof(float), SEEK_SET );

/*  配列aの各要素の値をファイルポインタfpで指定されたファイルから読みこむ。
    項目のサイズは浮動小数点のサイズsizeof(float)を指定する。
    項目の数はmmax*lmaxを指定する。  */
  fread( u, sizeof(float), mmax*lmax, fp );

/*  ファイルを閉じる。  */
  fclose( fp );

/*  南北風データのファイルを開く。
    モードは"rb"(バイナリファイルの読みこみ)を指定する。  */
  fp = fopen( "V500_201104.dat", "rb" );

/*  ファイルポインタfpで指定されたファイルの所定の位置に移動する。
    第3引数にSEEK_SETを指定し、ファイルの先頭を位置の基準とする。
    データファイルには4月1日0時UTCから6時間ごとのデータが記録されていて、
    ここでは、93番目のレコード(24日0時UTCのデータ)を読みこむので、
    まず、ファイルの先頭から
    (93-1)*mmax*lmax*sizeof(float)バイトだけ前方に移動する。  */
  fseek( fp, (93-1)*mmax*lmax*sizeof(float), SEEK_SET );

/*  配列aの各要素の値をファイルポインタfpで指定されたファイルから読みこむ。
    項目のサイズは浮動小数点のサイズsizeof(float)を指定する。
    項目の数はmmax*lmaxを指定する。  */
  fread( v, sizeof(float), mmax*lmax, fp );

/*  ファイルを閉じる。  */
  fclose( fp );


/*  ==== ここから相対渦度の計算 ====  */

/*  相対渦度を計算する。  */

/*  東西、南北格子間隔の値を計算する。
    東西間隔をdla、南北格子間隔をdphとする。
    単位はラジアンである。  */
  dla = 2.*pi/lmax;
  dph = pi/(mmax-1);

/*  ここから相対渦度の計算のためのforループが始まる。
    北極と南極を除いて、lmax*(mmax-1)回だけ反復する。  */
  for (m=1; m<=mmax-2; m++)
  {

/*  南北格子番号mに対応する緯度を計算する。
    その格子点の緯度をph、ひとつ北の格子点の緯度をphp、
    ひとつ南の格子点の緯度をphmとする。
    m=0が北極、m=mmax-1が南極である点に注意する。  */
    ph = 0.5*pi - pi*m/(mmax-1);
    php = 0.5*pi - pi*(m-1)/(mmax-1);
    phm = 0.5*pi - pi*(m+1)/(mmax-1);

    for (l=0; l<=lmax-1; l++)
    {

/*  となりの格子における東西格子番号を計算する。
    ひとつ東の格子における東西格子番号をlp、
    ひとつ西の格子における東西格子番号をlmとする。  */
      lp = l+1;
      if (l == lmax-1)
      {
	lp = 0;
      }
      lm = l-1;
      if (l == 0)
      {
	lm = lmax-1;
      }

/*  ひとつ東(LP)と西(LM)の格子における南北風Vの値を求める。
    ひとつ東の格子(LP,M)における南北風Vの値をVP、
    ひとつ西の格子(LM,M)における南北風Vの値をVMとする。  */
      vp = v[m][lp];
      vm = v[m][lm];

/*  ひとつ北(M-1)と南(M+1)の格子における東西風Uの値を求める。
    ひとつ北の格子(L,M-1)における東西風Uの値をUP、
    ひとつ南の格子(L,M+1)における東西風Uの値をUMとする。  */
      up = u[m-1][l];
      um = u[m+1][l];

/*  相対渦度を計算する。  */
      xi[m][l] = (vp-vm) / (2.*dla) / (re*cos(ph))
                 - (up*cos(php)-um*cos(phm)) / (2.*dph) / (re*cos(ph)) ;

/*  ここでforループが終了する。  */
    }
  }

/*  北極における相対渦度を計算する。  */
  usum = 0.;
  for (l=0; l<=lmax-1; l++)
  {
    usum = usum + u[1][l];
  }
  uave = usum / lmax; 
  for (l=0; l<=lmax-1; l++)
  {
    xi[0][l] = 2. * uave / (re * dph);
  }

/*  南極における相対渦度を計算する。  */
  usum = 0.;
  for (l=0; l<=lmax-1; l++)
  {
    usum = usum + u[mmax-2][l];
  }
  uave = usum / lmax; 
  for (l=0; l<=lmax-1; l++)
  {
    xi[mmax-1][l] = 2. * uave / (re * dph);
  }

/*  ==== ここまで相対渦度の計算 ====  */

/*  データを書き出す。  */

/*  ファイルを開く。
    モードは"wb"(バイナリファイルへの書き出し)を指定する。  */
  fp = fopen( "output.dat", "wb" );

/*  配列xiの各要素の値をファイルポインタfpで指定されたファイルに書き出す。
    項目のサイズは浮動小数点のサイズsizeof(float)を指定する。
    項目の数はmmax*lmaxを指定する。  */
  fwrite( xi, sizeof(float), mmax*lmax, fp );

/*  ファイルを閉じる。  */
  fclose( fp );

  return 0;
}