////////////////////////////////////////////////////////////////
//
//	人工的にコード位相を生成
//
////////////////////////////////////////////////////////////////

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

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

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

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

extern FILE *dd_code,*mp_code;

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,MP_FLAG,Correct_wide_integer[SATMAX],Correct_uwide_integer[SATMAX];
extern int Integer_r2[SATMAX],Integer_m2[SATMAX],DD_integer2[SATMAX];
extern int Integer_r5[SATMAX],Integer_m5[SATMAX],DD_integer5[SATMAX];


extern double Cp1ref[SATMAX],Cp1mov[SATMAX],Cp2ref[SATMAX],Cp2mov[SATMAX],Cp5ref[SATMAX],Cp5mov[SATMAX];
extern double DD_carrier1[SATMAX],DD_carrier2[SATMAX],DD_carrier5[SATMAX];
extern double Pr1ref[SATMAX],Pr1mov[SATMAX],Pr2ref[SATMAX],Pr2mov[SATMAX];
extern double Pr5ref[SATMAX],Pr5mov[SATMAX];
extern double DD_code1[SATMAX],DD_code2[SATMAX],DD_code5[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],Cn2[SATMAX],GPSTIME,MP_rov[SATMAX],MPG_rov[SATMAX],MPG_ref[SATMAX];
extern double MPG_rov5[SATMAX],MPG_ref5[SATMAX];
extern double MP_rov[SATMAX],MP_rov2[SATMAX],MP_rov3[SATMAX];
extern double MP_rov5[SATMAX],MP_rov5_2[SATMAX],MP_rov5_3[SATMAX];
extern double OBST1_MASK,Ultra_wide_carrier[SATMAX];

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

#define pow(x, y) pow((double)x, (double)y)

void make_code()
{
	int i,k=0;
	int sat,prn;
	int arr[1];
	double wn,cn[SATMAX];
	double p_sigma[SATMAX];
	double mp[SATMAX]={0};
	double noise_bw_code=0.05;
	double correlator_chip=0.1;
	double correlator_factor=0.5;
	double mp1[SATMAX]={0},mp2[SATMAX]={0};
	double wide_carrier[SATMAX];
	double noise_factor = 1.5;

	static int incr_code1[36]={1,1,-10,10,-20,20,-30,30,-40,40,-50,50,-60,60,-70,70,-80,80,-90,90,-100,100,-110,110,-120,120,-130,130,-140,140,-150,150,-160,160,-170,170};
	static int incr_code2[36]={1,10,-100,100,-200,200,-300,300,-400,400,-500,500,-600,600,-700,700,-800,800,-900,900,-1000,1000,-1100,1100,-1200,1200,-1300,1300,-1400,1400,-1500,1500,-1600,1600,-1700,1700};
	static int incr_code3[36]={1,-1,10,-10,20,-20,30,-30,40,-40,50,-50,60,-60,70,-70,80,-80,90,-90,100,-100,110,-110,120,-120,130,-130,140,-140,150,-150,160,-160,170,-170};
	static int incr_code4[36]={1000,-900,800,700,-600,-500,400,-300,200,-100,1000,-900,800,-700,600,-500,-400,300,-200,100,1000,-900,800,-700,600,-500,400,-300,-200,100,234,567,345,789,543,321};
	static int incr_code5[36]={34,-23,12,-45,56,-67,23,-54,83,-15,49,-42,58,-59,32,-19,43,-45,10,-98,103,-984,234,-765,298,-613,579,-982,346,-129,873,-459,285,-578,235,-653};
	static int incr_code6[36]={1030,-940,860,710,-660,-580,430,-350,280,-110,1450,-980,820,-760,630,-570,-410,330,-280,130,1450,-960,810,-700,660,-530,480,-320,-240,180,224,577,375,719,573,381};

////////////////////////   L1   ////////////////////////
	//基準局のコード位相の生成
	for(i=0;i<Enable[0][0];i++){
		sat = Enable_sat[0][0][i];

		//multipathを無視する場合
//		MP_rov[sat]=0;
//		MPG_ref[sat]=0;
//		MPG_rov[sat]=0;

//		MP_rov5[sat]=0;
//		MPG_ref5[sat]=0;
//		MPG_rov5[sat]=0;

		//DLLの1σ値を利用
		//白色雑音の計算に問題あり（衛星ごとに単独で発生させる）

		incr_code1[sat]++;
		if(incr_code1[sat] == 0)
			incr_code1[sat]++;

		arr[0] = incr_code1[sat];
		cn[sat]=Cn[sat]/10;
		cn[sat]=pow(10,cn[sat]);
		p_sigma[sat] = 293.05*sqrt(correlator_chip*noise_bw_code*correlator_factor/cn[sat]);
		p_sigma[sat] = p_sigma[sat]*noise_factor;
		wn = gauss_gen(0,p_sigma[sat],arr);
		mp1[sat]=MPG_ref[sat];

		//L1のマルチパス誤差を求める

		Pr1ref[sat] = (Sat_xpos[sat]-POSx[0][0])*(Sat_xpos[sat]-POSx[0][0])
			+ (Sat_ypos[sat]-POSy[0][0])*(Sat_ypos[sat]-POSy[0][0])
			+ (Sat_zpos[sat]-POSz[0][0])*(Sat_zpos[sat]-POSz[0][0]);
		Pr1ref[sat] = sqrt(Pr1ref[sat]);
		Pr1ref[sat] = Pr1ref[sat]+wn+MPG_ref[sat];
	}

	//移動局のコード位相の生成
	for(i=0;i<Enable[DIV1][DIV1];i++){
		sat = Enable_sat[DIV1][DIV1][i];

		//DLLの1σ値を利用
		incr_code2[sat]++;
		if(incr_code2[sat] == 0)
			incr_code2[sat]++;
		arr[0] = incr_code2[sat];
		cn[sat]=Cn[sat]/10;
		cn[sat]=pow(10,cn[sat]);
		p_sigma[sat] = 293.05*sqrt(correlator_chip*noise_bw_code*correlator_factor/cn[sat]);
		p_sigma[sat] = p_sigma[sat]*noise_factor;
		wn = gauss_gen(0,p_sigma[sat],arr);

		Pr1mov[sat] = (Sat_xpos[sat]-POSx[DIV1][DIV1])*(Sat_xpos[sat]-POSx[DIV1][DIV1])
			+ (Sat_ypos[sat]-POSy[DIV1][DIV1])*(Sat_ypos[sat]-POSy[DIV1][DIV1])
			+ (Sat_zpos[sat]-POSz[DIV1][DIV1])*(Sat_zpos[sat]-POSz[DIV1][DIV1]);
		Pr1mov[sat] = sqrt(Pr1mov[sat]);
		Pr1mov[sat] = Pr1mov[sat]+wn+MP_rov[sat]+MPG_rov[sat]+MP_rov2[sat]+MP_rov3[sat];

//		mp[sat]=wn+MP_rov[sat]+MPG_rov[sat];
//		mp[sat]=wn;
		mp[sat]=MPG_rov[sat];
	}

	//二重位相差の計算
	for(i=0;i<Com_satn;i++){
		prn = Com_sat[i];
		if(Com_sat[i] != Ref_sat){
			DD_code1[prn] = (Pr1mov[prn] - Pr1ref[prn])
				-(Pr1mov[Ref_sat] - Pr1ref[Ref_sat]);
		}
	}

////////////////////////   L2   ////////////////////////
	//基準局のコード位相の生成(L2)
	for(i=0;i<Enable[0][0];i++){
		sat = Enable_sat[0][0][i];

		//DLLの1σ値を利用
		incr_code3[sat]++;
		if(incr_code3[sat] == 0)
			incr_code3[sat]++;
		arr[0] = incr_code3[sat];
		cn[sat]=Cn2[sat]/10;
		cn[sat]=pow(10,cn[sat]);
		p_sigma[sat] = 293.05*sqrt(correlator_chip*noise_bw_code*correlator_factor/cn[sat]);
		p_sigma[sat] = p_sigma[sat]*noise_factor;
		wn = gauss_gen(0,p_sigma[sat],arr);

		Pr2ref[sat] = (Sat_xpos[sat]-POSx[0][0])*(Sat_xpos[sat]-POSx[0][0])
			+ (Sat_ypos[sat]-POSy[0][0])*(Sat_ypos[sat]-POSy[0][0])
			+ (Sat_zpos[sat]-POSz[0][0])*(Sat_zpos[sat]-POSz[0][0]);
		Pr2ref[sat] = sqrt(Pr2ref[sat]);
		Pr1ref[sat] = Pr2ref[sat]+wn+MPG_ref[sat];
	}

	//移動局のコード位相の生成
	for(i=0;i<Enable[DIV1][DIV1];i++){
		sat = Enable_sat[DIV1][DIV1][i];

		//DLLの1σ値を利用
		incr_code4[sat]++;
		if(incr_code4[sat] == 0)
			incr_code4[sat]++;

		arr[0] = incr_code4[sat];
		cn[sat]=Cn2[sat]/10;
		cn[sat]=pow(10,cn[sat]);
		p_sigma[sat] = 293.05*sqrt(correlator_chip*noise_bw_code*correlator_factor/cn[sat]);
		p_sigma[sat] = p_sigma[sat]*noise_factor;
		wn = gauss_gen(0,p_sigma[sat],arr);

//		mp[sat]=wn;

		Pr2mov[sat] = (Sat_xpos[sat]-POSx[DIV1][DIV1])*(Sat_xpos[sat]-POSx[DIV1][DIV1])
			+ (Sat_ypos[sat]-POSy[DIV1][DIV1])*(Sat_ypos[sat]-POSy[DIV1][DIV1])
			+ (Sat_zpos[sat]-POSz[DIV1][DIV1])*(Sat_zpos[sat]-POSz[DIV1][DIV1]);
		Pr2mov[sat] = sqrt(Pr2mov[sat]);
		Pr2mov[sat] = Pr2mov[sat]+wn+MPG_rov[sat];
	}
	//二重位相差の計算
	for(i=0;i<Com_satn;i++){
		prn = Com_sat[i];
		if(Com_sat[i] != Ref_sat){
			DD_code2[prn] = (Pr2mov[prn] - Pr2ref[prn])
				-(Pr2mov[Ref_sat] - Pr2ref[Ref_sat]);
		}
	}

////////////////////////   L5   ////////////////////////
	//基準局のコード位相の生成(L5)
	for(i=0;i<Enable[0][0];i++){
		sat = Enable_sat[0][0][i];

		//DLLの1σ値を利用
		incr_code5[sat]++;
		if(incr_code5[sat] == 0)
			incr_code5[sat]++;
		arr[0] = incr_code5[sat];
		cn[sat]=Cn2[sat]/10;
		cn[sat]=pow(10,cn[sat]);
		p_sigma[sat] = 29.305*sqrt(correlator_chip*10*noise_bw_code*correlator_factor/cn[sat]);
		p_sigma[sat] = p_sigma[sat]*noise_factor;
		wn = gauss_gen(0,p_sigma[sat],arr);

		mp2[sat]=wn;

		Pr5ref[sat] = (Sat_xpos[sat]-POSx[0][0])*(Sat_xpos[sat]-POSx[0][0])
			+ (Sat_ypos[sat]-POSy[0][0])*(Sat_ypos[sat]-POSy[0][0])
			+ (Sat_zpos[sat]-POSz[0][0])*(Sat_zpos[sat]-POSz[0][0]);
		Pr5ref[sat] = sqrt(Pr5ref[sat]);
		Pr5ref[sat] = Pr5ref[sat]+wn+MPG_ref5[sat];
	}

	//移動局のコード位相の生成
	for(i=0;i<Enable[DIV1][DIV1];i++){
		sat = Enable_sat[DIV1][DIV1][i];

		//DLLの1σ値を利用
		incr_code6[sat]++;
		if(incr_code6[sat] == 0)
			incr_code6[sat]++;

		arr[0] = incr_code6[sat];
		cn[sat]=Cn2[sat]/10;
		cn[sat]=pow(10,cn[sat]);
		p_sigma[sat] = 29.305*sqrt(correlator_chip*10*noise_bw_code*correlator_factor/cn[sat]);
		p_sigma[sat] = p_sigma[sat]*noise_factor;
		wn = gauss_gen(0,p_sigma[sat],arr);

		Pr5mov[sat] = (Sat_xpos[sat]-POSx[DIV1][DIV1])*(Sat_xpos[sat]-POSx[DIV1][DIV1])
			+ (Sat_ypos[sat]-POSy[DIV1][DIV1])*(Sat_ypos[sat]-POSy[DIV1][DIV1])
			+ (Sat_zpos[sat]-POSz[DIV1][DIV1])*(Sat_zpos[sat]-POSz[DIV1][DIV1]);
		Pr5mov[sat] = sqrt(Pr5mov[sat]);
		Pr5mov[sat] = Pr5mov[sat]+wn+MPG_rov5[sat]+MP_rov5[sat]+MP_rov5_2[sat]+MP_rov5_3[sat];
	}
	//二重位相差の計算
	for(i=0;i<Com_satn;i++){
		prn = Com_sat[i];
		if(Com_sat[i] != Ref_sat){
			DD_code5[prn] = (Pr5mov[prn] - Pr5ref[prn])
				-(Pr5mov[Ref_sat] - Pr5ref[Ref_sat]);
		}
	}



	//widelaneの二重位相差もここで求めている
	//二重位相差の計算
	for(i=0;i<Com_satn;i++){
		prn = Com_sat[i];
		if(Com_sat[i] != Ref_sat){
			wide_carrier[prn] = DD_carrier1[prn]*F1/CS - DD_carrier2[prn]*F2/CS-DD_code1[prn]*(F1-F2)/CS;
			Ultra_wide_carrier[prn] = DD_carrier2[prn]*F2/CS - DD_carrier5[prn]*F5/CS-DD_code2[prn]*(F2-F5)/CS;
		}
	}
	//整数値バイアスの二重差を計算
	for(i=0;i<Com_satn;i++){
		prn = Com_sat[i];
		if(Com_sat[i] != Ref_sat){
			//正しい整数値バイアスは？
			Correct_wide_integer[prn] = DD_integer1[prn]-DD_integer2[prn];
			Correct_uwide_integer[prn] = DD_integer2[prn]-DD_integer5[prn];
		}
	}



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