ef18_sp3.c

使用C语言编写的计算GPS卫星在任意时刻的三维速度的程序

/*
   fprintf(fpsp3, "%3s%4d%16.8lf%15.8lf%6ld%16.13lf\n",
            "## ", ef18_hdr.GPS_week_strt, ef18_hdr.second_of_week,
            ef18_hdr.epoch_interval, ef18_hdr.mjd_strt, ef18_hdr.fmjd_strt);
*/

   fprintf(fpsp3, "%3s%4ld%16.8lf%15.8lf%6ld%16.13lf\n",
            "## ", run.gpswk_strt, run.secofwk_strt,
            run.delta_t, run.mjd_strt, run.fmjd_strt);

   /*  PRINT HEADER LINES 3-7 TO SP3ASCII FILE  */
   /*+++++++++++++++++++++++++++++++++++++++++++*/

   fprintf(fpsp3, "%2s%4d%6d",
            "+ ", run.num_prns, run.idsv[0]);

   for(i=1; i<=16; i++) fprintf(fpsp3, "%3d", run.idsv[i]);
   fprintf(fpsp3, "\n");

   k = 17;
   for(i=1; i<=4; i++)
      { fprintf(fpsp3, "+        ");
        for(j=0; j<=16; j++) fprintf(fpsp3, "%3d", run.idsv[k++]);
        fprintf(fpsp3, "\n");
      }


   /*  PRINT HEADER LINES 8-12 TO SP3ASCII FILE */
   /*+++++++++++++++++++++++++++++++++++++++++++*/

   k = 0;
   for(i=1; i<=5; i++)
      { fprintf(fpsp3, "++       ");
        for(j=0; j<=16; j++) fprintf(fpsp3, "%3d", run.SV_acc[k++]);
        fprintf(fpsp3, "\n");
      }


   /*  PRINT HEADER LINES 13-14 TO SP3ASCII FILE */
   /*++++++++++++++++++++++++++++++++++++++++++++*/

   fprintf(fpsp3, "%.2s %.2s %.2s %.3s %.3s %.4s %.4s %.4s %.4s %.5s %.5s %.5s %.5s\n",
            "%c", ef18_hdr.ca, ef18_hdr.cb, ef18_hdr.cc, ef18_hdr.cd,
            ef18_hdr.ce, ef18_hdr.cf, ef18_hdr.cg, ef18_hdr.ch,
            ef18_hdr.ci, ef18_hdr.cj, ef18_hdr.ck, ef18_hdr.cl);

   fprintf(fpsp3, "%.2s %.2s %.2s %.3s %.3s %.4s %.4s %.4s %.4s %.5s %.5s %.5s %.5s\n",
            "%c", ef18_hdr.cm, ef18_hdr.cn, ef18_hdr.co, ef18_hdr.cp,
            ef18_hdr.cq, ef18_hdr.cr, ef18_hdr.cs, ef18_hdr.ct,
            ef18_hdr.cu, ef18_hdr.cv, ef18_hdr.cw, ef18_hdr.cx);


   /*  PRINT HEADER LINES 15-22 TO SP3ASCII FILE */
   /*++++++++++++++++++++++++++++++++++++++++++++*/

   fprintf(fpsp3, "%.2s %10.7lf %12.9lf %14.11lf %18.15lf\n",
            "%f", ef18_hdr.fa, ef18_hdr.fb, ef18_hdr.fc, ef18_hdr.fd);

   fprintf(fpsp3, "%.2s %10.7lf %12.9lf %14.11lf %18.15lf\n",
            "%f", ef18_hdr.fe, ef18_hdr.ff, ef18_hdr.fg, ef18_hdr.fh);

   fprintf(fpsp3, "%.2s %4d %4ld %4ld %4ld %6ld %6ld %6ld %6ld %9ld\n",
            "%i", ef18_hdr.ia2, ef18_hdr.ib4, ef18_hdr.ic4, ef18_hdr.id4,
            ef18_hdr.ie4, ef18_hdr.if4, ef18_hdr.ig4, ef18_hdr.ih4, ef18_hdr.ii4);

   fprintf(fpsp3, "%.2s %4d %4ld %4ld %4ld %6ld %6ld %6ld %6ld %9ld\n",
            "%i", ef18_hdr.ij2, ef18_hdr.ik4, ef18_hdr.il4, ef18_hdr.im4,
            ef18_hdr.in4, ef18_hdr.io4, ef18_hdr.ip4, ef18_hdr.iq4, ef18_hdr.ir4);

   fprintf(fpsp3, "%.2s %.18s%.18s%.18s%.3s\n%.2s %.18s%.18s%.18s%.3s\n%.2s %.18s%.18s%.18s%.3s\n%.2s %.18s%.18s%.18s%.3s\n",
            "/*", ef18_hdr.line29, ef18_hdr.line30,
                  ef18_hdr.line31, ef18_hdr.line32,
            "/*", ef18_hdr.line33, ef18_hdr.line34,
                  ef18_hdr.line35, ef18_hdr.line36,
            "/*", ef18_hdr.line37, ef18_hdr.line38,
                  ef18_hdr.line39, ef18_hdr.line40,
            "/*", ef18_hdr.line41, ef18_hdr.line42,
                  ef18_hdr.line43, ef18_hdr.line44);

} /* end of function "change_ef18_headr_lines_to_sp3" */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void gps_to_mjd( long gpswk, double secofwk, long *pmjd, double *pfmjd)

      /****************************************************|
      | void gps_to_mjd() converts a (gpsweek, seconds_of_week)
      |                   couple to a (mjd, fmjd) couple.
      |                   Speed was the primary consideration.
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
static long jan680=44244;
static double days_per_second = 0.00001157407407407407;
long dayofwk;

dayofwk = secofwk*days_per_second;
*pfmjd = (secofwk - dayofwk*86400)*days_per_second;
*pmjd = jan680 + gpswk*7 + dayofwk;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
 void gps_to_ydhms( long gpswk, double secofwk, long *pyear,
                    long *pyday, long *phour, long *pminute, double *psecond)


      /****************************************************|
      | void gps_to_ydhms() converts a (gpswk, secofwk)
      |                     couple to a (year, yday, hour,
      |                     minute, second) quintuple.  Speed
      |                     and clarity were the main goals.
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/
{
 void gps_to_mjd( long, double, long *, double *);
 void fmjd_to_hms( double, long *, long *, double *);
 static long mjd_jan1_1901 = 15385;
 long day_diff, remain, mjday;
 long elections, kyear;
 double fmjd;

 gps_to_mjd( gpswk, secofwk, &mjday, &fmjd);
 day_diff = mjday - mjd_jan1_1901;
 elections = day_diff/1461;
 *pyear = 1901 + 4*elections;
 remain = day_diff - elections*1461;     /*  For remain = 1460 we    */
 kyear = remain/365 - remain/1460;       /*  have remain/365 is 4.   */
 *pyear += kyear;                        /*  -remain/1460 fixes this.*/
 *pyday = remain - 365*kyear + 1;
 fmjd_to_hms( fmjd, phour, pminute, psecond);
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void fmjd_to_hms( double fmjd, long *phour, long *pmin, double *psec)

      /****************************************************|
      | void fmjd_to_hms() converts a fractional day to
      |                    hours, minutes, seconds. We
      |                    require that 0<=fmjd<1.
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
*phour = fmjd*24;
*pmin  = fmjd*1440  - *phour*60;
*psec  = fmjd*86400 - *phour*3600 - *pmin*60;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void month_and_day( long year, long yday, long *pmonth, long *pmday)

      /****************************************************|
      | void month_and_day() converts day_of_year to month
      |                      and day day_of_month.  The main
      |                      objective was speed...
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
 static long month_day[2][13] = {
            {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365},
            {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366}
                                };
 long leap, guess, more;

 leap = (year%4 == 0);  /* leap = true = 1 for a leap year.*/
 guess = yday*0.032;    /* 0.032 was derived empirically. */
 more = (  ( yday - month_day[leap][guess+1] )  > 0  );
        /* more computes to 0 or 1 depending on veracity. */
 *pmonth = guess + more + 1;
 *pmday = yday - month_day[leap][guess+more];

}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void mjd_to_gps( long mjd, double fmjd, long *pgpswk, double *psecofwk)

      /****************************************************|
      | void mjd_to_gps() converts a (mjd, fmjd) couple to
      |                   a (gpsweek, seconds_of_week)
      |                   couple.  Speed was the primary
      |                   consideration...
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
static long jan680=44244;
long idays, dayofwk;

idays = mjd - jan680;
*pgpswk = idays/7;
dayofwk = idays - (*pgpswk)*7;
*psecofwk = (dayofwk + fmjd)*86400;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void mjd_to_ydhms(long mjday, double fmjd, long *pyear,
               long *pyday, long *phour, long *pminute, double *psecond)

      /****************************************************|
      | void mjd_to_ydhms() converts a (mjd, fmjd) couple
      |                        to a (year, yday, hour,
      | mjdate to yearday      minute, second) quintupal.
      | -----------------      I believe it is efficient..
      |
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
 void fmjd_to_hms(double, long *, long *, double *);
 static long mjd_jan1_1901 = 15385;
 long day_diff, remain;
 long elections, kyear;

 day_diff = mjday - mjd_jan1_1901;
 elections = day_diff/1461;
 *pyear = 1901 + 4*elections;
 remain = day_diff - elections*1461;     /*  For remain = 1460 we    */
 kyear = remain/365 - remain/1460;       /*  have remain/365 is 4.   */
 *pyear += kyear;                        /*  -remain/1460 fixes this.*/
 *pyday = remain - 365*kyear + 1;
 fmjd_to_hms(fmjd, phour, pminute, psecond);
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void mjd_to_ymdhms(long mjday, double fmjd, long *pyear, long *pmonth,
                   long *pmday, long *phour, long *pminute, double *psecond)

      /****************************************************|
      | void mjd_to_ymdhms() converts a (mjd, fmjd) couple
      |                      to a (year, month, mday, hour,
      | mjdate to greg.      minute, second) sextupal.  I
      | ---------------      believe it is efficient..
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
 long yday;
 void mjd_to_ydhms(long, double, long *, long *, long *, long *, double *);
 void month_and_day(long, long, long *, long *);

 mjd_to_ydhms(mjday, fmjd, pyear, &yday, phour, pminute, psecond);
 month_and_day(*pyear, yday, pmonth, pmday);
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void ydhms_to_gps( long year, long yday, long hour, long minute,
                   double second, long *pgpswk, double *psecofwk)

      /****************************************************|
      |  void ydhms_to_gps() converts a (year, yday, hour,
      |                      minute, second) quintuple to a
      |                      (gpswk, secofwk) couple.
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
void ydhms_to_mjd( long, long, long, long, double, long *, double *);
void mjd_to_gps( long, double, long *, double *);
long mjday;
double fmjd;

ydhms_to_mjd( year, yday, hour, minute, second, &mjday, &fmjd);
mjd_to_gps( mjday, fmjd, pgpswk, psecofwk);
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
 void ydhms_to_mjd(long year, long yday, long hour, long minute,
    double second, long *pmjd, double *pfmjd)

      /****************************************************|
      | void ydhms_to_mjd() converts a (year, day, hour,
      |                     minute, second) quintuple to a
      |                     (mjday, fmjd) couple.  Speed and
      |                     clarity were the main goals...
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
 static long mjd_jan1_1901 = 15385;
 static double days_per_second = 0.00001157407407407407;
 long years_into_election, elections;

 years_into_election = (year-1)%4;
 elections = (year - 1901)/4;
 *pfmjd = ( hour*3600 + minute*60 + second)*days_per_second;
 *pmjd = mjd_jan1_1901 + elections*1461 + years_into_election*365 + yday -1;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void ymdhms_to_gps(long year, long month, long mday, long hour,
             long minute, double second, long *pgpswk, double *psecofwk)

      /****************************************************|
      | void ymdhms_to_gps() converts a (year, month, mday,
      |                      hour, minute, second) sextuple
      |                      to a (gpswk, secofwk) couple.
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
 static double days_per_second = 0.00001157407407407407;
 long grday_to_mjday(long, long, long);
 void mjd_to_gps(long, double, long *, double *);
 long mjday;
 double fmjd;

 mjday = grday_to_mjday(year, month, mday);
 fmjd = ( hour*3600 + minute*60 + second)*days_per_second;
 mjd_to_gps(mjday, fmjd, pgpswk, psecofwk);
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
long grday_to_mjday( long year, long month, long mdays)

      /****************************************************|
      | long grday_to_mjday converts year, month and day_
      |                     of_month to modified julian
      |                     days.  The main objective
      |                     was speed...
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
 long day_of_year(long, long, long);
 static long mjd_jan1_1901 = 15385;
 long days_from_jan1, years_into_election, elections;

 days_from_jan1 = day_of_year( year, month, mdays) - 1;
 years_into_election = (year-1)%4;
 elections = (year - 1901)/4;
 return( mjd_jan1_1901 + elections*1461 +
                         years_into_election*365 + days_from_jan1 );
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
long day_of_year(long year, long month, long mday)

      /****************************************************|
      | long day_of_year() converts month and day_of_month
      |                   to day_of_year.  The main objective
      |                   was speed...
      |
      |            [ B.W.Remondi - October 1987 ]
      |****************************************************/

{
 static long month_day[2][12] = {
            {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334},
            {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}
                                };
 long leap;
 leap = (year%4 == 0);  /* leap = true = 1 for a leap year.*/
 return( month_day[leap][month-1] + mday );
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void round_60_seconds( long *yr, long *mon, long *mday,
                       long *hr, long *min, double *second)
{
 void ymdhms_to_gps(long year, long month, long mday, long hour,
      long minute, double second, long *pgpswk, double *psecofwk);
 void gps_to_mjd( long gpswk, double secofwk, long *pmjd, double *pfmjd);
 void ajmjd(long mjdayo, double fmjdayo, double seconds,
           long *pmjdn, double *pfmjdn);
 void mjd_to_ymdhms(long mjd, double fmjd, long *pyear, long *pmonth,
               long *pmday, long *phour, long *pminute, double *psecond);
 long gpswk, mjd, mjdn;
 double secofwk, fmjd, fmjdn;

 ymdhms_to_gps(*yr, *mon, *mday, *hr, *min, (double)59.0, &gpswk, &secofwk);
 gps_to_mjd( gpswk, secofwk, &mjd, &fmjd);
 ajmjd( mjd, fmjd, (double)2.0, &mjdn, &fmjdn);
 mjd_to_ymdhms(mjdn, fmjdn, yr, mon, mday, hr, min, second);
 *second = 0.0;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void ajmjd(long mjdayo, double fmjdayo, double seconds,
           long *pmjdn, double *pfmjdn)

      /****************************************************|
      |
      |  (mjdold,fmjdold) + seconds/86400 = (mjdnew,fmjdnew)
      |
      |  This function allows one to increment
      |  a (mjday, fmjday) couple by amount seconds.
      |  seconds can be positive or negative.
      |
      |  This subroutine should always be used to