////////////////////////////////////////////////////////////////
//
//	人工的にコードマルチパス誤差を生成
//　（地面からの反射分を計算）
////////////////////////////////////////////////////////////////

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

const double PI       = 3.1415926535898;//円周率
const double CS       = 299792458.0;//光速
const double F1		  = 1575420000.0;
const double F2       = 1227600000.0;

const int SATMAX             = 36;//GPSの最大総数
const int GPS                = 32;//GPSの最大総数

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

extern FILE *dd_code,*temp;

extern int Enable[DIVMAX][DIVMAX],Enable_sat[DIVMAX][DIVMAX][SATMAX];//使用可能衛星
extern int Actual_satn;//実際のＧＰＳ，ＧＥＯ衛星の運用数を入れる
extern int Enable_sat_stock[DIVMAX][DIVMAX][SATMAX];//2000/11/04
//予備
extern int Enable2[DIVMAX][DIVMAX];
extern int Enable_sat2[DIVMAX][DIVMAX][SATMAX];//使用可能衛星
extern int Total_sat;//トータル衛星数
extern int Ref_sat;//最大仰角衛星のこと
extern int UNUSE_GPS[16];//アルマナック情報のない衛星をストックする
extern int Com_sat[SATMAX],Com_satn;
extern int Integer_r1[SATMAX],Integer_m1[SATMAX],DD_integer1[SATMAX];
extern int QZSS;

extern double Pr1ref[SATMAX],Pr1mov[SATMAX],Pr2ref[SATMAX],Pr2mov[SATMAX];
extern double DD_code1[SATMAX],DD_code2[SATMAX];
extern double Sat_xpos[SATMAX],Sat_ypos[SATMAX],Sat_zpos[SATMAX];//地球中心座標系
extern double POSx[DIVMAX][DIVMAX],POSy[DIVMAX][DIVMAX],POSz[DIVMAX][DIVMAX];//解析位置
extern double Cn[SATMAX],GPSTIME;
extern double Sat_ele[DIVMAX][DIVMAX][SATMAX],Sat_azi[DIVMAX][DIVMAX][SATMAX];//衛星の仰角、方位角
extern double MP_rov[SATMAX],MPG_rov[SATMAX],MPG_ref[SATMAX];
extern double Ant_height_ref,Ant_height_rover;

double gauss_gen(double,double,int *);
double random(int *);

void ground_multipath_l1()
{
	int i;
	int sat;

	double dl,dm[SATMAX]={0};
	double amp,amph,ampv,permt;
	double repermt=7;
	double cond=0.04;
	double lm=CS/F1;
	double ph;
	double dis = 30.0;//アンテナと壁との距離
	double corr=0.1;
	double chip=CS/1023000;
	double td=corr*chip/2.0;
	double defuse1=2.5;

//	static int kaisu=0;

	//反射係数の計算
	permt = repermt-60.0*CS/F1*cond;

//	kaisu++;

	//基準局のマルチパス誤差生成
	for(i=0;i<Enable[DIV1][DIV1];i++){
		sat = Enable_sat[DIV1][DIV1][i];
	}

	//基準局のマルチパス誤差生成
	//地面からの反射分
	for(i=0;i<Enable[DIV1][DIV1];i++){
		sat = Enable_sat[DIV1][DIV1][i];
		MPG_ref[sat]=0;

		if(Sat_azi[DIV1][DIV1][sat]>=0 && Sat_azi[DIV1][DIV1][sat]<=360){			
			dl = 2.0*Ant_height_ref*cos(PI*(90.0-Sat_ele[DIV1][DIV1][sat])/180);		

			amph = cos(PI*Sat_ele[DIV1][DIV1][sat]/180) - sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180));
			amph = amph/(cos(PI*Sat_ele[DIV1][DIV1][sat]/180) + sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180)));
			amph = fabs(amph);

			ampv = permt*cos(PI*Sat_ele[DIV1][DIV1][sat]/180) - sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180));
			ampv = ampv/(permt*cos(PI*Sat_ele[DIV1][DIV1][sat]/180) + sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180)));
			ampv = fabs(ampv);

			amp = (amph+ampv)/2.0;

			//アンテナの下面から入射されるので、アンテナによる減衰効果があり
			//60-90 30dB以上減衰 0.0316
			//30 20dB減衰 0.1
			//0 13dB減衰 0.224

			double keisu;
	//		dec = 0.23*Sat_ele[0][0][sat]*2;//GPS600
	//		keisu = 1/(pow(defuse1,dec/20));
			if(Sat_ele[0][0][sat]>=60)
				keisu=0.03;
			if(Sat_ele[0][0][sat]>=30 && Sat_ele[0][0][sat]<=60)
				keisu=-0.07/30*Sat_ele[0][0][sat]+0.17;
			if(Sat_ele[0][0][sat]>=0 && Sat_ele[0][0][sat]<=30)
				keisu=-0.124/30*Sat_ele[0][0][sat]+0.224;

			amp = amp*keisu;

			//遅延距離による場合分けが必要
			if(dl>=0 && dl<(1.0-amp)*td){
				ph = dl/lm-int(dl/lm);
				ph = 2*PI*ph;
				dm[sat] = dl*amp*cos(ph)/(1+amp*cos(ph));
				MPG_ref[sat] = dm[sat];
			}
			else if(dl>=(1.0-amp)*td && dl<(1.0+amp)*td){
				ph = dl/lm-int(dl/lm);
				//同相側の場合
				if(ph<=0.25 || ph>=0.75){
					ph = 2*PI*ph;
					dm[sat] = dl*amp*cos(ph)/(1+amp*cos(ph));
					MPG_ref[sat] = dm[sat];
				}
				else{//逆相側の場合
					ph = 2*PI*ph;
					dm[sat] = amp*cos(ph)*td;
					MPG_ref[sat] = dm[sat];
				}

			}
			else if(dl>=(1.0+amp)*td && dl<chip-(1.0+amp)*td){
				ph = dl/lm-int(dl/lm);
				ph = 2*PI*ph;
				dm[sat] = amp*cos(ph)*td;
				MPG_ref[sat] = dm[sat];
			}
			else if(dl>=chip-(1.0+amp)*td && dl<chip-(1.0-amp)*td){
				ph = dl/lm-int(dl/lm);
				//同相側の場合
				if(ph<=0.25 || ph>=0.75){
					ph = 2*PI*ph;
					dm[sat] = amp*cos(ph)*td;
					MPG_ref[sat] = dm[sat];
				}
				else{//逆相側の場合
					ph = 2*PI*ph;
					dm[sat] = amp*(chip+td-dl)*cos(ph);
					dm[sat] = dm[sat]/(2-amp*cos(ph));
					MPG_ref[sat] = dm[sat];
				}
			}
			else if(dl>=chip-(1.0-amp)*td && dl<chip+td){
				ph = dl/lm-int(dl/lm);
				ph = 2*PI*ph;
				dm[sat] = amp*(chip+td-dl)*cos(ph);
				dm[sat] = dm[sat]/(2-amp*cos(ph));
				MPG_ref[sat] = dm[sat];
			}
			else{
				dm[sat]=0;
				MPG_ref[sat] = dm[sat];
			}
//			MPG_rov[sat] = 0;
			MPG_ref[sat] = MPG_ref[sat];
		}//if

	}//for



	//移動局のマルチパス誤差生成
	//地面からの反射分
	for(i=0;i<Enable[DIV1][DIV1];i++){
		sat = Enable_sat[DIV1][DIV1][i];
		MPG_rov[sat]=0;

		if(Sat_azi[DIV1][DIV1][sat]>=0 && Sat_azi[DIV1][DIV1][sat]<=360){			
			dl = 2.0*Ant_height_rover*cos(PI*(90.0-Sat_ele[DIV1][DIV1][sat])/180);		

			amph = cos(PI*Sat_ele[DIV1][DIV1][sat]/180) - sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180));
			amph = amph/(cos(PI*Sat_ele[DIV1][DIV1][sat]/180) + sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180)));
			amph = fabs(amph);

			ampv = permt*cos(PI*Sat_ele[DIV1][DIV1][sat]/180) - sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180));
			ampv = ampv/(permt*cos(PI*Sat_ele[DIV1][DIV1][sat]/180) + sqrt(permt-sin(PI*Sat_ele[DIV1][DIV1][sat]/180)*sin(PI*Sat_ele[DIV1][DIV1][sat]/180)));
			ampv = fabs(ampv);

			amp = (amph+ampv)/2.0;

			//アンテナの下面から入射されるので、アンテナによる減衰効果があり
			//60-90 30dB以上減衰 0.0316
			//30 20dB減衰 0.1
			//0 13dB減衰 0.224

			double keisu;
	//		dec = 0.23*Sat_ele[0][0][sat]*2;//GPS600
	//		keisu = 1/(pow(defuse1,dec/20));
			if(Sat_ele[0][0][sat]>=60)
				keisu=0.03;
			if(Sat_ele[0][0][sat]>=30 && Sat_ele[0][0][sat]<=60)
				keisu=-0.07/30*Sat_ele[0][0][sat]+0.17;
			if(Sat_ele[0][0][sat]>=0 && Sat_ele[0][0][sat]<=30)
				keisu=-0.124/30*Sat_ele[0][0][sat]+0.224;

			amp = amp*keisu;

			//遅延距離による場合分けが必要
			if(dl>=0 && dl<(1.0-amp)*td){
				ph = dl/lm-int(dl/lm);
				ph = 2*PI*ph;
				dm[sat] = dl*amp*cos(ph)/(1+amp*cos(ph));
				MPG_rov[sat] = dm[sat];
			}
			else if(dl>=(1.0-amp)*td && dl<(1.0+amp)*td){
				ph = dl/lm-int(dl/lm);
				//同相側の場合
				if(ph<=0.25 || ph>=0.75){
					ph = 2*PI*ph;
					dm[sat] = dl*amp*cos(ph)/(1+amp*cos(ph));
					MPG_rov[sat] = dm[sat];
				}
				else{//逆相側の場合
					ph = 2*PI*ph;
					dm[sat] = amp*cos(ph)*td;
					MPG_rov[sat] = dm[sat];
				}

			}
			else if(dl>=(1.0+amp)*td && dl<chip-(1.0+amp)*td){
				ph = dl/lm-int(dl/lm);
				ph = 2*PI*ph;
				dm[sat] = amp*cos(ph)*td;
				MPG_rov[sat] = dm[sat];
			}
			else if(dl>=chip-(1.0+amp)*td && dl<chip-(1.0-amp)*td){
				ph = dl/lm-int(dl/lm);
				//同相側の場合
				if(ph<=0.25 || ph>=0.75){
					ph = 2*PI*ph;
					dm[sat] = amp*cos(ph)*td;
					MPG_rov[sat] = dm[sat];
				}
				else{//逆相側の場合
					ph = 2*PI*ph;
					dm[sat] = amp*(chip+td-dl)*cos(ph);
					dm[sat] = dm[sat]/(2-amp*cos(ph));
					MPG_rov[sat] = dm[sat];
				}
			}
			else if(dl>=chip-(1.0-amp)*td && dl<chip+td){
				ph = dl/lm-int(dl/lm);
				ph = 2*PI*ph;
				dm[sat] = amp*(chip+td-dl)*cos(ph);
				dm[sat] = dm[sat]/(2-amp*cos(ph));
				MPG_rov[sat] = dm[sat];
			}
			else{
				dm[sat]=0;
				MPG_rov[sat] = dm[sat];
			}
//			MPG_rov[sat] = 0;
			MPG_rov[sat] = MPG_rov[sat];
		}//if

	}//for

//	if((i%1)==0){
//		fprintf(temp,"%f",GPSTIME);
//		for(i=1;i<=31;i++)
//			fprintf(temp,",%f",MPG_ref[i]);
//		fprintf(temp,"\n");
//	}
}