time   calls  line
                  1 function OutCoo=horiz_coo(InCoo,JD,TopoPos,Direction);
                  2 %--------------------------------------------------------------------
                  3 % horiz_coo function        Horizontal coordinates conversion
                  4 %                       Converting from equatorial coordinates
                  5 %                       to Horizontal coordinates and visa
                  6 %                       versa.
                  7 % input  : - Two columns matrix of coordinates,
                  8 %            (Long & Lat) | (Az & Alt) in radians
                  9 %          - vector of JDs + UT fraction of day,
                 10 %            if scalar value is given, then it duplicated
                 11 %            for all coordinates.
                 12 %          - Geodetic Coordinates, east long & north lat in radians
                 13 %            if scalar value is given then it is duplicate for
                 14 %            all coordinates.
                 15 %          - Direction,
                 16 %            'h' - from equatorial to horizontal (default).
                 17 %            'e' - from horizontal to equatorial.
                 18 % output : - two column matrix of output coordinates.
                 19 % Tested : matlab 5.3
                 20 %     By : Eran O. Ofek           August 1999
                 21 %    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
                 22 %--------------------------------------------------------------------
             1   23 if (nargin==3), 
                 24    Direction = 'h';
             1   25 elseif (nargin==4), 
                 26    % no default
                 27 else
                 28    error('Illigal number of input arguments');
                 29 end
                 30 
  0.01       1   31 N = length(InCoo(:,1)); 
             1   32 if (length(JD)==1), 
                 33    JD = JD.*ones(N,1);
                 34 end
                 35 
             1   36 if (length(TopoPos(:,1))==1), 
  0.01       1   37    TopoPos = ones(N,1)*TopoPos; 
             1   38 end 
                 39 
                 40 
                 41 % Don't convert Geodetic latitude to Geocentric.
                 42 %GeodLat = TopoPos(:,2);
                 43 %GeocLatTemp = geod2geoc([TopoPos(:,1),GeodLat,zeros(N,1)],'WGS84');
                 44 %GeocLat     = GeocLatTemp(:,2);
                 45 %TopoPos(:,2) = GeocLat;
                 46 
                 47 
                 48 % calculating Local Mean Sidereal Time
  0.13       1   49 LST=lst(JD,TopoPos(:,1),'m'); 
                 50 
                 51 
             1   52 if (Direction=='h'), 
                 53    % convert equatorial to horizontal
                 54 
                 55    % calculate the Hour Angle
                 56    HA  = LST.*2.*pi - InCoo(:,1);
                 57    Dec = InCoo(:,2);
                 58    Lat = TopoPos(:,2);
                 59 
                 60    SinAlt = sin(Dec).*sin(Lat) + cos(Dec).*cos(HA).*cos(Lat);
                 61    CosAlt = sqrt(1-SinAlt.*SinAlt);
                 62 
                 63    SinAz  = (-cos(Dec).*sin(HA))./CosAlt;
                 64    CosAz  = (sin(Dec).*cos(Lat) - cos(Dec).*cos(HA).*sin(Lat))./CosAlt;
                 65 
                 66    Az     = atan2(SinAz, CosAz);
                 67    Alt    = asin(SinAlt);
                 68 
                 69    I      = find(Az<0);
                 70    Az(I)  = 2.*pi+Az(I);
                 71 
                 72    OutCoo = [Az, Alt];
             1   73 elseif (Direction=='e'), 
  0.01       1   74    Az     = InCoo(:,1); 
             1   75    Alt    = InCoo(:,2); 
  0.01       1   76    Lat = TopoPos(:,2); 
                 77 
  0.09       1   78    SinDec = sin(Alt).*sin(Lat) + cos(Alt).*cos(Az).*cos(Lat); 
  0.01       1   79    CosDec = sqrt(1 - SinDec.*SinDec); 
                 80 
  0.05       1   81    SinHA  = (-cos(Alt).*sin(Az))./CosDec; 
  0.09       1   82    CosHA  = (sin(Alt).*cos(Lat) - cos(Alt).*cos(Az).*sin(Lat))./CosDec; 
  0.04       1   83    HA     = atan2(SinHA, CosHA); 
                 84 
  0.02       1   85    RA     = LST.*2.*pi - HA; 
  0.04       1   86    Dec    = asin(SinDec); 
                 87 
                 88    % converting to range [0,1)
  0.02       1   89    RA     = 2.*pi.*(RA./(2.*pi) - floor(RA./(2.*pi))); 
                 90 
  0.01       1   91    OutCoo = [RA, Dec]; 
                 92 else
                 93    error('Illigal type of conversion');
                 94 end

Other subfunctions in this file are not included in this listing.