libapri.h

实现ipvlbi数据记录

/**************************************************************/
/*                                                            */
/*   俬俹亅倁俴俛俬梡僀儞僋儖乕僪僼傽僀儖 libapri.h           */
/*                                                            */
/*      Version:    2002-06-07                                */
/*                                                            */
/*      Copyright (c) 2002 T.Kondo/CRL All Right Reserved     */
/*                                                            */
/**************************************************************/
/*                                   */
/* 俬俹亅倁俴俛俬梡僿僢僟乕僼傽僀儖  */
/*    傾僾儕僆儕抣寁嶼梡              */
/*         by T.Kondo 2002.6.7       */
/*                    2002.5.26      */
/*                                   */
#ifndef _LIBAPRI_
#define _LIBAPRI_

#include <math.h>

#ifndef CONVA
#define CONVA 4.848136811E-6    // PI*3600/180 arcsec to radian
#endif
#ifndef CONVS
#define CONVS 7.272205216E-5    // PI/(3600*12) time sec to radian
#endif
#ifndef CONVD
#define CONVD 1.745329252E-2    // PI/180  degree to radian
#endif
#ifndef CONVH
#define CONVH 2.617993878E-1    // PI/12 hour to radian
#endif
#ifndef PIRA
#define PIRA  1.745329252E-2    // PI/180  degree to radian
#endif
#ifndef PIDE
#define PIDE  57.29577951       // 180/PI radian to degree
#endif
#ifndef TWOPI
#define TWOPI 6.283185307       // 2PI
#endif
#ifndef PI
#define PI    3.141592653
#endif
#ifndef VLIGHT
#define VLIGHT 2.99792458E+8    // speed of light (m/sec)
#endif

#ifndef MAXSTAR
#define MAXSTAR 1500            // Maximum star numbers in schedule file
#endif

#ifndef MAXOBS
#define MAXOBS 2000            // Maximum observation (scan) numbers is schedule file
#endif

#ifndef FILEMAXSIZE
#define FILEMAXSIZE  2000000000   //  File maximum limit size in bytes
#endif

void mjd2ut(int *year, int *month, int *day, int *hour,
           int *min, int *sec, long j, double d);
long ut2mjd(double *d, int year, int month, int day,
           int hour, int min, int sec);

void diurnal_motion(double T0, double Ut1, double D_phai,
			double Epsiron, double *Gast, double Org[][3]);
void nutation(double t, double *D_phai, double *epsiron0,
              double *sunlon, double Org[][3]);
void precession(double t, double Org[][3]);
void rotation(double theta, double org[][3], char* axis);
void wobble(double wobx, double woby, double Org[][3]);
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);
void m_mat33(double A[][3], double B[][3], double C[][3]) ;
void m_mat31(double A[][3], double *B, double *C); 
void Atm(double Hight, double El, double *Xatm);
void Uv_cal(double *T_base, double Rar, double Decr, 
			 double *Xuv);
void Xyz_to_llh(double *Xyz, double *Longi, 
				double *Latit, double *Hight);
void Az_el(double Longi, double Latit, double T[][3], 
		   double *Star_v, double *Az, double *El);

typedef struct apriinfo
{
     char star[11];     // star name
     double ra;         // star position RA  in degree
     double dec;        //       "       DEC in degree
     double epo;        //       "       Epoch

     char stations[2][11];  // X,Y station name
     double xyz[2][3];     // X,Y station position (x,y,z) in m

	 double coffset;    // Clock offset at PRT  sec
	 double crate;      // Clock rate at PRT    s/s
	 double ut1_c;
	 double wobbx;
	 double wobby;

	 int    year;       // PRT year
	 int    month;      // PRT month
	 int    mday;       // PRT day of month
	 int    ddd;        // PRT day of year
	 int    hh;         // PRT hour
	 int    mm;         // PRT minutes
	 int    ss;         // PRT sec

	 double prt;        // PRT in mjd in seconds
	 int dura;       // obs duration in sec   
	 double tauap[4];   // apriori at PRT
	 int ight[5];       // greenich siderial time at PRT

     double lon[2];     // X,Y station longitude in deg
     double lat[2];     // X,Y station latitude in deg
     double h[2];     // X,Y station hight in m
	 double az[2];     // star direction (Azimuth) in deg at PRT
	 double el[2];     // star direction (Elevation) in deg at PRT

 } apri_t;
void apriori(apri_t *apri);


/* UTC偐傜弨儐儕僂僗擔傊偺曄姺 */
long ut2mjd(double *d, int year, int month, int day,
           int hour, int min, int sec)
/*
	 ut2mjd   娭悢抣丂弨儐儕僂僗擔
     d      乮弌椡乯侾擔偺帪娫傪 1 偲偟偨偲偒偺惓屵偐傜偺帪娫
     year   乮擖椡乯擭
     month  乮擖椡乯寧 (1 乣 12)
     day    乮擖椡乯擔
     hour   乮擖椡乯帪 (0 乣 23)
     min    乮擖椡乯暘 (0 乣 59)
     sec    乮擖椡乯昩 (0 乣 59)
*/
{
    int  bc, gregory, mnth, y4; 
    long j;

	mnth=month;
	y4=year;
    bc = year <= 0;
    
    /* First day at Gregory's calendar */
    gregory = year > 1582 || (year == 1582 && month > 10) ||
              (year == 1582 && month == 10 && day >= 15);
            
    if (mnth <= 2) {
        y4--;
        mnth += 12;
    }
    
    if (hour < 12) {
        j = 0;
        *d = 0.5;
    } else {
        j = 1;
        *d = -0.5;
    }
    *d += (hour*3600L + min*60 + sec)/86400.0;
    
    j += bc ? (y4-3)/4 : y4/4;
    if (gregory) j+= 2 - y4/100 + y4/400;
    j += 1720994L + y4*365L + (mnth+1)*30 + (mnth+1)*3/5 + day;
    
	/* convert to MJD */
	if (hour < 12) {
		*d -= 0.5;
	} else {
		*d += 0.5;
		j -= 1L;
	}
    j -= 2400000L;
    return j;
}

/* 弨儐儕僂僗擔傪倀俿俠偵曄姺丂*/
void mjd2ut(int *year, int *month, int *day,
		   int *hour, int *min, int *sec, long j, double d)
/*
     year   乮弌椡乯擭
     month  乮弌椡乯寧
     day    乮弌椡乯擔
     hour   乮弌椡乯帪 (0 乣 23)
     min    乮弌椡乯暘 (0 乣 59)
     sec    乮弌椡乯昩 (0 乣 59)
     j      乮擖椡乯弨儐儕僂僗擔偺惍悢晹暘
     d      乮擖椡乯弨儐儕僂僗擔偺彫悢晹暘
*/
{
    long c, k, e, s;

    /* convert to JD */
    if( d >= 0.5) {
		d -= 0.5;
		j++;
	} else {
		d += 0.5;
	}
	j += 2400000L;

	/* JD to UTC */
    if (d >= 0.5) {
        j++;
        d -= 0.5;
    } else d += 0.5;

    if (j >= 2299161) j = j + 1 + (j-1867216)/36524 - (j-1867216)/146096;
    j += 1524;

    c = (long)((j - 122.1)/365.25);
    k = 365L*c + c/4;
    e = (long)((j - k)/30.6001);

    *year = c - 4716;
    *month = e - 1;
    if (*month > 12) {
       (*month) -= 12;
       (*year)++;
    }
    *day = j - k - (long)(30.6*e);

    s = (long)(d * 86400 + 0.5);
    *hour = s / 3600;
    *min = (s % 3600) / 60;
    *sec = s % 60;
}


void nutation(double t, double *D_phai, double *epsiron0,
              double *sunlon, double Org[][3])
{
    static int n_tab[13][7]=                      /* Nutation Table */
               {{-171996,92025, 0, 0, 0, 0, 1},
                { -13187, 5736, 0, 0, 2,-2, 2},
                {  -2274,  977, 0, 0, 2, 0, 2},
                {   2062, -895, 0, 0, 0, 0, 2},
                {   1426,   54, 0, 1, 0, 0, 0},
                {    712,   -7, 1, 0, 0, 0, 0},
                {   -517,  224, 0, 1, 2,-2, 2},
                {   -386,  200, 0, 0, 2, 0, 1},
                {   -301,  129, 1, 0, 2, 0, 2},
                {    217,  -95, 0,-1, 2,-2, 2},
                {   -158,   -1, 1, 0, 0,-2, 0},
                {    129,  -70, 0, 0, 2,-2, 1},
                {    123,  -53,-1, 0, 2, 0, 2}};
     int i;
	 double t2,t3,L,Ld,F,D,Omega,D_epsiron;
	 double S_omega;

	 /* 揤暥婎杮妏 arcsec */
	 t2=pow(t,2);
	 t3=pow(t,3);
     L=485866.733+fmod(1325*t,1.0)*1296000+715922.633*t+31.31*t2+.064*t3;
     Ld=1287099.804+fmod(99*t,1.0)*1296000+1292581.224*t-.577*t2-.012*t3;
     F=335778.877+fmod(1342*t,1.0)*1296000+295263.137*t-13.257*t2+.011*t3;
     D=1072261.307+fmod(1236*t,1.0)*1296000+1105601.328*t-6.891*t2+.019*t3;
     Omega=450160.280-fmod(5*t,1.0)*1296000-482890.539*t+7.455*t2+.008*t3;

	 
	 /*   J.Wahr    */

     *D_phai=0.0;
     D_epsiron=0.0;
     for (i=0; i<13; i++)
	 {
        S_omega=(n_tab[i][2]*L+n_tab[i][3]*Ld+n_tab[i][4]*F
		    +n_tab[i][5]*D+n_tab[i][6]*Omega)*CONVA;
        *D_phai=*D_phai+n_tab[i][0]/1.0e4*sin(S_omega);
        D_epsiron=D_epsiron+n_tab[i][1]/1.e4*cos(S_omega);
	 }

	 *D_phai=*D_phai*CONVA;
     D_epsiron=D_epsiron*CONVA;
     *epsiron0=(8.4381448E+4-46.815*t-5.9E-4*t2+1.813E-3*t3)*CONVA;
     *sunlon=(F-D+Omega)*CONVA;
     //printf("epsiron0,sunlon %f %f\n",*epsiron0,*sunlon);
	 //printf("org before %f %f %f\n",Org[0][0],Org[0][1],Org[0][2]);
	 //printf("           %f %f %f\n",Org[1][0],Org[1][1],Org[1][2]);
	 //printf("           %f %f %f\n",Org[2][0],Org[2][1],Org[2][2]);

     rotation(((*epsiron0)+D_epsiron),Org,"X");   // X axis rotation
     rotation(*D_phai,Org,"Z");                // Z axis rotation
     rotation(-*epsiron0,Org,"X");    // X axis rotation
  
	 //printf("org after  %f %f %f\n",Org[0][0],Org[0][1],Org[0][2]);
	 //printf("           %f %f %f\n",Org[1][0],Org[1][1],Org[1][2]);
	 //printf("           %f %f %f\n",Org[2][0],Org[2][1],Org[2][2]);

}

void rotation(double theta, double Org[][3], char* axis)
{
     double Rot[3][3],Ans[3][3];
	 double S,C;
	 int i,k;
	 for (i=0; i<3; i++)
	 {
		 for(k=0; k<3; k++)
		 {
			 Rot[i][k]=0.0;
		 }
	 }
	 Rot[0][0]=1.0; Rot[1][1]=1.0; Rot[2][2]=1.0;
     //printf("Theta (in) = %f\n",theta);
	 S=sin(theta);
	 C=cos(theta);
	 //printf("C,S %f %f\n",C,S);
     if(strncmp(axis,"X",1)==0)    //                   X Axis Rotation
	 {
            //printf(" X selected\n");
		    Rot[1][1]=C;           //
            Rot[1][2]=S;           //   1   0   0
            Rot[2][1]=-S;          //   0   C   S
            Rot[2][2]=C;           //   0  -S   C
	 }
	 else if(strncmp(axis,"Y",1)==0)  //                Y Axis Rotation
	 {
            //printf(" Y selected\n");
            Rot[0][0]=C;           //
            Rot[2][0]=S;           //   C   0  -S
            Rot[0][2]=-S;          //   0   1   0
            Rot[2][2]=C;           //   S   0   C
	 }
     else if(strncmp(axis,"Z",1)==0)  //                Z Axis Rotation
	 {
            //printf(" Z selected\n");
            Rot[0][0]=C;           //
            Rot[0][1]=S;           //   C   S   0
            Rot[1][0]=-S;          //  -S   C   0
            Rot[1][1]=C;           //   0   0   1
	 }
     else if(strncmp(axis,"DX",2)==0)  //                X Axis Rotation
	 {
            //printf(" DX selected\n");
            Rot[0][0]=0;           //   0   0   0
            Rot[1][1]=-S;          //   0  -S   C
            Rot[1][2]=C;           //   0  -C  -S
            Rot[2][1]=-C; 
            Rot[2][2]=-S;
	 }
     else if(strncmp(axis,"DY",2)==0)  //                Y Axis Rotation
	 {
            //printf(" DY selected\n");
            Rot[1][1]=0;           //  -S   0  -C
            Rot[0][0]=-S;          //   0   0   0
            Rot[0][2]=-C;          //   C   0  -S
            Rot[2][0]=C;
            Rot[2][2]=-S;
	 }
     else if(strncmp(axis,"DZ",2)==0)  //                Y Axis Rotation