Add Vallado's site.m to the toolbox

[mihalybaci]

One function that I think could be added is David Vallado's site.m routine. The description from the code is:

%  this function finds the position and velocity vectors for a site.  the
%    answer is returned in the geocentric equatorial (ecef) coordinate system.
%    note that the velocity is zero because the coordinate system is fixed to
%    the earth.

I have typed up a Julia version, though I have not been able to test/verify it as I do not own a Matlab license. Here is my initial attempt:

function site(latgd, lon, alt)
    # Constants, names in parantheses are the original Matlab variable names
    Rₑ = 6378.137         # km, (re)
    FLAT = 1.0/298.257223563  # (flat)
    eₑ = sqrt(2.0*FLAT - FLAT^2)  # (eccearth)
    eₑ₂ = eₑ^2  # (eccearthsqrd)

    # -------------------------  implementation   -----------------
    sinlat = sin(latgd)

    # ------  find rdel and rk components of site vector  ---------
    cearth = Rₑ/sqrt(1.0 - (eₑ₂*sinlat^2))
    rdel = (cearth + alt)*cos(latgd)
    rk = ((1.0 - eₑ₂)*cearth + alt)*sinlat

    # ---------------  find site position vector  -----------------
    rs = [rdel * cos(lon),
            rdel * sin(lon),
            rk]

    # ---------------  find site velocity vector  -----------------
    vs = [0.0, 0.0, 0.0];

    return rs, vs
end

Here is the original Matlab code:

function [rs,vs] = site ( latgd,lon,alt );

        constastro;

        % -------------------------  implementation   -----------------
        sinlat      = sin( latgd );

        % ------  find rdel and rk components of site vector  ---------
        cearth= re / sqrt( 1.0 - ( eccearthsqrd*sinlat*sinlat ) );
        rdel  = ( cearth + alt )*cos( latgd );
        rk    = ( (1.0-eccearthsqrd)*cearth + alt )*sinlat;

        % ---------------  find site position vector  -----------------
        rs(1) = rdel * cos( lon );
        rs(2) = rdel * sin( lon );
        rs(3) = rk;
        rs = rs';

        % ---------------  find site velocity vector  -----------------
        [vs] = [0.0; 0.0; 0.0];

Thanks!

[ronisbr]

Hi @mihalybaci !

I was thinking, since everything is represented in the ECEF, the site velocity is always zero. In this case, isn't this function exactly equal to GeodetictoECEF?

[mihalybaci]

ahhhh, you're right.They give identical results (except site.m is in km not meters). I still have a lot of learning to do (about the code and the subject itself). This issue can be closed.