libapri.h

实现ipvlbi数据记录

	 {
            //printf(" DZ selected\n");
            Rot[2][2]=0;           //  -S   C   0
            Rot[0][0]=-S;          //  -C  -S   0
            Rot[0][1]=C;           //   0   0   0
            Rot[1][0]=-C;
            Rot[1][1]=-S;
	 }

	 for (i=0; i<3; i++)
	 {
		 for(k=0; k<3; k++)
		 {
             Ans[i][k]=Rot[i][0]*(Org[0][k])+Rot[i][1]*(Org[1][k])+Rot[i][2]*(Org[2][k]);
		 }
	 }
	 for (i=0; i<3; i++)
	 {
		 for(k=0; k<3; k++)
		 {
             Org[i][k]=Ans[i][k];
		 }
	 }
}


void precession(double t, double Org[][3])
{
    /* precession matrix */
	double t2,t3, Guzai, Ieta, Theta;
	t2=pow(t,2);
	t3=pow(t,3);
	Guzai=(2306.2181*t+.30188*t2+.017998*t3)*CONVA;
	Ieta=(2306.2181*t+1.09468*t2+.018203*t3)*CONVA;
	Theta=(2004.3109*t-.42665*t2-.041833*t3)*CONVA;
	rotation(Ieta,Org,"Z");           // Z axis rotation
	rotation(-Theta,Org,"Y");         // Y axis rotation
	rotation(Guzai,Org,"Z");           // Z axis rotation
}

void wobble(double wobx, double woby, double Org[][3])
{
	/* Wobble Matrix */
        rotation(wobx,Org,"Y");            //! Y axis rotation
        rotation(-woby,Org,"X");           //! X axis rotation
}

void diurnal_motion(double T0, double Ut1, double D_phai,
					double Epsiron, double *Gast, double Org[][3])
{
/* **************************
  Convs TIME SEC ----> RAD
  Convd DEG      ----> RAD
  Conva arcsec   ----> RAD
  Convh HOUR     ----> RAD
 **************************
*/
      /* Diurnal Motion Matrix  */
         //  GMST : Greenwich Mean Sidereal Time
         //  GAST : Greenwich Apparent Sidereal Time
	    double T02, T03, Omega, Gmst_0h, Gmst;
		//printf("T0,Ut1,D_phai,Epsiron %g %g %g %g\n",T0,Ut1,D_phai,Epsiron);
        T02=pow(T0,2.0);
        T03=pow(T0,3.0);
        Omega=(1.00273790935080+5.9006E-11*T0-5.9E-15*T02)*CONVS;
        Gmst_0h=(24110.54841+8640184.812866*T0+.093104*T02-6.2E-6*T03)*CONVS;
        Gmst_0h=fmod(Gmst_0h, TWOPI);
        Gmst=Gmst_0h+Ut1*Omega;           // GMST (Time sec->RAD)
        *Gast=Gmst+D_phai*cos(Epsiron);    // GAST
        rotation(-(*Gast),Org,"Z");          // Z axis rotation
        //printf("Omega %g\n",Omega);
}

void m_mat33(double A[][3], double B[][3], double C[][3]) 
/* 3x3峴楍偺偐偗嶼 */
{
	//double C[3][3];
	int i,k;
	for (i=0; i<3; i++)
	{
	   for(k=0; k<3; k++)
	   {
          C[i][k]=A[i][0]*B[0][k]+A[i][1]*B[1][k]+A[i][2]*B[2][k];
	   }
	}
}

void m_mat31(double A[][3], double *B, double *C) 
/* 3x1峴楍偺偐偗嶼 */
{
	//double C[3][3];
	int i;
	for (i=0; i<3; i++)
	{
          C[i]=A[i][0]*B[0]+A[i][1]*B[1]+A[i][2]*B[2];
	}
}

void Az_el(double Longi, double Latit, double T[][3], 
		   double *Star_v, double *Az, double *El)
{
        double Tm_t[3][3], S_ans[3], R_ans[3], Rot_v[3][3];
		double d1,d2;
        int i,k;
	    for (i=0; i<3; i++)
		{
		   for(k=0; k<3; k++)
		   {
			 Rot_v[i][k]=0.0;
		   }
		}
	    Rot_v[0][0]=1.0; Rot_v[1][1]=1.0; Rot_v[2][2]=1.0;
        rotation(Longi,Rot_v,"Z");        // ! Z Axis
        rotation(-Latit,Rot_v,"Y");       //    ! Y Axis
	    for (i=0; i<3; i++)
		{
		   for(k=0; k<3; k++)
		   {
			 Tm_t[i][k]=T[k][i];
		   }
		}

	    for (i=0; i<3; i++)
		{
           R_ans[i]=Tm_t[i][0]*Star_v[0]+Tm_t[i][1]*Star_v[1]+Tm_t[i][2]*Star_v[2];
		}
	    for (i=0; i<3; i++)
		{
           S_ans[i]=Rot_v[i][0]*R_ans[0]+Rot_v[i][1]*R_ans[1]+Rot_v[i][2]*R_ans[2];
		}
        d1=S_ans[0];
        *El=asin(d1);
		d1=S_ans[1];
		d2=S_ans[2];
        *Az=atan2(d1,d2);
}

void Atm(double Hight, double El, double *Xatm)
{
/*
   !--------------------------------------------------------------
   ! **** ATM model is CHAO                                      !
   !      Zenith excess path is caluculated standard model       !
   !           Hight=0   P=1013.25 mb     15逤                   !
   !           P=P0-0.115*Hight                                  !
   !           WV PRESS = 8.5                                    !
   !           ZP=0.002277*(P+37.5)=2.393-2.62E-4*Hight (m)      !
   !             =7.98E-9-8.74E-13*Hight (sec)                   !
   !-------------------------------------------------------------
*/
	double Zp;
    Zp=7.98E-9-8.74E-13*Hight;
    *Xatm=Zp/(sin(El)+0.00143/(tan(El)+0.0445));
}

void  Uv_cal(double *T_base, double Rar, double Decr, double *Xuv)
{

     Xuv[0]=-(-T_base[0]*sin(Rar)+T_base[1]*cos(Rar))*CONVA;
     Xuv[1]=-((-T_base[0]*cos(Rar)-T_base[1]*sin(Rar))*sin(Decr)
		     +T_base[2]*cos(Decr))*CONVA;
}


void cal_apri(double Time, double *Base_v_c, double *Star_v, 
			  double *Latit, double *Longi, double *Hight,
              double Rar, double Decr, double Xwob, double Ywob,
              double Ut1_c, double Offset, double Rate,
              double Pr_time,
			  double *Xtau, double *Xaz, double *Xel,
              int *Ight, double *Xuv)
/*    Original is  KROSS    
   !
   !-------------------------------------------------------------------------
   ! INPUT PARAMETERS ;
   !  Time              : Time on a priori value (mjd in sec) [sec]
   !  Base_v_c(1;3)     : Baseline vector / C      X,Y,Z [sec]
   !  Star_V(1;3)       : Source vector
   !  Latit(1;2)        : Latitude  of each station (rad)
   !  Longi(1;2)        : Longitude of each station (rad)
   !  Hight(1;2)        : Hight     of each station (rad)
   !  Rar               : Right Ascention (rad)
   !  Decr              : Declination (rad)
   !  Xwob              : Wobbling X parameter [rad]
   !  Ywob              : Wobbling Y parameter [rad]
   !  Ut1_c             : UT1-UTC [sec]
   !  Offset,Rate       : Clock parameter ([sec],[sec/sec]) at Pr_time
   !  Pr_time           : Reference time of clock (mjd in sec) (PRT) [sec]
   ! OUTPUT PARAMETERS;
   !  Xtau              : Delay (Sec)
   !  Xaz(1;2),Xel(1;2) : Az,El (rad)
   !  Ight(1;5)         : Greenwich Hour Time
   !  Xuv(1;2)          : U,V    (Fringe/arcsec)
   !-------------------------------------------------------------------------
   !
*/
/*  [T] = [Precession] [Nutation] [Diurnal motion] [Wobbling] */
{
    //double Ini[3][3];
    double Prec[3][3],Nuta[3][3],Diur[3][3],Wobb[3][3];
    double T[3][3],T_base[3],T1[3][3],T2[3][3];
    double Star_v_cr[3];
    double Atm_del[2];
	double Utc_t, Ut1, Bbk, Bbell, Omgat, Comg, Somg;
	double Utc_0, T0, Utc, D_phai, Epsiron, Sunlon, Gast;
	double Se, Ce, Sl, Cl;
	double Xvc, Yvc, Zvc, Normal;
	double Ra_cr, Ght, Lon, Lat, Hi, Az, El, Xatm;
	double d1, dw1;
	long mjd1;
	int i,k;
	/* matrix initialize */
	for (i=0; i<3; i++)
	{
		 for(k=0; k<3; k++)
		 {
			 Nuta[i][k]=0.0;
			 if (k==i) { Nuta[i][k]=1.0; }
			 Diur[i][k]=0.0;
			 if (k==i) { Diur[i][k]=1.0; }
			 Prec[i][k]=0.0;
			 if (k==i) { Prec[i][k]=1.0; }
			 Wobb[i][k]=0.0;
			 if (k==i) { Wobb[i][k]=1.0; }		 
		 }
	}
     /* time set */
	 mjd1=ut2mjd(&d1, 2000, 1, 1, 0, 0, 0);
     dw1=((double)mjd1 + d1)*86400.0;
     Utc_t=(Time-dw1)/86400.0;
     Utc_0=(int)(Utc_t) -0.5 ;                            // 2000y0.5d
     T0=Utc_0/36525.0  ;                            // (TDB-J2000.0)/CENT.
     Ut1=fmod(Utc_t, 1.0) ;                        // UT1 FROM 00:00:00UTC
     //printf("Ut1 (before) %g\n",Ut1);
     if (Ut1<0) {
		  Ut1=1.0+Ut1;
	 }
     //printf("Ut1 (middle) Ut1_c %g %g\n",Ut1,Ut1_c);
     Ut1=Ut1*86400.0+Ut1_c;
     //printf("Ut1 (after) %20.10e\n",Ut1);
     Utc=(Utc_t-0.5)/36525.0; // (UTC-J2000.0)/CENT.
     //printf("Utc_0, Utc_t %20.10e %20.10e\n",Utc_0,Utc_t);
     //printf("mjd1, d1, dw1 %d %f %f\n",mjd1,d1,dw1);
     //printf("Time %f\n",Time);
     //printf("Utc_t, Utc %20.10e %20.10e\n",Utc_t,Utc);
     //printf("T0, Ut1 %20.10e %20.10e\n",T0,Ut1);
	 /* Precession,Nutation,Diurnal motion,Wobbling */
     precession(Utc,Prec);
     nutation(Utc,&D_phai,&Epsiron,&Sunlon,Nuta);
     diurnal_motion(T0,Ut1,D_phai,Epsiron,&Gast,Diur);
     wobble(Xwob,Ywob,Wobb);
	 /*
     printf("Utc %f\n",Utc);
	 printf("Prec %f %f %f\n",Prec[0][0],Prec[0][1],Prec[0][2]);
	 printf("     %f %f %f\n",Prec[1][0],Prec[1][1],Prec[1][2]);
	 printf("     %f %f %f\n",Prec[2][0],Prec[2][1],Prec[2][2]);
	 printf("Nuta %f %f %f\n",Nuta[0][0],Nuta[0][1],Nuta[0][2]);
	 printf("     %f %f %f\n",Nuta[1][0],Nuta[1][1],Nuta[1][2]);
	 printf("     %f %f %f\n",Nuta[2][0],Nuta[2][1],Nuta[2][2]);
     */
	 /*--------- Transfer matrix  T */
     m_mat33( Prec, Nuta ,T1) ;                           // Precession*Nutation
     m_mat33( T1, Diur, T2)  ;                            // P*N*Diurnal motion
     m_mat33( T2, Wobb ,T) ;                            // MATRIX T=PNDW

     //printf("Base_v_c %f %f %f\n",Base_v_c[0],Base_v_c[1],Base_v_c[2]);
	 m_mat31( T, Base_v_c , T_base);                  // T *Baseline Vector
     //printf("T_base %f %f %f\n",T_base[0],T_base[1],T_base[2]);

	/* ***** Abraviation and normarized source vector
              Relativistic transform from TDB to TAI
              and movement of Y station during delay
    */
     Bbk=20.49525*CONVA;
     Bbell=Bbk*(0.01675104-0.00004180*(1.0+Utc));
     Omgat=(281.22083+0.000047084*(1.0+Utc))*CONVD-PI;
     Sunlon=fmod(Sunlon, TWOPI);
	 //printf("Sunlon %f\n",Sunlon);
     Comg=cos(Omgat);
     Somg=sin(Omgat);
     Se=sin(Epsiron);
     Ce=cos(Epsiron);
     Sl=sin(Sunlon);
     Cl=cos(Sunlon);
/*
   !---------------------------------------
   ! *****  Velocity of earth (30km/s)    !
   !---------------------------------------
*/
     Xvc=Bbk*Sl-Bbell*Somg;
     Yvc=-Bbk*Cl*Ce+Bbell*Comg*Ce;
     Zvc=-Bbk*Cl*Se+Bbell*Comg*Se;
     //printf("Xvc,Yvc,Zvc %g %g %g\n",Xvc,Yvc,Zvc);
	 //printf("Star_v  %g %g %g\n",Star_v[0],Star_v[1],Star_v[2]);
     Star_v_cr[0]=Star_v[0]+Xvc;
     Star_v_cr[1]=Star_v[1]+Yvc;
     Star_v_cr[2]=Star_v[2]+Zvc;

     Normal=sqrt(pow(Star_v_cr[0],2)+
		         pow(Star_v_cr[1],2)+
				 pow(Star_v_cr[2],2));
	for (i=0; i<3; i++)
	{
       Star_v_cr[i]= Star_v_cr[i]/Normal;   // Normalized
	}
	 //printf("Star_v_cr  %g %g %g\n",Star_v_cr[0],Star_v_cr[1],Star_v_cr[2]);

     Ra_cr=atan2(Star_v_cr[1],Star_v_cr[0]) ; //Correct Right Ascention
     //printf("Ra_cr %f\n",Ra_cr);
/*
  !------------------------------------------------
  !  Greenwich hour angle   Ight  HH,MM,SS,MSS,0  !
  !------------------------------------------------
  !
*/
     Ght=Gast-Ra_cr;
     Ght=fmod(Ght, TWOPI);
     if (Ght < 0.0) { Ght=Ght+TWOPI; }
     Ght=Ght/CONVH;
     Ight[0]=(int)Ght;
     Ght=fmod(Ght,1.0)*60.0;
     Ight[1]=(int)Ght;
     Ght=fmod(Ght, 1.0)*60.0;
     Ight[2]=(int)Ght;
     Ight[3]=(int)(fmod(Ght,1.0)*1000.0);
     Ight[4]=0;
	 //printf("Gast,Ght %g %g\n",Gast,Ght);
/*
   !
   !--------------------
   !    Az,El          !
   !--------------------
   !
*/
	 for (i=0; i<2; i++)
	 {
          Lon=Longi[i];             //  Longitude
          Lat=Latit[i];             //  Latitude
          Hi=Hight[i];              //  Hight
          Az_el(Lon,Lat,T,Star_v_cr,&Az,&El); // Az,El
          Atm(Hi,El,&Xatm);         // Atmosphere Correction
          Xaz[i]=Az;
          Xel[i]=El;
          Atm_del[i]=Xatm;
		  //printf("Lon,Lat,Hi %g %g %g\n",Lon,Lat,Hi);
		  //printf("AZ,El,Xatm %g %g %g\n",Az,El,Xatm);
	 }

/*
   !
   !----------------------
   !      Delay          !
   !----------------------
   !
*/
  /* **** Geometric DELAY  (B