ef18_sp3.c

GPS卫星速度计算程序

      |  assure that forward or backward day crossovers
      |  have been properly computed.
      |
      |            [ B.W.Remondi - November 87 ]
      *****************************************************/

{
  static  long idays;
  static  double days;
  static double days_per_second = 0.00001157407407407407;

  days = fmjdayo + seconds*days_per_second;
  idays = days;
  *pmjdn = mjdayo + idays;
  *pfmjdn = days - idays;
  if( *pfmjdn < 0 )
     {
      *pmjdn  -= 1;
      *pfmjdn += 1;
     }
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void change_ef18_epoch_lines_to_sp3(void)
{
    void ajmjd(long, double, double, long *, double *);
    void mjd_to_ymdhms(long, double, long *, long *, long *, long *, long *, double *);
    int ef13_bwr17(int, long, double, double recf[], double vecf[]);
    long mjd, yr, mon, mday, hr, min;
    double fmjd, sec, seconds_from_start;
    double recf[4], vecf[4];
    int ii, jj, error_message;
    register long i;

    for(i=0; i<run.num_epochs; i++)
    { 
		printf("\r[%lu]", i+1);
        seconds_from_start = i * run.delta_t;
        ajmjd(run.mjd_strt, run.fmjd_strt, seconds_from_start, &mjd, &fmjd);
        if(fabs(fmjd - 1.00) < 1.0e-12 ) { mjd++; fmjd = 0.0; }
        mjd_to_ymdhms(mjd, fmjd, &yr, &mon, &mday, &hr, &min, &sec);

        if(fabs(60.0-sec) < 1.0e-8)
			round_60_seconds( &yr, &mon, &mday, &hr, &min, &sec);

        fprintf(fpsp3, "*  %4ld %2ld %2ld %2ld %2ld %11.8lf\n", yr, mon, mday, hr, min, sec);

        for(ii=0; ii<run.num_prns; ii++)
        { 
             error_message = ef13_bwr17( run.idsv[ii], mjd, fmjd, recf, vecf);
             if(error_message != 0)
             { 
				 for(jj=0; jj<3; jj++) 
				 {
					 recf[jj] = 0.0;
					 vecf[jj] = 0.0; 
				 }
				 recf[3] = vecf[3] = 999999.999999;
             }
             fprintf(fpsp3, "P%3d%14.6lf%14.6lf%14.6lf%14.6lf\n",
                      run.idsv[ii], recf[0], recf[1], recf[2], recf[3]);

			 if(run.pva=='V' && ef18_hdr.pva=='P') 
			 {
					fprintf(fpsp3, "V%3d%14.6lf%14.6lf%14.6lf%14.6lf\n",
					run.idsv[ii], 10000*vecf[0],10000*vecf[1],
					10000*vecf[2],10000*vecf[3] );
			 }

        }
    }

}  /*  End of function  */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
int ef13_bwr17(int isv, long mjdt, double fmjdt, double recf[], double vecf[])
{
  void bwr_17th(double tin, double *xout, double *vxout, double xt[],
                 double dtmin, char comp_coef, double coef[]);
  char   comp_coef[85];
  static double xx[17], yy[17], zz[17], clock[17];
  register int i, j;
  static int lread1[85], lerr[85];
  static int initl = 0;
  static double tzerol[85];
  static double coefx[85][17], coefy[85][17], coefz[85][17], coefclk[85][17];
  char   orbit_filename[80];
  double tfroms, tnorm, trefno, tzero;
  static double dtmin, dnorm, arcl, half, perday;
  int    ierr=0, itype, ksv, irefep, iread1;
  long   irn;
  static int iorder, iom1, iod2;

  if(initl == 0 )
  {
	  initl = 1;

       for(i=0; i<85; i++)
       { 
		   lread1[i] = -999;
           lerr[i]   = 0;
           tzerol[i] = -999.0;
       }

       iorder = 17;
       dtmin  = ef18_hdr.epoch_interval/60.0;
       dnorm  = 1.0;
       arcl   = (ef18_hdr.nepoch - 1)*dnorm;
       iom1   = iorder - 1;
       iod2   = iorder/2;
       half   = iod2*dnorm;
       perday = 86400.0/ef18_hdr.epoch_interval;
  }

  tfroms = ((mjdt - ef18_hdr.mjd_strt) + (fmjdt - ef18_hdr.fmjd_strt))*perday;

  if(tfroms < (-dnorm))       return (ierr = 1);
  if(tfroms > (arcl + dnorm)) return (ierr = 2);

  itype = 1;
  if(tfroms < half) itype = 2;
  if(tfroms >= (arcl - half)) itype = 3;

  for(ksv=0; ksv<ef18_hdr.number_of_SVs; ksv++)
     { if(ef18_hdr.idsv[ksv] == isv) goto I_KNOW_SAT; }
  return(ierr = 3);

  I_KNOW_SAT:;

  switch(itype) 
  {
     case 1: irefep = tfroms + 1 + dnorm*0.001;
             trefno = irefep - 1;
             tzero  = trefno - iod2;
             iread1 = ksv + (irefep - iod2 - 1)*ef18_hdr.number_of_SVs;
             break;
     case 2: iread1 = ksv;
             tzero  = 0;
             break;
     case 3: iread1 = ksv + (ef18_hdr.nepoch - iom1 - 1)*ef18_hdr.number_of_SVs;
             tzero  = ef18_hdr.nepoch - iorder;
             break;
  }

  comp_coef[ksv] = 't';

  if(lerr[ksv] == 0)
  { 
	  tnorm = tfroms - tzerol[ksv];
      if(abs(iread1-lread1[ksv]) == 0)
      {
		  comp_coef[ksv] = 'f';
          goto NO_NEED_2_READ;
	  }
  }

  tnorm = tfroms - tzero;
  tzerol[ksv] = tzero;
  lread1[ksv] = iread1;

  for(i=0; i<iorder; i++)
  {
	   irn = iread1 + ef18_hdr.number_of_SVs*i;
       fseek(fpef18, (long)(irn*sizeof(ef18)+sizeof(ef18_hdr)), 0);
       fread(&ef18, sizeof(ef18), 1, fpef18);
       xx[i]    = ef18.x/20000.0;
       yy[i]    = ef18.y/20000.0;
       zz[i]    = ef18.z/20000.0;
       clock[i] = ef18.clock/10000.0;
       if(ef18.xyz_good_or_bad == 1) ierr = 4;
       if(ef18.clock_good_or_bad == 1) ierr = 5;
       lerr[ksv] = ierr;
  }

  NO_NEED_2_READ:;

  bwr_17th(tnorm, &recf[0], &vecf[0],    xx, dtmin, comp_coef[ksv], coefx[ksv]);
  bwr_17th(tnorm, &recf[1], &vecf[1],    yy, dtmin, comp_coef[ksv], coefy[ksv]);
  bwr_17th(tnorm, &recf[2], &vecf[2],    zz, dtmin, comp_coef[ksv], coefz[ksv]);
  bwr_17th(tnorm, &recf[3], &vecf[3], clock, dtmin, comp_coef[ksv], coefclk[ksv]);

  return(ierr = 0);

} /* End of function "ef13_bwr17" */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void bwr_17th( double tin, double *fout, double *vfout, double x[],
               double dtmin, char comp_coef, double a[] )
{
  void compute_odd_or_even_coeff( double *alpha, double *a);
  double alpha[8], beta[8], t1, t2, t3, t4, t5, t6, t7, t8, t9, t10,
         t11, t12, t13, t14, t15, t16;
  register int ii;

  if(comp_coef == 't')
  {
	  for(ii=1; ii <=8; ii++)
      {
		  alpha[ii-1] = ( x[8+ii]+x[8-ii]-2.0*x[8] )/(2.0*ii*ii);
          beta[ii-1]  = ( x[8+ii]-x[8-ii] )/(2.0*ii);
	  }
      a[0] = x[8];
      compute_odd_or_even_coeff(alpha, a+9);
      compute_odd_or_even_coeff(beta , a+1);
  }


  t1  = tin-8.0;
  t2  = t1*t1;
  t3  = t2*t1;
  t4  = t2*t2;
  t5  = t4*t1;
  t6  = t3*t3;
  t7  = t4*t3;
  t8  = t4*t4;
  t9  = t4*t5;
  t10 = t5*t5;
  t11 = t5*t6;
  t12 = t6*t6;
  t13 = t7*t6;
  t14 = t7*t7;
  t15 = t7*t8;
  t16 = t8*t8;

  *fout = a[0] + a[1]*t1 + a[9]*t2 + a[2]*t3 + a[10]*t4 + a[3]*t5 + a[11]*t6
          + a[4]*t7 + a[12]*t8 + a[5]*t9 + a[13]*t10 + a[6]*t11 + a[14]*t12
          + a[7]*t13 + a[15]*t14 + a[8]*t15 + a[16]*t16;

  *vfout = a[1] + 2.0*a[9]*t1 + 3.0*a[2]*t2 + 4.0*a[10]*t3 + 5.0*a[3]*t4 +
           6.0*a[11]*t5 + 7.0*a[4]*t6 + 8.0*a[12]*t7 + 9.0*a[5]*t8 +
           10.0*a[13]*t9 + 11.0*a[6]*t10 + 12.0*a[14]*t11 + 13.0*a[7]*t12 +
           14.0*a[15]*t13 + 15.0*a[8]*t14 + 16.0*a[16]*t15;

  *vfout /= dtmin*60.0;
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void compute_odd_or_even_coeff(double p[], double a[])
{
   double b8, b7, b6, b5, b4, b3, b2;
   double c8, c7, c6, c5, c4, c3;
   double d8, d7, d6, d5, d4;
   double e8, e7, e6, e5;
   double f8, f7, f6;
   double g8, g7;

   b8 = (p[7] - p[0])/63.0;
   b7 = (p[6] - p[0])/48.0;
   b6 = (p[5] - p[0])/35.0;
   b5 = (p[4] - p[0])/24.0;
   b4 = (p[3] - p[0])/15.0;
   b3 = (p[2] - p[0])/8.0;
   b2 = (p[1] - p[0])/3.0;

   c8 = (b8 - b2)/60.0;
   c7 = (b7 - b2)/45.0;
   c6 = (b6 - b2)/32.0;
   c5 = (b5 - b2)/21.0;
   c4 = (b4 - b2)/12.0;
   c3 = (b3 - b2)/5.0;

   d8 = (c8 - c3)/55.0;
   d7 = (c7 - c3)/40.0;
   d6 = (c6 - c3)/27.0;
   d5 = (c5 - c3)/16.0;
   d4 = (c4 - c3)/7.0;

   e8 = (d8 - d4)/48.0;
   e7 = (d7 - d4)/33.0;
   e6 = (d6 - d4)/20.0;
   e5 = (d5 - d4)/9.0;

   f8 = (e8 - e5)/39.0;
   f7 = (e7 - e5)/24.0;
   f6 = (e6 - e5)/11.0;

   g8 = (f8 - f6)/28.0;
   g7 = (f7 - f6)/13.0;

   a[7] = (g8-g7)/15.0;
   a[6] = g7   - 140.0*a[7];
   a[5] = f6   - 5278.0*a[7] - 91.0*a[6];
   a[4] = e5   - 61490.0*a[7] - 2002.0*a[6] - 55.0*a[5];
   a[3] = d4   - 196053.0*a[7] - 11440.0*a[6] - 627.0*a[5] - 30.0*a[4];
   a[2] = c3   - 118482.0*a[7] - 13013.0*a[6] - 1408.0*a[5] - 147.0*a[4] - 14.0*a[3];
   a[1] = b2   - 5461.0*a[7] - 1365.0*a[6] - 341.0*a[5] - 85.0*a[4] - 21.0*a[3] - 5.0*a[2];
   a[0] = p[0] - a[7] - a[6] - a[5] - a[4] - a[3] - a[2] - a[1];

}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void compute_end_time(long *yr_end, long *mon_end, long *mday_end, long *hr_end,     long *min_end, double *sec_end)
{
 void mjd_to_ymdhms(long mjd, double fmjd, long *pyear, long *pmonth,
               long *pmday, long *phour, long *pminute, double *psecond);
 void ajmjd(long mjdayo, double fmjdayo, double seconds,
           long *pmjdn, double *pfmjdn);
 long mjd_end;
 double fmjd_end;

ajmjd(ef18_hdr.mjd_strt, ef18_hdr.fmjd_strt,
      ef18_hdr.epoch_interval*(ef18_hdr.nepoch-1), &mjd_end, &fmjd_end);
mjd_to_ymdhms(mjd_end, fmjd_end, yr_end, mon_end, mday_end, hr_end, min_end, sec_end);
if(fabs(60.0-*sec_end) < 1.0e-8)
       round_60_seconds( yr_end, mon_end, mday_end, hr_end, min_end, sec_end);
    
}
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void extract_args (char *fname, int argc, char *argv[])
{
	char drive[MAXDRIVE];
	char directory[MAXDIR];
	char extnt[MAXEXT];

	clrscr();
	fnsplit(argv[0], drive, directory, prog_name, extnt);
	SW_version (prog_name, SW_VERSION, WS_RELEASED);

	printf("                          PROGRAM EF18_SP3.EXE                                 \n");
	printf("                                                                               \n");
	printf("  Program EF18_SP3.EXE converts a DOS binary EF18BIN orbital file to a DOS     \n");
	printf("ASCII SP3ASCII file.  The program begins by providing key header information   \n");
	printf("extracted from the header block of the EF18BIN file.  The filenames EF18BIN    \n");
	printf("and SP3ASCII have been hardcoded.                                              \n");
	printf("                                                                               \n");
	printf("  Using the header information, the user answers the prompts.  The program     \n");
	printf("permits the user to enter the FROM date and the TO date, the epoch spacing,    \n");
	printf("a list of satellites or all (-1 option), and whether he wants just position-   \n");
	printf("clock data (option P) or position-clock plus position-clock rates (option V).  \n");
	delay(INFO_DELAY);
	clrscr();
	printf("                                                                               \n");
	printf("                                                                               \n");
	printf("Contact Mr. Paul Spofford (301)-713-3205 with queries and suggestions.         \n");
	printf("                                                                               \n");
	delay(INFO_DELAY);
	clrscr();
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
SW_version (char *prog_name, char *ver_num, char *rel_date)
/*
 This routine simply generates a display containing a program
   name, version number, and version release date.  Currently,
   the maximum string width (for printing) are as follows:
   prog_name = 8    ver_num = 8    rel_date  = 8
*/
{
  int ii;
  char t_string[TEX_WDTH];
  void center(char *str);

  putchar (LU_ELBOW);
  for (ii=0; ii<TEX_WDTH; ii++) putchar(HORIZONTAL);
  putchar (RU_ELBOW);
  putchar ('\n');

  center ("\0");
  sprintf(t_string, "Program: %-8.8s     Version: %-8.8s     Released: %-8.8s",
          prog_name, ver_num, rel_date);
  center (t_string);
  center ("\0");
  center ("by");
  center ("\0");
  center ("National Geodetic Survey (NOAA)");
  center (" ");
  center ("May 20, 1994  (B.W.Remondi)");
  center ("\0");

  putchar (LL_ELBOW);
  for (ii=0; ii<TEX_WDTH; ii++) putchar(HORIZONTAL);
  putchar (RL_ELBOW);
  putchar ('\n');
  putchar ('\n');
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
void center(char *str)
{
  int i;                        /* Loop control and index variable.   */
  int str_len;                  /* Lenght of string passed            */

  str_len = strlen (str);

  putchar (VERTICAL);
  for (i=1; i < TEX_WDTH/2 - str_len/2; i++) putchar(' ');
  cputs (str);
  for (i += str_len; i <= TEX_WDTH;i++) putchar(' ');
  putchar (VERTICAL);
  putchar ('\n');
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/