orbit.Computations (Java2HTML)

package orbit;

/**
 * <p>Title: ISS Transit</p>
 * <p>Description: Predict transits of the ISS across the moon and sun</p>
 * <p>Copyright: Copyright (c) 2002</p>
 * <p>Company: J2EE Consultants</p>
 * @author Thomas Fly
 * @version 1.0
 *
 * @see also http://cow.physics.wisc.edu/~craigm/idl/down/hprstatn.txt
 */

public class Computations {
/*
  Values for EARTH_RADIUS, EQUATORIAL_RADIUS, CK2, CK4, xke, and qoms2t are given on page 76 of
  http://www.celestrak.com/NORAD/documentation/spacetrk.pdf

where er = Earth radii


*/
  static final double EARTH_RADIUS = 6378.135;  // kilometers
  // This is kind of odd- it's used basically to indicate that distances are measured in units of EARTH_RADIUS_KM:
  static final double EQUATORIAL_RADIUS = 1;

//  static final double J2 = 1.0826158e-3;  /* J2 harmonic (WGS '72) */
//  static final double J3 = -2.53881e-6;   /* J3 harmonic (WGS '72) */
//  static final double J4 = -1.65597e-6;   /* J4 harmonic (WGS '72) */

  // from http://celestrak.com/NORAD/documentation/spacetrk.pdf
  static final double CK2 = 5.413080e-4;    // equals J2 / 2 * EQUATORIAL_RADIUS^2
  static final double CK4 = .62098875e-6;   // equals -3 * J4 / 8 * EQUATORIAL_RADIUS^4

  static final double S = EQUATORIAL_RADIUS + 78 / EARTH_RADIUS;

/*
  from SGP4SDP4.PAS:
  Procedure Define_Derived_Constants;
    begin
    xke    := Sqrt(3600*ge/Cube(xkmper));   {Sqrt(ge) ER^3/min^2}
    qoms2t := Sqr(Sqr(qo-s));               {(qo-s)^4 ER^4}
    end; {Define_Derived_Constants}

  from SGP_INTF.PAS:
  "ke" = Math.sqrt(G/M) where G is Newton's universal gravitational constant and M is the mass of the Earth
         ge = 398600.8; {Earth gravitational constant (WGS '72)}

  ae       = 1;
  xkmper   = 6378.135;        {Earth equatorial radius - kilometers (WGS '72)}
  qo       = ae + 120/xkmper;
  s        = ae + 78/xkmper;

  so... qo - s = 42 / EARTH_RADIUS = 6.584997e-3
  qomss2t is evidently "qo minus s squared 2 times"
*/
  static double xke = .743669161e-1;
  static double qoms2t = 1.88027916e-9;
  static double xj3 = -.253881e-5;

/*  Specified by Two-Line Elements:
  static double dragCoefficient = 0;

  static double rightAscensionAscendingNode = 0;
  static double orbitInclination = 0;
  static double orbitEccentricity = 0;
  static double orbitPerigee = 0;
  static double meanAnomaly = 0;
  static double meanMotion = 0;
*/

  static double xmdot;
  static double omgdot;
  static double xnodot;
  static double xnodcf;
  static double c1;
  static double c4;
  static double t2cof;
  static double omgcof;
  static double xmcof;
  static double eta;
  static double delmo;
  static double c5;
  static double sinmo;
  static double semimajorAxis;
  static double meanMotionAtEpoch;
  static double xlcof;
  static double aycof;
  static double temp2;
  static double temp3;
  static double temp1;

  static double x3thm1;
  static double x1mth2;
  static double x7thm1;
  static double cosio;
  static double sinio;

  static boolean perigeeLessThan220km;

  static double d2 = 0;
  static double d3 = 0;
  static double d4 = 0;

  static double t3cof = 0;
  static double t4cof = 0;
  static double t5cof = 0;

  public Computations(TwoLineElement tle) {
/*
    Equation 1.

*/
    double a1 = Math.pow(xke / tle.meanMotion, 2.0 / 3);

/*
    Equation 2.

*/
    cosio = Math.cos(tle.orbitInclination);
    double theta2 = cosio * cosio;
    x3thm1 = 3 * theta2 - 1;

    double betao2 = 1 - tle.orbitEccentricity * tle.orbitEccentricity;

    double betao = Math.sqrt(betao2);


    double del1 = 1.5 * CK2 * x3thm1 / (a1 * a1 * betao * betao2);
/*
    Equation 3.

*/
    double ao = a1 * (1 - del1 * (1.0 / 3 + del1 * (1 + 134.0 / 81 * del1)));
/*
    Equation 4 (nearly identical to Equation 2).

*/
    double delo = 1.5 * CK2 * x3thm1 / (ao * ao * betao * betao2);
/*
    Equation 5.

*/
    meanMotionAtEpoch = tle.meanMotion / (1 + delo);
/*
    Equation 6.

*/
    semimajorAxis = ao / (1 - delo);
/*
    { Initialization }
    { For perigee less than 220 kilometers, the isimp flag is set and
      the equations are truncated to linear variation in sqrt a and
      quadratic variation in mean anomaly.  Also, the c3 term, the
      delta omega term, and the delta m term are dropped. }

    isimp := 0;
*/
    perigeeLessThan220km = false;
/*
    IF((AODP*(1.-EO)/AE) .LT. (220./XKMPER+AE)) ISIMP=1
    if (aodp*(1 - eo)/ae) < (220/xkmper + ae) then
      isimp := 1;
*/
    if ((semimajorAxis * (1 - tle.orbitEccentricity) / EQUATORIAL_RADIUS) < (220.0 / EARTH_RADIUS + EQUATORIAL_RADIUS))
      perigeeLessThan220km = true;
/*
    * FOR PERIGEE BELOW 156 KM, THE VALUES OF
    * S AND QOMS2T ARE ALTERED
    { For perigee below 156 km, the values of s and qoms2t are altered. }

    S4=S
    s4 := s;
*/
    double s4 = S;
/*
    QOMS24=QOMS2T
    qoms24 := qoms2t;
*/
    double qoms24 = qoms2t;
/*
    PERIGE=(AODP*(1.-EO)-AE)*XKMPER
    perige := (aodp*(1 - eo) - ae)*xkmper;
*/
    double perigee = (semimajorAxis * (1 - tle.orbitEccentricity) - EQUATORIAL_RADIUS) * EARTH_RADIUS;
/*
    IF(PERIGE .GE. 156.) GO TO 10
    S4=PERIGE-78.
    IF(PERIGE .GT. 98.) GO TO 9
    S4=20.
  9 QOMS24=((120.-S4)*AE/XKMPER)**4
    S4=S4/XKMPER+AE
 10 PINVSQ=1./(AODP*AODP*BETAO2*BETAO2)

    if (perige >= 156) then
      goto 10;
    s4 := perige - 78;
    if (perige > 98) then
      goto 9;
    s4 := 20;
    9:
    qoms24 := Power((120 - s4)*ae/xkmper,4);
    s4 := s4/xkmper + ae;
    10:
    pinvsq := 1/(aodp*aodp*betao2*betao2);



*/
    if (perigee < 156) {
      s4 = perigee - 78;
      if (perigee <= 98) {
        s4 = 20;
      }
      else {
        qoms24 = Math.pow((120 - s4) * EQUATORIAL_RADIUS / EARTH_RADIUS, 4);
        s4 = s4 / EARTH_RADIUS + EQUATORIAL_RADIUS;
      }
    }
/*
    TSI=1./(AODP-S4)
    tsi := 1/(aodp - s4);

*/
    double tsi = 1 / (semimajorAxis - s4);
/*
    ETA=AODP*EO*TSI
    eta := aodp*eo*tsi;

*/
    eta = semimajorAxis * tle.orbitEccentricity * tsi;
/*
    ETASQ=ETA*ETA
    etasq := eta*eta;
*/
    double etasq = eta * eta;
/*
    EETA=EO*ETA
    eeta := eo*eta;
*/
    double eeta = tle.orbitEccentricity * eta;
/*
    PSISQ=ABS(1.-ETASQ)
    psisq := abs(1 - etasq);
*/
    double psisq = Math.abs(1 - etasq);
/*
    COEF=QOMS24*TSI**4
    coef := qoms24*Power(tsi,4);
*/
    double coef = qoms24 * Math.pow(tsi, 4);
/*
    COEF1=COEF/PSISQ**3.5
    coef1 := coef/Power(psisq,3.5);
*/
    double coef1 = coef / Math.pow(psisq, 3.5);
/*
    C2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*
  1 CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)))

  c2 := coef1*xnodp*(aodp*(1 + 1.5*etasq + eeta*(4 + etasq))
      + 0.75*ck2*tsi/psisq*x3thm1*(8 + 3*etasq*(8 + etasq)));

*/
    double c2 = coef1 * meanMotionAtEpoch * (semimajorAxis * (1 + 1.5 * etasq + eeta * (4 + etasq)) +
                                0.75 * CK2 * tsi / psisq * x3thm1 * (8 + 3 * etasq * (8 + etasq)));
/*
    C1=BSTAR*C2
    c1 := bstar*c2;

*/
    c1 = tle.dragCoefficient * c2;
/*
    SINIO=SIN(XINCL)
    sinio := Sin(xincl);
*/
    sinio = Math.sin(tle.orbitInclination);
/*
    A3OVK2=-XJ3/CK2*AE**3
    a3ovk2 := -xj3/ck2*Power(ae,3);
*/
    double a3ovk2 = - xj3 / CK2 * Math.pow(EQUATORIAL_RADIUS, 3);
/*
    C3=COEF*TSI*A3OVK2*XNODP*AE*SINIO/EO
    c3 := coef*tsi*a3ovk2*xnodp*ae*sinio/eo;

*/
    double c3 = coef * tsi * a3ovk2 * meanMotionAtEpoch * EQUATORIAL_RADIUS * sinio / tle.orbitEccentricity;
/*
    X1MTH2=1.-THETA2
    x1mth2 := 1 - theta2;
*/
    x1mth2 = 1 - theta2;
/*
    C4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*
  1 (2.+.5*ETASQ)+EO*(.5+2.*ETASQ)-2.*CK2*TSI/
  2 (AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*
  3 (1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*
  4 (1.+ETASQ))*COS(2.*OMEGAO)))

  c4 := 2*xnodp*coef1*aodp*betao2*(eta*(2 + 0.5*etasq)
      + eo*(0.5 + 2*etasq) - 2*ck2*tsi/(aodp*psisq)
      *(-3*x3thm1*(1 - 2*eeta + etasq*(1.5 - 0.5*eeta))
      + 0.75*x1mth2*(2*etasq - eeta*(1 + etasq))*Cos(2*omegao)));


*/
    c4 = 2 * meanMotionAtEpoch *coef1 * semimajorAxis * betao2 * (eta * (2 + 0.5 * etasq) +
      tle.orbitEccentricity * (0.5 + 2 * etasq) - 2 * CK2 * tsi / (semimajorAxis * psisq) *
      (-3 * x3thm1 * (1 - 2 * eeta + etasq * (1.5 - 0.5 * eeta)) +
      0.75 * x1mth2 * (2 * etasq - eeta * (1 + etasq)) * Math.cos(2 * tle.orbitPerigee)));
/*
    C5=2.*COEF1*AODP*BETAO2*(1.+2.75*(ETASQ+EETA)+EETA*ETASQ)
    c5 := 2*coef1*aodp*betao2*(1 + 2.75*(etasq + eeta) + eeta*etasq);

*/
    c5 = 2 * coef1 * semimajorAxis * betao2 * (1 + 2.75 * (etasq + eeta) + eeta * etasq);
/*
    THETA4=THETA2*THETA2
    theta4 := theta2*theta2;
*/
    double theta4 = theta2 * theta2;
/*
    10 PINVSQ=1./(AODP*AODP*BETAO2*BETAO2)
    pinvsq := 1/(aodp*aodp*betao2*betao2);

    TEMP1=3.*CK2*PINVSQ*XNODP
    temp1 := 3*ck2*pinvsq*xnodp;
*/
    double pinvsq = 1 / (semimajorAxis * semimajorAxis * betao2 * betao2);
    temp1 = 3 * CK2 * pinvsq * meanMotionAtEpoch;
/*
    TEMP2=TEMP1*CK2*PINVSQ
    temp2 := temp1*ck2*pinvsq;
*/
    temp2 = temp1 * CK2 * pinvsq;
/*
    TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP
    temp3 := 1.25*ck4*pinvsq*pinvsq*xnodp;
*/
    temp3 = 1.25 * CK4 * pinvsq * pinvsq * meanMotionAtEpoch;
/*
    XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*
  1 (13.-78.*THETA2+137.*THETA4)

  xmdot := xnodp + 0.5*temp1*betao*x3thm1
         + 0.0625*temp2*betao*(13 - 78*theta2 + 137*theta4);
*/
    xmdot = meanMotionAtEpoch + 0.5 * temp1 * betao *x3thm1 +
      0.0625 * temp2 * betao * (13 - 78 * theta2 + 137 * theta4);
/*
    X1M5TH=1.-5.*THETA2
    x1m5th := 1 - 5*theta2;
*/
    double x1m5th = 1 - 5 * theta2;
/*
    OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+
  1 395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4)

  omgdot := -0.5*temp1*x1m5th + 0.0625*temp2*(7 - 114*theta2 +395*theta4)
          + temp3*(3 - 36*theta2 + 49*theta4);

*/
    omgdot = -0.5 * temp1 * x1m5th + 0.0625 * temp2 * (7 - 114 * theta2 + 395 * theta4) +
      temp3 * (3 - 36 * theta2 + 49 * theta4);
/*
    XHDOT1=-TEMP1*COSIO
    xhdot1 := -temp1*cosio;
*/
    double xhdot1 = -temp1 * cosio;
/*
    XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-
  1 7.*THETA2))*COSIO

  xnodot := xhdot1 + (0.5*temp2*(4 - 19*theta2)
          + 2*temp3*(3 - 7*theta2))*cosio;

*/
    xnodot = xhdot1 + (0.5 * temp2 * (4 - 19 * theta2) + 2 * temp3 * (3 - 7 * theta2)) * cosio;
/*
    OMGCOF=BSTAR*C3*COS(OMEGAO)
    omgcof := bstar*c3*Cos(omegao);

*/
    omgcof = tle.dragCoefficient * c3 * Math.cos(tle.orbitPerigee);
/*
    XMCOF=-TOTHRD*COEF*BSTAR*AE/EETA
    xmcof := -tothrd*coef*bstar*ae/eeta;
*/
    xmcof = -2.0 / 3 * coef * tle.dragCoefficient * EQUATORIAL_RADIUS / eeta;
/*
    XNODCF=3.5*BETAO2*XHDOT1*C1
    xnodcf := 3.5*betao2*xhdot1*c1;
*/
    xnodcf = 3.5 * betao2 * xhdot1 * c1;
/*
    T2COF=1.5*C1
    t2cof := 1.5*c1;
*/
    t2cof = 1.5 * c1;
/*
    XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO)
    xlcof := 0.125*a3ovk2*sinio*(3 + 5*cosio)/(1 + cosio);

*/
    xlcof = 0.125 * a3ovk2 * sinio * (3 + 5 * cosio) / (1 + cosio);
/*
    AYCOF=.25*A3OVK2*SINIO
    aycof := 0.25*a3ovk2*sinio;
*/
    aycof = 0.25 * a3ovk2 * sinio;
/*
    DELMO=(1.+ETA*COS(XMO))**3
    delmo := Power(1 + eta*Cos(xmo),3);
*/
    delmo = Math.pow(1 + eta * Math.cos(tle.meanAnomaly), 3);
/*
    SINMO=SIN(XMO)
    sinmo := Sin(xmo);
*/
    sinmo = Math.sin(tle.meanAnomaly);
/*
    X7THM1=7.*THETA2-1.
    x7thm1 := 7*theta2 - 1;
*/
    x7thm1 = 7 * theta2 - 1;
/*
    IF(ISIMP .EQ. 1) GO TO 90

    if (isimp = 1) then
      goto 90;
*/
    if (!perigeeLessThan220km) {
/*
      C1SQ=C1*C1
      c1sq := c1*c1;
*/
      double c1sq = c1 * c1;
/*
      D2=4.*AODP*TSI*C1SQ
      d2 := 4*aodp*tsi*c1sq;

*/
      d2 = 4 * semimajorAxis * tsi * c1sq;
/*
      TEMP=D2*TSI*C1/3.
      temp := d2*tsi*c1/3;
*/
      double temp = d2 * tsi * c1 / 3;
/*
      D3=(17.*AODP+S4)*TEMP
      d3 := (17*aodp + s4)*temp;

*/
      d3 = (17 * semimajorAxis + s4) * temp;
/*
      D4=.5*TEMP*AODP*TSI*(221.*AODP+31.*S4)*C1
      d4 := 0.5*temp*aodp*tsi*(221*aodp + 31*s4)*c1;

*/
      d4 = 0.5 * temp * semimajorAxis * tsi * (221 * semimajorAxis + 31 * s4) * c1;
/*
      T3COF=D2+2.*C1SQ
      t3cof := d2 + 2*c1sq;
*/
      t3cof = d2 + 2 * c1sq;
/*
      T4COF=.25*(3.*D3+C1*(12.*D2+10.*C1SQ))
      t4cof := 0.25*(3*d3 + c1*(12*d2 + 10*c1sq));
*/
      t4cof = 0.25 * (3 * d3 + c1 * (12 * d2 + 10 * c1sq));
/*
      T5COF=.2*(3.*D4+12.*C1*D3+6.*D2*D2+15.*C1SQ*(
    1 2.*D2+C1SQ))

      t5cof := 0.2*(3*d4 + 12*c1*d3 + 6*d2*d2 + 15*c1sq*(2*d2 + c1sq));
*/
      t5cof = 0.2 * (3 * d4 + 12 * c1 * d3 + 6 * d2 * d2 + 15 * c1sq * (2 * d2 + c1sq));
    }

/* "90"
    IFLAG=0
    iflag := 0;
*/
//    double iflag = 0;
  }

  public static void sgp4(double elapsedTime, TwoLineElement tle) {
/*  "100"
    * UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
    { Update for secular gravity and atmospheric drag. }
    XMDF=XMO+XMDOT*TSINCE
    xmdf := xmo + xmdot*tsince;

*/
    double xmdf = tle.meanAnomaly + xmdot * elapsedTime;
/*
    OMGADF=OMEGAO+OMGDOT*TSINCE
    omgadf := omegao + omgdot*tsince;

*/
    double omgadf = tle.orbitPerigee + omgdot * elapsedTime;
/*
    XNODDF=XNODEO+XNODOT*TSINCE
    xnoddf := xnodeo + xnodot*tsince;

*/
    double xnoddf = tle.rightAscensionAscendingNode + xnodot * elapsedTime;
/*
    OMEGA=OMGADF
    omega := omgadf;
*/
    double omega = omgadf;
/*
    XMP=XMDF
    xmp := xmdf;
*/
    double xmp = xmdf;
/*
    TSQ=TSINCE*TSINCE
    tsq := tsince*tsince;
*/
    double tsq = elapsedTime * elapsedTime;
/*
    XNODE=XNODDF+XNODCF*TSQ
    xnode := xnoddf + xnodcf*tsq;
*/
    double xnode = xnoddf + xnodcf * tsq;
/*
    TEMPA=1.-C1*TSINCE
    tempa := 1 - c1*tsince;
*/
    double tempa = 1 - c1 * elapsedTime;
/*
    TEMPE=BSTAR*C4*TSINCE
    tempe := bstar*c4*tsince;

*/
    double tempe = tle.dragCoefficient * c4 * elapsedTime;
/*
    TEMPL=T2COF*TSQ
    templ := t2cof*tsq;
*/
    double templ = t2cof * tsq;
/*
    IF(ISIMP .EQ. 1) GO TO 110

    if (isimp = 1) then
      goto 110;
*/
    if (!perigeeLessThan220km) {
/*
      DELOMG=OMGCOF*TSINCE
      delomg := omgcof*tsince;

*/
      double delomg = omgcof * elapsedTime;
/*
      DELM=XMCOF*((1.+ETA*COS(XMDF))**3-DELMO)
      delm := xmcof*(Power(1 + eta*Cos(xmdf),3) - delmo);

*/
      double delm = xmcof * (Math.pow(1 + eta * Math.cos(xmdf), 3) - delmo);
/*
      TEMP=DELOMG+DELM
      temp := delomg + delm;
*/
      double temp = delomg + delm;
/*
      XMP=XMDF+TEMP
      xmp := xmdf + temp;

*/
      xmp = xmdf + temp;
/*
      OMEGA=OMGADF-TEMP
      omega := omgadf - temp;

*/
      omega = omgadf - temp;
/*
      TCUBE=TSQ*TSINCE
      tcube := tsq*tsince;

*/
      double tcube = tsq * elapsedTime;
/*
      TFOUR=TSINCE*TCUBE
      tfour := tsince*tcube;
*/
      double tfour = elapsedTime * tcube;
/*
      TEMPA=TEMPA-D2*TSQ-D3*TCUBE-D4*TFOUR
      tempa := tempa - d2*tsq - d3*tcube - d4*tfour;
*/
      tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour;
/*
      TEMPE=TEMPE+BSTAR*C5*(SIN(XMP)-SINMO)
      tempe := tempe + bstar*c5*(Sin(xmp) - sinmo);

*/
      tempe = tempe + tle.dragCoefficient * c5 * (Math.sin(xmp) - sinmo);
/*
      TEMPL=TEMPL+T3COF*TCUBE+
    1 TFOUR*(T4COF+TSINCE*T5COF)

      templ := templ + t3cof*tcube + tfour*(t4cof + tsince*t5cof);
*/
      templ = templ + t3cof * tcube + tfour * (t4cof + elapsedTime * t5cof);

    }   /* End if (!perigeeLessThan220km) */

/*  "110"
    A=AODP*TEMPA**2
    a := aodp*Sqr(tempa);

*/
    double a = semimajorAxis * Math.pow(tempa, 2);
/*
    E=EO-TEMPE
    e := eo - tempe;

*/
    double e = tle.orbitEccentricity - tempe;
/*
    XL=XMP+OMEGA+XNODE+XNODP*TEMPL
    xl := xmp + omega + xnode + xnodp*templ;

*/
    double xl = xmp + omega + xnode + meanMotionAtEpoch * templ;
/*
    BETA=SQRT(1.-E*E)
    beta := sqrt(1 - e*e);

*/
    double beta = Math.sqrt(1 - e * e);
/*
    XN=XKE/A**1.5
    xn := xke/Power(a,1.5);

*/
    double xn = xke / Math.pow(a, 1.5);
/*
    * LONG PERIOD PERIODICS
    { Long period periodics }
    AXN=E*COS(OMEGA)
    axn := e*Cos(omega);

*/
    double axn = e * Math.cos(omega);
/*
    TEMP=1./(A*BETA*BETA)
    temp := 1/(a*beta*beta);
*/
    double temp = 1 / (a * beta * beta);
/*
    XLL=TEMP*XLCOF*AXN
    xll := temp*xlcof*axn;

*/
    double xll = temp * xlcof * axn;
/*
    AYNL=TEMP*AYCOF
    aynl := temp*aycof;

*/
    double aynl = temp * aycof;
/*
    XLT=XL+XLL
    xlt := xl + xll;

*/
    double xlt = xl + xll;
/*
    AYN=E*SIN(OMEGA)+AYNL
    ayn := e*Sin(omega) + aynl;

*/
    double ayn = e * Math.sin(omega) + aynl;
/*
    * SOLVE KEPLERS EQUATION
    { Solve Kepler's Equation }
    CAPU=FMOD2P(XLT-XNODE)
    capu := fmod2p(xlt - xnode);

*/
    double capu = (xlt - xnode) % (2 * Math.PI);
/*
    TEMP2=CAPU
    temp2 := capu;

*/
    temp2 = capu;
/*
    DO 130 I=1,10

    for i := 1 to 10 do
      begin
*/
    double temp4 = 0;
    double temp5 = 0;
    double temp6 = 0;

    double cosepw = 0;
    double sinepw = 0;

    for (int i = 1;  i <= 10;  i++) {
/*
      SINEPW=SIN(TEMP2)
      sinepw := Sin(temp2);
*/
      sinepw = Math.sin(temp2);
/*
      COSEPW=COS(TEMP2)
      cosepw := Cos(temp2);
*/
      cosepw = Math.cos(temp2);
/*
      TEMP3=AXN*SINEPW
      temp3 := axn * sinepw;
*/
      temp3 = axn * sinepw;
/*
      TEMP4=AYN*COSEPW
      temp4 := ayn*cosepw;
*/
      temp4 = ayn * cosepw;
/*
      TEMP5=AXN*COSEPW
      temp5 := axn*cosepw;
*/
      temp5 = axn * cosepw;
/*
      TEMP6=AYN*SINEPW
      temp6 := ayn*sinepw;
*/
      temp6 = ayn * sinepw;
/*
      EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2
      epw := (capu - temp4 + temp3 - temp2)/(1 - temp5 - temp6) + temp2;



*/
      double epw = (capu - temp4 + temp3 - temp2) / (1 - temp5 - temp6) + temp2;
/*
      BUTT UGLY!!!

      IF(ABS(EPW-TEMP2) .LE. E6A) GO TO 140

      if (abs(epw - temp2) <= e6a) then
        goto 140;
*/
      if (Math.abs(epw - temp2) <= 1e-6)
        break;
/*  "130"
    TEMP2=EPW

    temp2 := epw;
    end; {for i}
*/
      temp2 = epw;
    }

/*  "140"
    * SHORT PERIOD PRELIMINARY QUANTITIES
    { Short period preliminary quantities }
    ECOSE=TEMP5+TEMP6
    ecose := temp5 + temp6;

*/
    double ecose = temp5 + temp6;
/*
    ESINE=TEMP3-TEMP4
    esine := temp3 - temp4;

*/
    double esine = temp3 - temp4;
/*
    ELSQ=AXN*AXN+AYN*AYN
    elsq := axn*axn + ayn*ayn;

*/
    double elsq = axn * axn + ayn * ayn;
/*
    TEMP=1.-ELSQ
    temp := 1 - elsq;
*/
    temp = 1 - elsq;
/*
    PL=A*TEMP
    pl := a*temp;

*/
    double pl = a * temp;
/*
    R=A*(1.-ECOSE)
    r := a*(1 - ecose);

*/
    double r = a * (1 - ecose);
/*
    TEMP1=1./R
    temp1 := 1/r;
*/
    temp1 = 1D / r;
/*
    RDOT=XKE*SQRT(A)*ESINE*TEMP1
    rdot := xke*sqrt(a)*esine*temp1;

*/
    double rdot = xke * Math.sqrt(a) * esine * temp1;
/*
    RFDOT=XKE*SQRT(PL)*TEMP1
    rfdot := xke*sqrt(pl)*temp1;

*/
    double rfdot = xke * Math.sqrt(pl) * temp1;
/*
    TEMP2=A*TEMP1
    temp2 := a*temp1;
*/
    temp2 = a * temp1;
/*
    BETAL=SQRT(TEMP)
    betal := sqrt(temp);
*/
    double betal = Math.sqrt(temp);
/*
    TEMP3=1./(1.+BETAL)
    temp3 := 1/(1 + betal);
*/
    temp3 = 1 / (1 + betal);
/*
    COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3)
    cosu := temp2*(cosepw - axn + ayn*esine*temp3);

*/
    double cosu = temp2 * (cosepw - axn + ayn * esine * temp3);
/*
    SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3)
    sinu := temp2*(sinepw - ayn - axn*esine*temp3);

*/
    double sinu = temp2 * (sinepw - ayn - axn * esine * temp3);
/*
    U=ACTAN(SINU,COSU)
    u := actan(sinu,cosu);

*/
    double u = Math.atan2(sinu, cosu);  // ... MAYBE
/*
    SIN2U=2.*SINU*COSU
    sin2u := 2*sinu*cosu;
*/
    double sin2u = 2 * sinu * cosu;
/*
    COS2U=2.*COSU*COSU-1.
    cos2u := 2*cosu*cosu - 1;
*/
    double cos2u = 2 * cosu * cosu - 1;
/*
    TEMP=1./PL
    temp := 1/pl;
*/
    temp = 1 / pl;
/*
    TEMP1=CK2*TEMP
    temp1 := ck2*temp;
*/
    temp1 = CK2 * temp;
/*
    TEMP2=TEMP1*TEMP
    temp2 := temp1*temp;
*/
    temp2 = temp1 * temp;
/*
    * UPDATE FOR SHORT PERIODICS
    { Update for short periodics }
    RK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U
    rk := r*(1 - 1.5*temp2*betal*x3thm1) + 0.5*temp1*x1mth2*cos2u;


*/
    double rk = r * (1 - 1.5 * temp2 * betal * x3thm1) + 0.5 * temp1 * x1mth2 * cos2u;
/*
    UK=U-.25*TEMP2*X7THM1*SIN2U
    uk := u - 0.25*temp2*x7thm1*sin2u;


*/
    double uk = u - 0.25 * temp2 * x7thm1 * sin2u;
/*
    XNODEK=XNODE+1.5*TEMP2*COSIO*SIN2U
    xnodek := xnode + 1.5*temp2*cosio*sin2u;


*/
    double xnodek = xnode + 1.5 * temp2 * cosio * sin2u;
/*
    XINCK=XINCL+1.5*TEMP2*COSIO*SINIO*COS2U
    xinck := xincl + 1.5*temp2*cosio*sinio*cos2u;


*/
    double xinck = tle.orbitInclination + 1.5 * temp2 * cosio * sinio * cos2u;
/*
    RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U
    rdotk := rdot - xn*temp1*x1mth2*sin2u;


*/
    double rdotk = rdot - xn * temp1 * x1mth2 * sin2u;
/*
    RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5*X3THM1)
    rfdotk := rfdot + xn*temp1*(x1mth2*cos2u + 1.5*x3thm1);


*/
    double rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + 1.5 * x3thm1);
/*
    * ORIENTATION VECTORS
    { Orientation vectors }
    SINUK=SIN(UK)
    sinuk := Sin(uk);
*/
    double sinuk = Math.sin(uk);
/*
    COSUK=COS(UK)
    cosuk := Cos(uk);
*/
    double cosuk = Math.cos(uk);
/*
    SINIK=SIN(XINCK)
    sinik := Sin(xinck);
*/
    double sinik = Math.sin(xinck);
/*
    COSIK=COS(XINCK)
    cosik := Cos(xinck);
*/
    double cosik = Math.cos(xinck);
/*
    SINNOK=SIN(XNODEK)
    sinnok := Sin(xnodek);
*/
    double sinnok = Math.sin(xnodek);
/*
    COSNOK=COS(XNODEK)
    cosnok := Cos(xnodek);
*/
    double cosnok = Math.cos(xnodek);
/*
    XMX=-SINNOK*COSIK
    xmx := -sinnok*cosik;

*/
    double xmx = -sinnok * cosik;
/*
    XMY=COSNOK*COSIK
    xmy := cosnok*cosik;
*/
    double xmy = cosnok * cosik;
/*
    UX=XMX*SINUK+COSNOK*COSUK
    ux := xmx*sinuk + cosnok*cosuk;

*/
    double ux = xmx * sinuk + cosnok * cosuk;
/*
    UY=XMY*SINUK+SINNOK*COSUK
    uy := xmy*sinuk + sinnok*cosuk;
*/
    double uy = xmy * sinuk + sinnok * cosuk;
/*
    UZ=SINIK*SINUK
    uz := sinik*sinuk;
*/
    double uz = sinik * sinuk;
/*
    VX=XMX*COSUK-COSNOK*SINUK
    vx := xmx*cosuk - cosnok*sinuk;

*/
    double vx = xmx * cosuk - cosnok * sinuk;
/*
    VY=XMY*COSUK-SINNOK*SINUK
    vy := xmy*cosuk - sinnok*sinuk;
*/
    double vy = xmy * cosuk - sinnok * sinuk;
/*
    VZ=SINIK*COSUK
    vz := sinik*cosuk;
*/
    double vz = sinik * cosuk;
/*
    * POSITION AND VELOCITY
    { Position and velocity }
    X=RK*UX
    x := rk*ux;  pos[1] := x;

*/
    double x = rk * ux;
/*
    Y=RK*UY
    y := rk*uy;  pos[2] := y;
*/
    double y = rk * uy;
/*
    Z=RK*UZ
    z := rk*uz;  pos[3] := z;
*/
    double z = rk * uz;
/*
    XDOT=RDOTK*UX+RFDOTK*VX
    xdot := rdotk*ux + rfdotk*vx;  vel[1] := xdot;

*/
    double xdot = rdotk * ux + rfdotk * vx;
/*
    YDOT=RDOTK*UY+RFDOTK*VY
    ydot := rdotk*uy + rfdotk*vy;  vel[2] := ydot;
*/
    double ydot = rdotk * uy + rfdotk * vy;
/*
    ZDOT=RDOTK*UZ+RFDOTK*VZ
    zdot := rdotk*uz + rfdotk*vz;  vel[3] := zdot;
*/
    double zdot = rdotk * uz + rfdotk * vz;

    double xkm = x * EARTH_RADIUS;
    double ykm = y * EARTH_RADIUS;
    double zkm = z * EARTH_RADIUS;
    double altitude = (Math.sqrt(Math.pow(x, 2) + Math.pow(y, 2) + Math.pow(z, 2)) - 1) * EARTH_RADIUS;

    double xdotkmps = xdot * EARTH_RADIUS / 60;
    double ydotkmps = ydot * EARTH_RADIUS / 60;
    double zdotkmps = zdot * EARTH_RADIUS / 60;
    double velocity = Math.sqrt(Math.pow(xdot, 2) + Math.pow(ydot, 2) + Math.pow(zdot, 2)) * EARTH_RADIUS / 60;

  }   /* End sgp4 */

}   /* End class Computations */