///////////////////////////////////////////////////////////////////////////////////////////////
///////DIMM観測結果のファイルから年月日時分秒(YMDhms)、観測開始時刻からの経過時間(t)、/////////
///////シーイング値12個を読み取り天頂角を補正した値を返すプログラム                   /////////
///////東北大学理学部宇宙地球物理学科天文学コース学部4年 沖田博文                     /////////
///////E-mail: h-okita@astr.tohoku.ac.jp                       Last Update:2008/03/13 /////////
///////                                                                               /////////
///////広大DIMMとの比較観測用に出力データを変更。出力するのは観測時間ではなく、観測時 /////////
///////刻(JST)とした。ただしUTC15時縲鰀24時では日付が1日ずれるので注意が必要。          /////////
///////                                                        last Update:2008/07/16 /////////
///////                                                                               /////////
///////JST30時間制にした。日付はこの場合ずれない(はず)         last Update:2008/10/06 /////////
///////////////////////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>
#define pi 3.14159265	//円周率

int main(void)
{

//////////以下の情報を入力/////////////////////////////////////////////////////////////////////
	char file[100]="seeing.001" ;	//処理するファイル名
	char city[100]="東北大学理学部物理A棟(宮城県仙台市)" ;	//観測地
	#define l 140.83333	//経度(東経は+,西経は-)
	#define f 38.266667	//緯度(北緯は+,南緯は-)
	#define k 0.0	//観測に用いた標準時の経度(東経は+,西経は-)(世界時の場合は0.0)
	char star[100]="カペラ(Alpha-Aur)" ;	//観測に用いた恒星名
	#define a 79.1725	//天体の赤経(°単位)
	#define d 45.998056	//天体の赤緯
	#define N 12433	//ループの回数(処理するファイルの行数)
////////////////////////////////////////////////////////////////////////////////////////////////

	int n=1, i, j ;
	float Ave=0.0, d14[N], Med ;
	float s=0.0, t=0.0 ;
	FILE *fi ;

	fi=fopen(file, "r") ;	//処理するファイル名
	if(fi == NULL){
		printf("No File. \n") ;
	}
	else{  
		printf("#Z Time(JST) Vla Vta Hla Hta Vlo Vto Hlo Hto Vle Vte Hle Hte Ave\n") ;
		printf("# ---------------------------------------------------------------------------\n") ;
		for(n=1; n<=N; n=n+1){
			int Y=0, M=0, D=0, h=0, m=0 ;

			float d1=0.0, d2=0.0, d3=0.0, d4=0.0, d5=0.0, d6=0.0, d7=0.0, d8=0.0, d9=0.0, d10=0.0, d11=0.0, d12=0.0, d13=0.0 ;
       			float ae=0.0, af=0.0, ag=0.0, ah=0.0, ai=0.0 ;
			char aa[1], ab[1], ac[1], ad[1] ;
			float ax=0.0, dx=0.0, fx=0.0 ;

			float DA=0.0 ;

			float A=0.0, B=0.0, C=0.0, E=0.0, F=0.0, G=0.0, H=0.0, z=0.0, Z=0.0, I=0.0 ;
			//ファイルの読み込み
			//        Y     M     D   h     m     s    t                                   d1   d2   d3   d4   d5   d6   d7   d8   d9   d10  d11  d12    
			fscanf(fi,"%4d%1s%2d%1s%2d %2d%1s%2d%1s%9f  %12f	%5f %5f %5f %5f %5f	%12f %12f %12f %12f %12f %12f %12f %12f %12f %12f %12f %12f \n", &Y, aa, &M, ab, &D, &h, ac, &m, ad, &s, &t, &ae, &af, &ag, &ah, &ai, &d1, &d2, &d3, &d4, &d5, &d6, &d7, &d8, &d9, &d10, &d11, &d12 ) ;

			//読み込み結果の表示
			  //      Y   M   D   h   m   s   t    d1   d2   d3   d4   d5   d6   d7   d8   d9   d10  d11  d12
			//printf("%4d %2d %2d %2d %2d %9f %14f %12f %12f %12f %12f %12f %12f %12f %12f %12f %12f %12f %12f \n", Y, M, D, h, m, s, t, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12) ;

			//20080716追加
			//時刻を日本標準時に変換
			DA = h + m/60.0 + s/3600.0 + 9 ;

			if(DA < 24){
				DA=DA ;
			}
			else{
				DA=DA-24.0 ;
			}

			//恒星時Gの計算
				//経過日数Aの計算
			if(M > 2){
				A=365*(Y-1900) + 30*M + D - 33.5 + floor(3*(M+1)/5.0) + floor((Y-1900)/4.0) ;
			}
			else{
				A=365*(Y-1901) + 30*(M+12) + D - 33.5 + floor(3*(M+13)/5.0) + floor((Y-1901)/4.0) ;
			}
			//経過日数Tuの値Bを求める
			B=A/36525.0 ;
			//観測日の世界時0時のグリニッジ視恒星時θ0の値Dを求める
			C=(6*3600+38*60+45.836)+8640184.542*B+0.0929*B*B ;
			E=(C/(3600*24)-floor(C/3600.0/24.0))*24*15 ;	
			//ある観測時刻、観測地での恒星時θの値Gを求める
			F=E + l + ((h + m/60.0 +s/3600.0)*15 -k)*1.00273791 ;
			G=F - floor(F/360.0)*360 ;
			//計算結果の表示
			//printf("%f %f %f %f %f %f\n", A, B, C, E, F, G) ;

			//degree単位からRadian単位に変換
			ax=a*pi /180.0 ;
			dx=d*pi /180.0 ;
			fx=f*pi /180.0 ;
			G=G*pi /180.0 ;

			//この観測地点、観測時刻での(緯度f,恒星時Gでの)赤経a,赤緯dの天体の高度H、天頂角zを求める
			H=asin(sin(fx)*sin(dx)+cos(fx)*cos(dx)*cos(G-ax)) ;
			z=pi/2.0-H ;
			//計算結果の表示
			//printf("%f %f\n", H, z) ;

			//DIMMでの天頂角補正の値Z	
			Z=pow(cos(z), 3/5.0) ;
			//計算結果の表示
			//printf("%f\n", Z) ;

			//観測結果に天頂角補正を行う
			d1=d1*Z ;
			d2=d2*Z ;
			d3=d3*Z ;
			d4=d4*Z ;
			d5=d5*Z ;
			d6=d6*Z ;
			d7=d7*Z ;
			d8=d8*Z ;
			d9=d9*Z ;
			d10=d10*Z ;
			d11=d11*Z ;
			d12=d12*Z ;

			//d5縲彭12の平均値
			d13= (d5+d6+d7+d8+d9+d10+d11+d12)/8.0 ;
			//シーイングの時間平均を求める計算
			Ave=Ave+d13 ;

			//30時間制へ変換
			if(DA <= 6){
				DA=DA+24 ;
			}

			//最終的な観測結果の表示
			//      Z   DA   d1  d2  d3  d4  d5  d6  d7  d8  d9  d10 d11 d12 d13
			printf("%8f %8f %8f %8f %8f %8f %8f %8f %8f %8f %8f %8f %8f %8f %8f \n", Z, DA, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12, d13) ;
			//Median用のデータ取得
			d14[n]=d13 ;
		}
		// Medianデータを昇順に並び替え
		for(i = 0 ; i < N-1 ; i = i+1){
			for(j = N-1 ; j > i ; j = j-1){
				if(d14[j] < d14[j-1] ){
					float tmp = d14[j];
					d14[j] = d14[j-1];
					d14[j-1] = tmp;
				}
			}
		}
		// Medianの計算
		if(N %2 == 1){
			Med = d14[N/2] ;
		}
		else{
			Med = (d14[N/2-1] + d14[N/2]) / 2.0 ;
		}
		//平均値の計算
		Ave=Ave/N ;
		fclose(fi) ;

		printf("# ---------------------------------------------------------------------------\n") ;
		printf("# 解析したファイル: %s\n", file) ;
		printf("# ---------------------------------------------------------------------------\n") ;
                printf("# 観測地: %s\n", city) ;
		printf("# 観測値の経度: %f°  緯度: %f°\n", l, f) ;	
                printf("# 観測に用いた恒星名: %s\n", star) ;
		printf("# 観測に用いた恒星の赤経: %f°  赤緯: %f°\n", a, d) ;
		fi=fopen(file, "r") ;  
		int Y1=0, M1=0, D1=0, h1=0, m1=0 ;
		float s1=0.0 ;
		char aa1[1], ab1[1], ac1[1], ad1[1] ;
		fscanf(fi,"%4d%1s%2d%1s%2d %2d%1s%2d%1s%9f\n", &Y1, aa1, &M1, ab1, &D1, &h1, ac1, &m1, ad1, &s1 ) ;
		printf("# 観測日: %4d年%2d月%2d日\n", Y1, M1, D1) ;
		printf("# 観測開始時刻: %2d時%2d分%9f秒(世界時)\n",  h1, m1, s1) ;
		printf("# 観測時間: %f秒\n",t) ;
		fclose(fi) ;
		printf("# ---------------------------------------------------------------------------\n") ;
		printf("# 観測結果\n") ;
		printf("# シーイング(最大値): %f\n", d14[N-1]) ; 
		printf("# シーイング(最小値): %f\n", d14[1]) ;
		printf("# シーイング(平均値): %f\n", Ave) ;
		printf("# シーイング(Median): %f\n", Med) ;
		printf("# ---------------------------------------------------------------------------\n") ;
		printf("# データ解説\n") ;
		printf("#  1列目   Z:補正値\n") ;
		printf("#  2列目   H:観測時刻(日本標準時)\n") ;
		printf("#  3列目 Vla:偶数行画像、奇数行画像の平均の縦方向のペアのlongitudinalシーイング\n") ;
		printf("#  4列目 Vta:偶数行画像、奇数行画像の平均の縦方向のペアのtransverseシーイング\n") ;
		printf("#  5列目 Hla:偶数行画像、奇数行画像の平均の縦方向のペアのlongitudinalシーイング\n") ;
		printf("#  6列目 Hta:偶数行画像、奇数行画像の平均の縦方向のペアのtransverseシーイング\n") ;
		printf("#  7列目 Vlo:奇数行画像の縦方向のペアのlongitudinalシーイング\n") ;
		printf("#  8列目 Vto:奇数行画像の縦方向のペアのtransverseシーイング\n") ;
		printf("#  9列目 Hlo:奇数行画像の縦方向のペアのlongitudinalシーイング\n") ;
		printf("# 10列目 Hto:奇数行画像の縦方向のペアのtransverseシーイング\n") ;
		printf("# 11列目 Vle:偶数行画像の縦方向のペアのlongitudinalシーイング\n") ;
		printf("# 12列目 Vte:偶数行画像の縦方向のペアのtransverseシーイング\n") ;
		printf("# 13列目 Hle:偶数行画像の縦方向のペアのlongitudinalシーイング\n") ;
		printf("# 14列目 Hte:偶数行画像の縦方向のペアのtransverseシーイング\n") ;
		printf("# 15列目 Ave:7列から14列のシーイングの平均\n") ;
		printf("# ---------------------------------------------------------------------------\n") ;

	}
	return 0 ;
}