////////////////////////////////////////////////////////////////
//	
//	各衛星(GPS)の仰角と方位角を計算(ユーザ位置からみた)
//      
////////////////////////////////////////////////////////////////

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

extern FILE *ele,*azi;

#define DIV1     1//何個(縦横)にグリッドを区切るか
#define DIVMAX  4//グリッドを格納する配列の最大値を設定

const int SATMAX             = 36;//GPSの最大総数
const int GPS                = 32;//GPS衛星総数
const double PI       = 3.1415926535898;//円周率

extern double Sat_xpos2[DIVMAX][DIVMAX][SATMAX],Sat_ypos2[DIVMAX][DIVMAX][SATMAX];
extern double Sat_zpos2[DIVMAX][DIVMAX][SATMAX];//地平直交座標系
extern double Sat_ele[DIVMAX][DIVMAX][SATMAX],Sat_azi[DIVMAX][DIVMAX][SATMAX];//衛星の仰角、方位角
extern double GPSTIME,START_TIME;
extern int UNUSE_GPS[16];//アルマナック情報のない衛星をストックする
extern int Ref_sat,QZSS;

void calc_angle()
{
	int i,j,k;
	double root,stock;
	
	for(i=0;i<=DIV1;i++){
		for(j=0;j<=DIV1;j++){
			for(k=1;k<=GPS+QZSS;k++){

				if(k==UNUSE_GPS[0] || k==UNUSE_GPS[1] || k==UNUSE_GPS[2] || 
					k==UNUSE_GPS[3] || k==UNUSE_GPS[4] || k==UNUSE_GPS[5] || k== UNUSE_GPS[6])
					continue;
		
		root=sqrt(Sat_xpos2[i][j][k]*Sat_xpos2[i][j][k]+
			Sat_ypos2[i][j][k]*Sat_ypos2[i][j][k]);		
		Sat_ele[i][j][k] = atan(Sat_zpos2[i][j][k]/root);
		Sat_ele[i][j][k] = 360.0*Sat_ele[i][j][k]/2.0/PI;
				
		//0<azi<90
		if(Sat_xpos2[i][j][k]>=0 && Sat_ypos2[i][j][k]>=0){
		Sat_azi[i][j][k] = atan(Sat_xpos2[i][j][k]/Sat_ypos2[i][j][k]);
		Sat_azi[i][j][k] = 180.0*Sat_azi[i][j][k]/PI;
		}
		//90<azi<180
		if(Sat_xpos2[i][j][k]>=0 && Sat_ypos2[i][j][k]<=0){
		Sat_azi[i][j][k] = PI+atan(Sat_xpos2[i][j][k]/Sat_ypos2[i][j][k]);
		Sat_azi[i][j][k] = 180.0*Sat_azi[i][j][k]/PI;
		}
		//180<azi<270
		if(Sat_xpos2[i][j][k]<=0 && Sat_ypos2[i][j][k]<=0){
		Sat_azi[i][j][k] = PI+atan(Sat_xpos2[i][j][k]/Sat_ypos2[i][j][k]);
		Sat_azi[i][j][k] = 180.0*Sat_azi[i][j][k]/PI;
		}
		//270<azi<360
		if(Sat_xpos2[i][j][k]<=0 && Sat_ypos2[i][j][k]>=0){
		Sat_azi[i][j][k] = 2.0*PI+atan(Sat_xpos2[i][j][k]/Sat_ypos2[i][j][k]);
		Sat_azi[i][j][k] = 180.0*Sat_azi[i][j][k]/PI;
		}
				
			}
		}
	}

	if((int(GPSTIME)%5)==0){
		fprintf(ele,"%f ",GPSTIME);
		for(i=1;i<=GPS+QZSS;i++)
			fprintf(ele,",%f",Sat_ele[0][0][i]);
		fprintf(ele,"\n");

		fprintf(azi,"%f ",GPSTIME);
		for(i=1;i<=GPS+QZSS;i++)
			fprintf(azi,",%f",Sat_azi[0][0][i]);
		fprintf(azi,"\n");
	}

	//移動局の最大仰角衛星を求めておく
	//最初のGPS時刻で求めておかないと基準衛星が決まらない
	if(GPSTIME==START_TIME){
		for(k=1;k<=GPS+QZSS;k++){
			if(Sat_ele[DIV1][DIV1][k]>=stock){
				stock = Sat_ele[DIV1][DIV1][k];
				Ref_sat = k;
			}
		}
	}
}