horiz_coo

--------------------------------------------------------------------
 horiz_coo function        Horizontal coordinates conversion
                       Converting from equatorial coordinates
                       to Horizontal coordinates and visa
                       versa.
 input  : - Two columns matrix of coordinates,
            (Long & Lat) | (Az & Alt) in radians
          - vector of JDs + UT fraction of day,
            if scalar value is given, then it duplicated
            for all coordinates.
          - Geodetic Coordinates, east long & north lat in radians
            if scalar value is given then it is duplicate for
            all coordinates.
          - Direction,
            'h' - from equatorial to horizontal (default).
            'e' - from horizontal to equatorial.
 output : - two column matrix of output coordinates.
 Tested : matlab 5.3
     By : Eran O. Ofek           August 1999
    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
--------------------------------------------------------------------

0001 function OutCoo=horiz_coo(InCoo,JD,TopoPos,Direction);
0002 %--------------------------------------------------------------------
0003 % horiz_coo function        Horizontal coordinates conversion
0004 %                       Converting from equatorial coordinates
0005 %                       to Horizontal coordinates and visa
0006 %                       versa.
0007 % input  : - Two columns matrix of coordinates,
0008 %            (Long & Lat) | (Az & Alt) in radians
0009 %          - vector of JDs + UT fraction of day,
0010 %            if scalar value is given, then it duplicated
0011 %            for all coordinates.
0012 %          - Geodetic Coordinates, east long & north lat in radians
0013 %            if scalar value is given then it is duplicate for
0014 %            all coordinates.
0015 %          - Direction,
0016 %            'h' - from equatorial to horizontal (default).
0017 %            'e' - from horizontal to equatorial.
0018 % output : - two column matrix of output coordinates.
0019 % Tested : matlab 5.3
0020 %     By : Eran O. Ofek           August 1999
0021 %    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
0022 %--------------------------------------------------------------------
0023 if (nargin==3),
0024    Direction = 'h';
0025 elseif (nargin==4),
0026    % no default
0027 else
0028    error('Illigal number of input arguments');
0029 end
0030 
0031 N = length(InCoo(:,1));
0032 if (length(JD)==1),
0033    JD = JD.*ones(N,1);
0034 end
0035 
0036 if (length(TopoPos(:,1))==1),
0037    TopoPos = ones(N,1)*TopoPos;
0038 end
0039 
0040 
0041 % Don't convert Geodetic latitude to Geocentric.
0042 %GeodLat = TopoPos(:,2);
0043 %GeocLatTemp = geod2geoc([TopoPos(:,1),GeodLat,zeros(N,1)],'WGS84');
0044 %GeocLat     = GeocLatTemp(:,2);
0045 %TopoPos(:,2) = GeocLat;
0046 
0047 
0048 % calculating Local Mean Sidereal Time
0049 LST=lst(JD,TopoPos(:,1),'m');
0050 
0051 
0052 if (Direction=='h'),
0053    % convert equatorial to horizontal
0054 
0055    % calculate the Hour Angle
0056    HA  = LST.*2.*pi - InCoo(:,1);
0057    Dec = InCoo(:,2);
0058    Lat = TopoPos(:,2);
0059 
0060    SinAlt = sin(Dec).*sin(Lat) + cos(Dec).*cos(HA).*cos(Lat);
0061    CosAlt = sqrt(1-SinAlt.*SinAlt);
0062 
0063    SinAz  = (-cos(Dec).*sin(HA))./CosAlt;
0064    CosAz  = (sin(Dec).*cos(Lat) - cos(Dec).*cos(HA).*sin(Lat))./CosAlt;
0065 
0066    Az     = atan2(SinAz, CosAz);
0067    Alt    = asin(SinAlt);
0068 
0069    I      = find(Az<0);
0070    Az(I)  = 2.*pi+Az(I);
0071 
0072    OutCoo = [Az, Alt];
0073 elseif (Direction=='e'),
0074    Az     = InCoo(:,1);
0075    Alt    = InCoo(:,2);
0076    Lat = TopoPos(:,2);
0077 
0078    SinDec = sin(Alt).*sin(Lat) + cos(Alt).*cos(Az).*cos(Lat);
0079    CosDec = sqrt(1 - SinDec.*SinDec);
0080 
0081    SinHA  = (-cos(Alt).*sin(Az))./CosDec;
0082    CosHA  = (sin(Alt).*cos(Lat) - cos(Alt).*cos(Az).*sin(Lat))./CosDec;
0083    HA     = atan2(SinHA, CosHA);
0084 
0085    RA     = LST.*2.*pi - HA;
0086    Dec    = asin(SinDec);
0087 
0088    % converting to range [0,1)
0089    RA     = 2.*pi.*(RA./(2.*pi) - floor(RA./(2.*pi)));
0090 
0091    OutCoo = [RA, Dec];
0092 else
0093    error('Illigal type of conversion');
0094 end
0095 
0096 
0097 
0098 
0099 
0100 
0101 
0102 
0103 
0104