jave-ppp.txt

精密单点定位程序

      if(tsec >= allBC[prn][mxrow].gBCjts()){ return mxrow; }  //*** last row

//*** increment until upper point passes the given time

      for(int i=1; i<=mxrow-1; ++i) {
         if(allBC[prn][i+1].gBCjts() >= tsec){
            if(Math.abs(tsec-allBC[prn][i  ].gBCjts())  < 
               Math.abs(tsec-allBC[prn][i+1].gBCjts()))  { return i; }
            else                                         { return i+1; }
         }
      }

//*** fell thru loop -- not allowed

      System.err.println("Loop prob in locBC() " + prn + "  " + tsec);
      System.exit(1);
      return -1;
   }
}
//************************************************************************
class BCstuf {

//*** broadcast orbit object

   private static final double gm    =  3.986005e14;
   private static final double we    =  7.2921151467e-5;        //*** rad/sec
   private static final double gpspi =  3.1415926535898e0;
   private static final double bigf  = -4.442807633e-10;

   private int    id;         //*  my id (just a serial number)
   private double jts;        //*  my time for lookup index
   private int    prn;        //*  prn id

//*** clock parameters

   private double t0c;     //*  clock corrector reference epoch            (s)
   private double af0;     //*  clock offset                               (s)
   private double af1;     //*  clock error rate                          (s/s)
   private double af2;     //*  clock error accelleration              (s/(s*s))

//*** kepler elements

   private double t0e;     //*  reference epoch of ephemeris             (sow)
   private double sqrta;   //*  square root of semimajor axis           sqrt(m)
   private double e;       //*  eccentricity
   private double eye0;    //*  inclination at ref. epoch                  (r)
   private double omg0;    //*  right ascension of asc. node at ref. epoch (r)
   private double w;       //*  argument of perigee                        (r)
   private double em0;     //*  mean anomaly at ref. epoch                 (r)

//*** perturbation parameters

   private double deln;    //*  mean motion difference from computed value (r/s)
   private double eyedot;  //*  rate of change of inclination              (r/s)
   private double omgdot;  //*  rate of change of right ascension          (r/s)

   private double cuc;     //*  cosine harmonic correction to arg. of lat. (r)
   private double cus;     //*  sine harmonic correction to arg. of lat.   (r)
   private double crc;     //*  cosine harmonic correction to radius       (m)
   private double crs;     //*  sine harmonic correction to radius         (m)
   private double cic;     //*  cosine harmonic correction to inclination  (r)
   private double cis;     //*  sine harmonic correction to inclination    (r)

//*** other data 

   private double iode;    //*  issue of data ephemeris
   private double l2cd;    //*  codes on l2 channel (binary, 01-p, 10-c/a)
   private double wkno;    //*  gps week # (to go with toe)               (wk)
   private double l2pd;    //*  l2 p data flag (binary, 0-norm, 1-set off)
   private double svac;    //*  URA index (converted)                      (m)
   private double svhl;    //*  sv health    (6 bits, 0-ok, 1-not)
   private double tgd;     //*  est. group delay   (tsv-tgd)               (s)
   private double iodc;    //*  issue of data clock
   private double tot;     //*  transmission time of message             (sow)

//*****
//***** methods
//*****

   public BCstuf(int myid)   {id=myid;}
   public int    getBC()    {return id;}
   public double gBCjts()   {return jts;}
   public String toString() {return "" + id;}

   public int    gBCprn()   {return prn;}
   public double gBCt0c()   {return t0c;}
   public double gBCt0e()   {return t0e;}
   public double gBCwkno()  {return wkno;}
   public double gBCsvac()  {return svac;}
   public double gBCsvhl()  {return svhl;}
   public double gBCtot()   {return tot;}

//*****
//***** load methods
//*****

   public void put0BC(String myString, int myid, Gtime tStuf) {

//*** load record (epoch and clock)
//*** y2k code fails after 2079

//*** clock parameters

      String subField;
      int     iyr,imo,idy,ihr,imn;
      double  sec,tsec;

      subField = myString.substring( 0,  2);
      prn      = Integer.parseInt(subField.trim());

      subField = myString.substring( 2,  5);
      iyr      = Integer.parseInt(subField.trim());
      if     (iyr <  80) {iyr=iyr+2000;}               //*** y2k
      else if(iyr <= 99) {iyr=iyr+1900;}               //*** fail 2080

      subField = myString.substring( 5,  8);
      imo      = Integer.parseInt(subField.trim());
      subField = myString.substring( 8, 11);
      idy      = Integer.parseInt(subField.trim());
      subField = myString.substring(11, 14);
      ihr      = Integer.parseInt(subField.trim());
      subField = myString.substring(14, 17);
      imn      = Integer.parseInt(subField.trim());
      subField = myString.substring(17, 22);
      sec      = Double.valueOf(subField.trim()).doubleValue();

//*** if first epoch, check to establish time offset

      if(myid <= 1){
         if(tStuf.getMjd0() <= -1){
            System.err.println("Time was set in nav file -- check obs file");
            tStuf.setjd0(iyr, imo, idy);
         }
      }

      jts = tStuf.civjts(iyr,imo,idy,ihr,imn,sec);   //*** also used for index
      t0c = tStuf.jtssow(jts);

      subField = myString.substring(22, 41).replace('D','e');
      af0      = Double.valueOf(subField.trim()).doubleValue();
      subField = myString.substring(41, 60).replace('D','e');
      af1      = Double.valueOf(subField.trim()).doubleValue();
      subField = myString.substring(60, 79).replace('D','e');
      af2      = Double.valueOf(subField.trim()).doubleValue();
   }

   public void put1BC(String myString) {

//*** load record (BC Orbit 1)

      String subField;

      subField = myString.substring( 3, 22).replace('D','e');
      iode     = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(22, 41).replace('D','e');
      crs      = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(41, 60).replace('D','e');
      deln     = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(60, 79).replace('D','e');
      em0      = Double.valueOf(subField.trim()).doubleValue();
   }

   public void put2BC(String myString) {

//*** load record (BC Orbit 2)

      String subField;

      subField = myString.substring( 3, 22).replace('D','e');
      cuc      = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(22, 41).replace('D','e');
      e        = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(41, 60).replace('D','e');
      cus      = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(60, 79).replace('D','e');
      sqrta    = Double.valueOf(subField.trim()).doubleValue();
   }

   public void put3BC(String myString) {

//*** load record (BC Orbit 3)

      String subField;

      subField = myString.substring( 3, 22).replace('D','e');
      t0e      = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(22, 41).replace('D','e');
      cic      = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(41, 60).replace('D','e');
      omg0     = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(60, 79).replace('D','e');
      cis      = Double.valueOf(subField.trim()).doubleValue();
   }

   public void put4BC(String myString) {

//*** load record (BC Orbit 4)

      String subField;

      subField = myString.substring( 3, 22).replace('D','e');
      eye0     = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(22, 41).replace('D','e');
      crc      = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(41, 60).replace('D','e');
      w        = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(60, 79).replace('D','e');
      omgdot   = Double.valueOf(subField.trim()).doubleValue();
   }

   public void put5BC(String myString) {

//*** load record (BC Orbit 5)

      String subField;

      subField = myString.substring( 3, 22).replace('D','e');
      eyedot   = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(22, 41).replace('D','e');
      l2cd     = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(41, 60).replace('D','e');
      wkno     = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(60, 79).replace('D','e');
      l2pd    = Double.valueOf(subField.trim()).doubleValue();
   }

   public void put6BC(String myString) {

//*** load record (BC Orbit 6)

      String subField;

      subField = myString.substring( 3, 22).replace('D','e');
//*** troublesome field in gait nav file
   try {
      svac     = Double.valueOf(subField.trim()).doubleValue();
   } 
   catch (NumberFormatException e) {
//****System.err.println("Reset acc. for " + subField.trim());
      svac = 384.0;
   }

      subField = myString.substring(22, 41).replace('D','e');
      svhl     = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(41, 60).replace('D','e');
      tgd      = Double.valueOf(subField.trim()).doubleValue();

      subField = myString.substring(60, 79).replace('D','e');
      iodc     = Double.valueOf(subField.trim()).doubleValue();
   }

   public void put7BC(String myString) {

//*** load record (BC Orbit 7)

      String subField;

      subField = myString.substring( 3, 22).replace('D','e');
      tot      = Double.valueOf(subField.trim()).doubleValue();
   }

//*****
//***** compute methods
//*****

   public PosT bcxyzt (double tsv, Modes myModes) {
      
//*** compute pos of broadcast orbit (BC) for a single time

//*** sense of the corrector (with relativity) is:   t=tsv-dt-dtr
//*** BC is in GPS system time (use clock corr, time in sec.of.week)

      double dt, dtr, bcclok;

//*** correct satellite time

      dt     = bcclok1(tsv);                       //*** icd-200 20.3.3.3.3.1
      dtr    = bcclok2(tsv,dt);
      bcclok = dt+dtr;

//*** icd-200 20.3.3.3.3.2  (not dual and not b-cast iono)

      if( (!myModes.isdFreq()) && (!myModes.isbIon())) bcclok = bcclok - tgd;

//*** compute broadcast position at corrected time

      return bcxyz(tsv-bcclok, bcclok, dtr);
   }

   public double bcclok1 (double t) {
      
//*** compute satellite clock correction (non-relativity term)
//*** sense of the corrector is:   t=tsv-dt

      double        dt;

      if( (t-t0c) > 302400.0 ) t=t-604800.0;        //*** handle week rollover
      if( (t-t0c) <-302400.0 ) t=t+604800.0;        //*** valid icd200

      dt=af0+af1*(t-t0c)+af2*(t-t0c)*(t-t0c);       //*** valid icd200
      return dt;
   }

   public double bcclok2 (double t, double dt) {
      
//*** compute satellite clock correction (relativity term)
//*** sense of the corrector is:   t=tsv-dt-dtr

      double dtr;       //*** relativity term
      double ek;

      if( (t-t0c) > 302400.0 ) t=t-604800.0;        //*** handle week rollover
      if( (t-t0c) <-302400.0 ) t=t+604800.0;        //*** valid icd200
      
//*** relativity term

      ek=eccano(t-dt);
      dtr=bigf*e*sqrta*Math.sin(ek);           //*** valid icd200
      return dtr;
   }

   public PosT bcxyz(double t, double bcclok, double dtr) {

//*** compute ECBF XYZ at time t from broadcast elements
//***               units: seconds, meters, radians
//*** note: t corrected for SV clock error, also passed in bcclok
//*** note: also passing relativitistic correction

      double a,en0,tk,en,emk,ek,sinek,cosek,f,u,cu2,su2,duk,drk,dik;
      double uk,rk,eyek,ceyek,xpk,ypk,omgk,comgk,somgk,xk,yk,zk;

      if( (t-t0e) > 302400.0 ) t=t-604800.0;            //*** handle week rollover
      if( (t-t0e) <-302400.0 ) t=t+604800.0;            //*** valid icd200

//*** mean motion

      a   = sqrta*sqrta;
      en0 = Math.sqrt(gm/(a*a*a));

//*** time since orbit reference epoch

      tk = t-t0e;

//*** corrected mean motion

      en = en0+deln;

//*** mean anomaly, M

      emk = em0+en*tk;

//*** solve kepler's equation for eccentric anomaly, E

      ek    = kepEQ(emk);
      sinek = Math.sin(ek);
      cosek = Math.cos(ek);

//*** true anomaly, f

      f = Math.atan2( Math.sqrt(1.0-e*e)*sinek, cosek-e);

//*** argument of latitude, u

      u   = f+w;
      cu2 = Math.cos(u+u);
      su2 = Math.sin(u+u);

//*** corrections to the arg. of lat., radius, and inclination

      duk = cuc*cu2+cus*su2;
      drk = crc*cu2+crs*su2;
      dik = cic*cu2+cis*su2;

//*** correct the arg. of lat., radius, and inclination

      uk    = u               + duk;
      rk    = a*(1.0-e*cosek) + drk;
      eyek  = eye0+eyedot*tk  + dik;
      ceyek = Math.cos(eyek);

//*** position in the orbital plane

      xpk = rk*Math.cos(uk);
      ypk = rk*Math.sin(uk);

//*** correct the longitude of the ascending node

      omgk  = omg0+(omgdot-we)*tk-we*t0e;
      comgk = Math.cos(omgk);
      somgk = Math.sin(omgk);