Astro-satpass

 view release on metacpan or  search on metacpan

lib/Astro/Coord/ECI/TLE.pm  view on Meta::CPAN

#*    zcosi  , zsini  , zcosil , zsinil ,
#*    zx          -
#*    zy          -
#*
#*  coupling      :
#*    none.
#*
#*  references    :
#*    hoots, roehrich, norad spacetrack report #3 1980
#*    hoots, norad spacetrack report #6 1986
#*    hoots, schumacher and glover 2004
#*    vallado, crawford, hujsak, kelso  2006
#*------------------------------------------------------------------------------

sub _r_dscom {
    my ($self, $tc) = @_;
    my $parm = $self->{&TLE_INIT}{TLE_sgp4r}
        or confess "Programming error - Sgp4r not initialized";
    my $init = $parm->{init}
        or confess "Programming error - Sgp4r initialization not in progress";

#* -------------------------- Local Variables --------------------------
    my ($c1ss, $c1l, $zcosis, $zsinis, $zsings, $zcosgs, $zes, $zel);

    my ($a1, $a2, $a3, $a4, $a5, $a6, $a7, $a8, $a9, $a10, $betasq, $cc,
        $ctem, $stem, $x1, $x2, $x3, $x4, $x5, $x6, $x7, $x8, $xnodce,
        $xnoi, $zcosg, $zcosgl, $zcosh, $zcoshl, $zcosi, $zcosil,
        $zsing, $zsingl, $zsinh, $zsinhl, $zsini, $zsinil, $zx, $zy);
#>>>>trw	INCLUDE 'ASTMATH.CMN'

#* ------------------------------ Constants ----------------------------
    $zes=  0.01675;
    $zel=  0.0549;
    $c1ss=  2.9864797e-06;
    $c1l=  4.7968065e-07;
    $zsinis=  0.39785416;
    $zcosis=  0.91744867;
    $zcosgs=  0.1945905;

    $zsings= -0.98088458;
#* ----------------- DEEP SPACE PERIODICS INITIALIZATION ---------------
    $init->{xn}= $parm->{meanmotion};
    $init->{eccm}= $parm->{eccentricity};
    $init->{snodm}= sin($parm->{ascendingnode});
    $init->{cnodm}= cos($parm->{ascendingnode});
    $init->{sinomm}= sin($parm->{argumentofperigee});
    $init->{cosomm}= cos($parm->{argumentofperigee});
    $init->{sinim}= sin($parm->{inclination});
    $init->{cosim}= cos($parm->{inclination});
    $init->{emsq}= $init->{eccm}*$init->{eccm};
    $betasq= 1-$init->{emsq};

    $init->{rtemsq}= sqrt($betasq);
#* --------------------- INITIALIZE LUNAR SOLAR TERMS ------------------
    $parm->{peo}= 0;
    $parm->{pinco}= 0;
    $parm->{plo}= 0;
    $parm->{pgho}= 0;
    $parm->{pho}= 0;
    $init->{day}= $self->{ds50}+ 18261.5 + $tc/1440;
    $xnodce= fmod(4.523602 - 0.00092422029*$init->{day}, &SGP_TWOPI);
    $stem= sin($xnodce);
    $ctem= cos($xnodce);
    $zcosil= 0.91375164 - 0.03568096*$ctem;
    $zsinil= sqrt(1 - $zcosil*$zcosil);
    $zsinhl= 0.089683511*$stem/ $zsinil;
    $zcoshl= sqrt(1 - $zsinhl*$zsinhl);
    $init->{gam}= 5.8351514 + 0.001944368*$init->{day};
    $zx= 0.39785416*$stem/$zsinil;
    $zy= $zcoshl*$ctem+ 0.91744867*$zsinhl*$stem;
    $zx= atan2($zx, $zy);
    $zx= $init->{gam}+ $zx- $xnodce;
    $zcosgl= cos($zx);

    $zsingl= sin($zx);
#* ---------------------------- DO SOLAR TERMS -------------------------
    $zcosg= $zcosgs;
    $zsing= $zsings;
    $zcosi= $zcosis;
    $zsini= $zsinis;
    $zcosh= $init->{cnodm};
    $zsinh= $init->{snodm};
    $cc= $c1ss;

    $xnoi= 1 / $init->{xn};
    foreach my $lsflg (1 .. 2) {
        $a1=   $zcosg*$zcosh+ $zsing*$zcosi*$zsinh;
        $a3=  -$zsing*$zcosh+ $zcosg*$zcosi*$zsinh;
        $a7=  -$zcosg*$zsinh+ $zsing*$zcosi*$zcosh;
        $a8=   $zsing*$zsini;
        $a9=   $zsing*$zsinh+ $zcosg*$zcosi*$zcosh;
        $a10=   $zcosg*$zsini;
        $a2=   $init->{cosim}*$a7+ $init->{sinim}*$a8;
        $a4=   $init->{cosim}*$a9+ $init->{sinim}*$a10;
        $a5=  -$init->{sinim}*$a7+ $init->{cosim}*$a8;

        $a6=  -$init->{sinim}*$a9+ $init->{cosim}*$a10;
        $x1=  $a1*$init->{cosomm}+ $a2*$init->{sinomm};
        $x2=  $a3*$init->{cosomm}+ $a4*$init->{sinomm};
        $x3= -$a1*$init->{sinomm}+ $a2*$init->{cosomm};
        $x4= -$a3*$init->{sinomm}+ $a4*$init->{cosomm};
        $x5=  $a5*$init->{sinomm};
        $x6=  $a6*$init->{sinomm};
        $x7=  $a5*$init->{cosomm};

        $x8=  $a6*$init->{cosomm};
        $init->{z31}= 12*$x1*$x1- 3*$x3*$x3;
        $init->{z32}= 24*$x1*$x2- 6*$x3*$x4;
        $init->{z33}= 12*$x2*$x2- 3*$x4*$x4;
        $init->{z1}=  3* ($a1*$a1+ $a2*$a2) +
            $init->{z31}*$init->{emsq};
        $init->{z2}=  6* ($a1*$a3+ $a2*$a4) +
            $init->{z32}*$init->{emsq};
        $init->{z3}=  3* ($a3*$a3+ $a4*$a4) +
            $init->{z33}*$init->{emsq};
        $init->{z11}= -6*$a1*$a5+ $init->{emsq}*
            (-24*$x1*$x7-6*$x3*$x5);
        $init->{z12}= -6* ($a1*$a6+ $a3*$a5) + $init->{emsq}* (
            -24*($x2*$x7+$x1*$x8) - 6*($x3*$x6+$x4*$x5) );
        $init->{z13}= -6*$a3*$a6+ $init->{emsq}*(-24*$x2*$x8-
            6*$x4*$x6);

lib/Astro/Coord/ECI/TLE.pm  view on Meta::CPAN

#*  Locals        :
#*    ainv2       -
#*    aonv        -
#*    cosisq      -
#*    eoc         -
#*    f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543        -
#*    g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533        -
#*    sini2       -
#*    temp, temp1 -
#*    Theta       -
#*    xno2        -
#*
#*  Coupling      :
#*    getgravconst-
#*
#*  references    :
#*    hoots, roehrich, norad spacetrack report #3 1980
#*    hoots, norad spacetrack report #6 1986
#*    hoots, schumacher and glover 2004
#*    vallado, crawford, hujsak, kelso  2006
#*------------------------------------------------------------------------------

sub _r_dsinit {
    my ($self, $t, $tc) = @_;
    my $parm = $self->{&TLE_INIT}{TLE_sgp4r}
        or confess "Programming error - Sgp4r not initialized";
    my $init = $parm->{init}
        or confess "Programming error - Sgp4r initialization not in progress";

#* -------------------------- Local Variables --------------------------
    my ($ainv2, $aonv, $cosisq, $eoc, $f220, $f221, $f311, $f321, $f322,
        $f330, $f441, $f442, $f522, $f523, $f542, $f543, $g200, $g201,
        $g211, $g300, $g310, $g322, $g410, $g422, $g520, $g521, $g532,
        $g533, $ses, $sgs, $sghl, $sghs, $shs, $shl, $sis, $sini2, $sls,
        $temp, $temp1, $theta, $xno2);

    my ($q22, $q31, $q33, $root22, $root44, $root54, $rptim, $root32,
        $root52, $znl, $zns, $emo, $emsqo);
#>>>>trw	INCLUDE 'ASTMATH.CMN'

    $q22= 1.7891679e-06;
    $q31= 2.1460748e-06;
    $q33= 2.2123015e-07;
    $root22= 1.7891679e-06;
    $root44= 7.3636953e-09;
    $root54= 2.1765803e-09;
    $rptim= 0.0043752690880113;
    $root32= 3.7393792e-07;
    $root52= 1.1428639e-07;
#>>>>trw	X2o3   = 2.0D0 / 3.0D0
    $znl= 0.00015835218;

    $zns= 1.19459e-05;

#>>>>trw	CALL getgravconst( whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2 )
#* ------------------------ DEEP SPACE INITIALIZATION ------------------
    $parm->{irez}= 0;
    if (($init->{xn} < 0.0052359877) && ($init->{xn} > 0.0034906585)) {
        $parm->{irez}= 1;
    }
    if (($init->{xn} >= 0.00826) && ($init->{xn} <= 0.00924) &&
        ($init->{eccm} >= 0.5)) {
        $parm->{irez}= 2;

    }
#* ---------------------------- DO SOLAR TERMS -------------------------
    $ses=  $init->{ss1}*$zns*$init->{ss5};
    $sis=  $init->{ss2}*$zns*($init->{sz11}+ $init->{sz13});
    $sls= -$zns*$init->{ss3}*($init->{sz1}+ $init->{sz3}- 14 -
        6*$init->{emsq});
    $sghs=  $init->{ss4}*$zns*($init->{sz31}+ $init->{sz33}- 6);
    $shs= -$zns*$init->{ss2}*($init->{sz21}+ $init->{sz23});
#c       sgp4fix for 180 deg incl
    if (($init->{inclm} < 0.052359877) || ($init->{inclm} >
        &SGP_PI-0.052359877)) {
        $shs= 0;
    }
    if ($init->{sinim} != 0) {
        $shs= $shs/$init->{sinim};
    }

    $sgs= $sghs- $init->{cosim}*$shs;
#* ----------------------------- DO LUNAR TERMS ------------------------
    $parm->{dedt}= $ses+ $init->{s1}*$znl*$init->{s5};
    $parm->{didt}= $sis+ $init->{s2}*$znl*($init->{z11}+ $init->{z13});
    $parm->{dmdt}= $sls- $znl*$init->{s3}*($init->{z1}+ $init->{z3}- 14
        - 6*$init->{emsq});
    $sghl= $init->{s4}*$znl*($init->{z31}+ $init->{z33}- 6);
    $shl= -$znl*$init->{s2}*($init->{z21}+ $init->{z23});
#c       sgp4fix for 180 deg incl
    if (($init->{inclm} < 0.052359877) || ($init->{inclm} >
        &SGP_PI-0.052359877)) {
        $shl= 0;
    }
    $parm->{domdt}= $sgs+$sghl;
    $parm->{dnodt}= $shs;
    if ($init->{sinim} !=  0) {
        $parm->{domdt}=
            $parm->{domdt}-$init->{cosim}/$init->{sinim}*$shl;
        $parm->{dnodt}= $parm->{dnodt}+$shl/$init->{sinim};

    }
#* --------------- CALCULATE DEEP SPACE RESONANCE EFFECTS --------------
    $init->{dndt}= 0;
    $theta= fmod($parm->{gsto}+ $tc*$rptim, &SGP_TWOPI);
    $init->{eccm}= $init->{eccm}+ $parm->{dedt}*$t;
    $init->{emsq}= $init->{eccm}**2;
    $init->{inclm}= $init->{inclm}+ $parm->{didt}*$t;
    $init->{argpm}= $init->{argpm}+ $parm->{domdt}*$t;
    $init->{nodem}= $init->{nodem}+ $parm->{dnodt}*$t;
    $init->{mm}= $init->{mm}+ $parm->{dmdt}*$t;
#c   sgp4fix for negative inclinations
#c   the following if statement should be commented out
#c           IF(Inclm .lt. 0.0D0) THEN
#c             Inclm  = -Inclm
#c             Argpm  = Argpm-PI
#c             nodem = nodem+PI
#c           ENDIF

#* ------------------ Initialize the resonance terms -------------------
    if ($parm->{irez} !=  0) {

lib/Astro/Coord/ECI/TLE.pm  view on Meta::CPAN

        $su= $su- 0.25*$temp2*$parm->{x7thm1}*$sin2u;
        $xnode= $nodep+ 1.5*$temp2*$cosip*$sin2u;
        $xinc= $xincp+ 1.5*$temp2*$cosip*$sinip*$cos2u;
        $mv= $rdotl- $xn*$temp1*$parm->{x1mth2}*$sin2u/ $parm->{xke};

        $rvdot= $rvdotl+ $xn*$temp1*
            ($parm->{x1mth2}*$cos2u+1.5*$parm->{con41}) / $parm->{xke};
#* ------------------------- ORIENTATION VECTORS -----------------------
        $sinsu=  sin($su);
        $cossu=  cos($su);
        $snod=  sin($xnode);
        $cnod=  cos($xnode);
        $sini=  sin($xinc);
        $cosi=  cos($xinc);
        $xmx= -$snod*$cosi;
        $xmy=  $cnod*$cosi;
        $ux=  $xmx*$sinsu+ $cnod*$cossu;
        $uy=  $xmy*$sinsu+ $snod*$cossu;
        $uz=  $sini*$sinsu;
        $vx=  $xmx*$cossu- $cnod*$sinsu;
        $vy=  $xmy*$cossu- $snod*$sinsu;

        $vz=  $sini*$cossu;
#* ----------------------- POSITION AND VELOCITY -----------------------
        $r[1] = $mr*$ux* $parm->{radiusearthkm};
        $r[2] = $mr*$uy* $parm->{radiusearthkm};
        $r[3] = $mr*$uz* $parm->{radiusearthkm};
        $v[1] = ($mv*$ux+ $rvdot*$vx) * $vkmpersec;
        $v[2] = ($mv*$uy+ $rvdot*$vy) * $vkmpersec;
        $v[3] = ($mv*$uz+ $rvdot*$vz) * $vkmpersec;

    }
#* --------------------------- ERROR PROCESSING ------------------------
#c     sgp4fix for decaying satellites
    if ($mr <  1) {
#c          write(*,*) '# decay condition ',mr
        $self->{model_error}= &SGP4R_ERROR_6;

    }
#c        INCLUDE 'debug7.for'

#>>>>trw	RETURN

    $self->__universal( $time );
    $self->eci (@r[1..3], @v[1..3]);
    $self->equinox_dynamical ($self->{epoch_dynamical});
    return $self;
}

=begin comment

The following code was converted from the Fortran reference
implementation, but is not used by this code.

#* -----------------------------------------------------------------------------
#*
#*                           FUNCTION GSTIME
#*
#*  This function finds the Greenwich SIDEREAL time.  Notice just the INTEGER
#*    part of the Julian Date is used for the Julian centuries calculation.
#*    We use radper Solar day because we're multiplying by 0-24 solar hours.
#*
#*  Author        : David Vallado                  719-573-2600    1 Mar 2001
#*
#*  Inputs          Description                    Range / Units
#*    JD          - Julian Date                    days from 4713 BC
#*
#*  OutPuts       :
#*    GSTIME      - Greenwich SIDEREAL Time        0 to 2Pi rad
#*
#*  Locals        :
#*    Temp        - Temporary variable for reals   rad
#*    TUT1        - Julian Centuries from the
#*                  Jan 1, 2000 12 h epoch (UT1)
#*
#*  Coupling      :
#*
#*  References    :
#*    Vallado       2007, 194, Eq 3-45
#* -----------------------------------------------------------------------------

sub _r_gstime {
    my ($jd) = @_;
    my $gstime;
#* ----------------------------  Locals  -------------------------------

    my ($temp, $tut1);
#>>>>trw	INCLUDE 'astmath.cmn'

    $tut1= ( $$jd- 2451545 ) / 36525;
    $temp= - 6.2e-06*$tut1*$tut1*$tut1+ 0.093104*$tut1*$tut1+
        (876600*3600 + 8640184.812866)*$tut1+ 67310.54841;

    $temp= fmod($temp*&SGP_DE2RA/240, &SGP_TWOPI);
    if ( $temp <  0 ) {
        $temp= $temp+ &SGP_TWOPI;

    }

    $gstime= $temp;
    return $gstime;
}

=end comment

=cut

#* -----------------------------------------------------------------------------
#*
#*                           function getgravconst
#*
#*  this function gets constants for the propagator. note that mu is identified to
#*    facilitiate comparisons with newer models.
#*
#*  author        : david vallado                  719-573-2600   21 jul 2006
#*
#*  inputs        :
#*    whichconst  - which set of constants to use  721, 72, 84
#*
#*  outputs       :
#*    tumin       - minutes in one time unit



( run in 0.738 second using v1.01-cache-2.11-cpan-a49fcb8fa48 )