lst

--------------------------------------------------------------------
 lst function         Local Sidereal Time, (mean or apparent),
                    for vector of JD's and a given East Longitude.
 input  : - Vector of JD, in UT1 time scale.
          - East Longitude in radians.
          - Sidereal Time Type,
            'm' - Mean (default).
            'a' - apparent.
 output : - vector of LST in fraction of day.
    By  Eran O. Ofek           August 1999
--------------------------------------------------------------------

0001 function LST=lst(JD,EastLong,STType);
0002 %--------------------------------------------------------------------
0003 % lst function         Local Sidereal Time, (mean or apparent),
0004 %                    for vector of JD's and a given East Longitude.
0005 % input  : - Vector of JD, in UT1 time scale.
0006 %          - East Longitude in radians.
0007 %          - Sidereal Time Type,
0008 %            'm' - Mean (default).
0009 %            'a' - apparent.
0010 % output : - vector of LST in fraction of day.
0011 %    By  Eran O. Ofek           August 1999
0012 %--------------------------------------------------------------------
0013 RAD = 57.29577951308232;
0014 
0015 if (nargin==2),
0016    STType = 'm';
0017 elseif (nargin==3),
0018    % do nothing
0019 else
0020    error('Illigal number of input arguments');
0021 end
0022 
0023 
0024 % convert JD to integer day + fraction of day
0025 TJD = floor(JD - 0.5) + 0.5;
0026 DayFrac = JD - TJD;
0027 
0028 T = (TJD - 2451545.0)./36525.0;
0029 
0030 GMST0UT = 24110.54841 + 8640184.812866.*T + 0.093104.*T.*T - 6.2e-6.*T.*T.*T;
0031 
0032 % convert to fraction of day in range [0 1)
0033 GMST0UT = GMST0UT./86400.0;
0034 
0035 GMST0UT = GMST0UT - floor(GMST0UT);
0036 LST = GMST0UT + 1.0027379093.*DayFrac + EastLong./(2.*pi);
0037 LST = LST - floor(LST);
0038 
0039 
0040 switch STType
0041  case {'m'}
0042     % do nothing
0043  case {'a'}
0044     % calculate nutation
0045     NutMat = nutation(JD);
0046     Obl    = obliquity(JD);
0047     EquationOfEquinox = (RAD.*3600).*NutMat(:,1).*cos(Obl)./15;
0048     LST = LST + EquationOfEquinox./86400;    
0049  otherwise
0050     error('Unknown sidereal time type');
0051 end