calcSourceDist

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function [dist,az,el] = calcSourceDist(d, source, antenna)


  calculates the distance between where the antenna is pointing and the
  source given by 'source'.  source can either be 'sun' or 'moon', and d
  must be the data structure with az/el

  sjcm

  MAS 23-Feb-2012: Added optional antenna parameter. '1' for North, 
 '2' for South.

  Paddy 10-Apr-2013: Re-wrote & simplified.
                     New ephem files with angles in degrees
                     no need any more to convert HMS DMS to degrees,
                     hence eliminating a couple of bugs.
                     Antenna '3' = Hartebeesthoek
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [dist,az,el] = calcSourceDist(d, source, antenna)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function [dist,az,el] = calcSourceDist(d, source, antenna)
0006 %
0007 %
0008 %  calculates the distance between where the antenna is pointing and the
0009 %  source given by 'source'.  source can either be 'sun' or 'moon', and d
0010 %  must be the data structure with az/el
0011 %
0012 %  sjcm
0013 %
0014 %  MAS 23-Feb-2012: Added optional antenna parameter. '1' for North,
0015 % '2' for South.
0016 %
0017 %  Paddy 10-Apr-2013: Re-wrote & simplified.
0018 %                     New ephem files with angles in degrees
0019 %                     no need any more to convert HMS DMS to degrees,
0020 %                     hence eliminating a couple of bugs.
0021 %                     Antenna '3' = Hartebeesthoek
0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0023 
0024 if ~exist('antenna','var')
0025     antenna = 1;
0026 end
0027 
0028 switch source
0029   case 'moon'
0030     % load the source ephemeris in
0031 %    moon = load('constants/moon_ephem2.m');
0032     if antenna == 1
0033         ephem = load('constants/moon_n3.ephem');
0034     else
0035         ephem = load('constants/moon_s.ephem');
0036     end
0037   case 'sun'
0038     if antenna == 1
0039         ephem = load('constants/sun_n3.ephem');
0040     else
0041         ephem = load('constants/sun_s.ephem');
0042     end
0043  otherwise
0044     error('source not recognized');
0045 end
0046 
0047 mjd = ephem(:,1) - 2400000.5;
0048 ra  = ephem(:,2)/15 ;  % in hours
0049 dec = ephem(:,3);  % in degrees
0050  
0051 
0052 % spline won't work unless RA is unwrapped.
0053 utc = d.antenna0.receiver.utc;
0054 start_utc = find(mjd < utc(1),1,'last');
0055 end_utc = find(mjd > utc(end),1,'first');
0056 obsrange = [start_utc end_utc];
0057 dra = ra(obsrange) - ra(obsrange-1);
0058 if max(dra) > 12.
0059     ra = wrap180(ra*15)/15;
0060 end 
0061 raInt  = interp1(mjd, ra, utc, 'linear');
0062 decInt = interp1(mjd, dec, utc, 'linear');
0063 
0064 % make sure we don't have repeated data in tracker.utc and tracker.lst
0065 if(length(d.antenna0.tracker.utc)~=length(unique(d.antenna0.tracker.utc)))
0066   indBadTime = zeros(size(d.antenna0.tracker.utc));
0067   aa = sort(d.antenna0.tracker.utc);
0068   f = find(diff(aa)==0);
0069   for m=1:length(f)
0070     ff = find(d.antenna0.tracker.utc == aa(f(m)));
0071     indBadTime(ff(2:length(ff))) = 1;
0072   end
0073 
0074   trackutc = d.antenna0.tracker.utc(~indBadTime);
0075   tracklst = d.antenna0.tracker.lst(~indBadTime);
0076 else  
0077   trackutc = d.antenna0.tracker.utc;
0078   tracklst = d.antenna0.tracker.lst;
0079 end
0080 
0081 lstInt = interp1(trackutc, tracklst, utc, 'linear', 'extrap');
0082 clear utc;
0083 
0084 % from RA/DEC to az/el
0085 lat = d.antenna0.tracker.siteActual(1,2);
0086 
0087 [az,el] = radec2azel(raInt, decInt, lstInt, lat);
0088 
0089 % next we find the distance to what we're currently pointing at
0090 dist = spaceangle(az, el, d.antenna0.servo.apparent(:,1), ...
0091                   d.antenna0.servo.apparent(:,2), 'deg');
0092 
0093 f = isnan(dist);
0094 dist(f) = 0;
0095 return;
0096