flagRFI_pos

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function  [flagShort flagLong sunFlag moonFlag] = flagRFI_pos(d, Antenna)

   March 17, 2011 (MAS) - Flags data based upon positional proximity to
   known sources of RFI.

   March 31, 2011 (MAS) - Also, the Sun and Moon because they mess with my
   RFI tests.

   April 15, 2011 (MAS) - Committed.
   August 8, 2011 (SJCM) - Modified for memory efficiency.  will now run on
    at least 7 hours of data where switchData is not present.

   August 11, 2011 (MAS) - Split the northern survey into two epochs:  
                           pre- and post-filters.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [horzFlag satFlag sunFlag moonFlag] = flagRFI_pos(d, Antenna)
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %  function  [flagShort flagLong sunFlag moonFlag] = flagRFI_pos(d, Antenna)
0004 %
0005 %   March 17, 2011 (MAS) - Flags data based upon positional proximity to
0006 %   known sources of RFI.
0007 %
0008 %   March 31, 2011 (MAS) - Also, the Sun and Moon because they mess with my
0009 %   RFI tests.
0010 %
0011 %   April 15, 2011 (MAS) - Committed.
0012 %   August 8, 2011 (SJCM) - Modified for memory efficiency.  will now run on
0013 %    at least 7 hours of data where switchData is not present.
0014 %
0015 %   August 11, 2011 (MAS) - Split the northern survey into two epochs:
0016 %                           pre- and post-filters.
0017 %
0018 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0019 
0020 
0021 %%%%%%%%%%%%%%
0022 % Parameters:
0023 %%%%%%%%%%%%%%
0024 
0025 epochN_1 = '15-Jul-2011:00:00:00'; % End of first Northern epoch (pre-filters).
0026 
0027 % Horizon positions:
0028 horzAzN_1 = [168.5 318];
0029 horzAzN_2 = [];
0030 horzAzS = [];
0031 
0032 % Horizon radii:
0033 horzRadN_1 = [40 40];
0034 horzRadN_2 = [];
0035 horzRadS = [];
0036 
0037 % Satellite positions:
0038 satAzN_1 = [110.5 119 138.5 146 157.5 170 176.5 183 186.5 190 199.5 202.5 205.5 208];
0039 satAzN_2 = [110.5];
0040 satAzS = [];
0041 
0042 % Satellite radii:
0043 satRadN_1 = 6;
0044 satRadN_2 = 15;
0045 satRadS = 6;
0046 
0047 % Sun and moon radii:
0048 moonRad = 7;
0049 sunRad = 9;
0050 
0051 
0052 % Observatory latitudes:
0053 obsLatN = 37.2333;
0054 obsLatS = -30.83;
0055 
0056 
0057 %%%%%%%%%%%%%%%%%%%%%%%%%%
0058 % Begin the calculations!
0059 %%%%%%%%%%%%%%%%%%%%%%%%%%
0060 
0061 
0062 % Get azimuth and elevation:
0063 az = d.antenna0.servo.apparent(:,1);
0064 el = d.antenna0.servo.apparent(:,2);
0065 
0066 % Get its length.
0067 dLength = length(az);
0068 
0069 
0070 % North
0071 if Antenna == 1
0072     
0073     % MAS:  Added second epoch of northern observing (ie, after filters
0074     % installed).
0075     if (d.antenna0.servo.utc < tstr2mjd(epochN_1))
0076 %        disp('flagRFI_pos:: Pre-');
0077         horzAz = horzAzN_1;
0078         horzRad = horzRadN_1;
0079         satAz = satAzN_1;
0080         satRad = satRadN_1;
0081         obsLat = obsLatN;
0082     else
0083 %        disp('flagRFI_pos:: Post-');
0084         horzAz = horzAzN_2;
0085         horzRad = horzRadN_2;
0086         satAz = satAzN_2;
0087         satRad = satRadN_2;
0088         obsLat = obsLatN;
0089     end
0090     
0091     % South
0092 elseif Antenna == 2
0093     horzAz = horzAzS;
0094     horzRad = horzRadS;
0095     satAz = satAzS;
0096     satRad = satRadS;
0097     obsLat = obsLatS;
0098 
0099 % Invalid choice
0100 else
0101     return;
0102 end
0103 
0104 % How many horizon positions?
0105 nHorz = length(horzAz);
0106 % How many satellites?
0107 nSat = length(satAz);
0108 
0109 
0110 
0111 % Calculate the geosynchronous elevations.
0112 satPhi = asind(sqrt(sind(obsLat)^2./(1-sind(satAz).^2*cosd(obsLat)^2)));
0113 satEl = acosd(sind(satPhi) ./ sqrt(1 + 0.178^2 - 2 * 0.178 * cosd(satPhi)));
0114 
0115 
0116 % Correct for refraction the satellite elevations.
0117 Avals = d.antenna0.tracker.refraction(:,1)*pi/180;
0118 Bvals = d.antenna0.tracker.refraction(:,2)*pi/180;
0119 A = interp1(d.array.frame.utc, Avals, d.antenna0.receiver.utc);
0120 B = interp1(d.array.frame.utc, Bvals, d.antenna0.receiver.utc);
0121 
0122 clear Avals;
0123 clear Bvals;
0124 
0125 
0126 % SJCM:  the following piece of code runs out of memory.  rewritting as a loop.
0127 %satElRep = repmat(satEl,dLength,1);
0128 %ARep = repmat(A,1,nSat);
0129 %BRep = repmat(B,1,nSat);
0130 %refracCorr = (ARep.*cosd(satElRep).*(sind(satElRep)).^3 + ...
0131 %    BRep.*sind(satElRep).*(cosd(satElRep)).^3 ) ./ ...
0132 %    ((sind(satElRep)).^4 + ARep.*(sind(satElRep)).^2 + ...
0133 %    3*BRep.*(cosd(satElRep)).^2) * 180 / pi;
0134 %satElCorr = satElRep + refracCorr;
0135 % SJCM: REWRITING IT AS A LOOP
0136 % MAS: Fixed bug in loop code.  Also added check to make sure we have at
0137 % least one satellite.
0138 if nSat > 0
0139     for m=1:nSat
0140         thisSatEl = satEl(m);
0141         thisRefracCorr = (A * cosd(thisSatEl) * sind(thisSatEl)^3 + ...
0142             B * sind(thisSatEl) * cosd(thisSatEl)^3 ) ./ ...
0143             (sind(thisSatEl)^4 + A * sind(thisSatEl)^2 + ...
0144             3*B * cosd(thisSatEl)^2) * 180 / pi;
0145 %         disp(thisSatEl);
0146 %         figure;
0147 %         plot(thisRefracCorr);
0148 %         input('sdfds');
0149 %         disp(nanmean(thisRefracCorr));
0150         satElCorr(:,m) = thisSatEl + thisRefracCorr;
0151     end
0152 
0153     clear thisSatEl;
0154     clear thisRefracCorr;
0155 end
0156 
0157 
0158 
0159 % Calculate the positions of the Sun and Moon.  Should take into account
0160 % parallax, but have not, yet.
0161 [junk, moonAz, moonEl] = calcSourceDist(d, 'moon');
0162 [junk, sunAz, sunEl] = calcSourceDist(d, 'sun');
0163 clear junk;
0164 
0165 moonRefracCorr = (A.*cosd(moonEl).*(sind(moonEl)).^3 + ...
0166     B.*sind(moonEl).*(cosd(moonEl)).^3 ) ./ ...
0167     ((sind(moonEl)).^4 + A.*(sind(moonEl)).^2 + ...
0168     3*B.*(cosd(moonEl)).^2) * 180 / pi;
0169 moonElCorr = moonEl + moonRefracCorr;
0170 
0171 % clear up some memory
0172 clear moonEl;
0173 clear moonRefracCorr;
0174 
0175 sunRefracCorr = (A.*cosd(sunEl).*(sind(sunEl)).^3 + ...
0176     B.*sind(sunEl).*(cosd(sunEl)).^3 ) ./ ...
0177     ((sind(sunEl)).^4 + A.*(sind(sunEl)).^2 + ...
0178     3*B.*(cosd(sunEl)).^2) * 180 / pi;
0179 
0180 sunElCorr = sunEl + sunRefracCorr;
0181 
0182 % clear up some memory
0183 clear sunEl;
0184 clear sunRefracCorr;
0185 
0186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0187 % OK, now do the calc:
0188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0189 
0190 % SJCM:  re-ordering flagging so we can conserve some memory.
0191 
0192 % Now the Sun and Moon!
0193 moonFlag = (spaceangle(moonAz, moonElCorr, az, el, 'deg') < moonRad);
0194 sunFlag = (spaceangle(sunAz, sunElCorr, az, el, 'deg') < sunRad);
0195 
0196 clear sunElCorr; 
0197 clear sunAz;
0198 clear moonElCorr;
0199 clear moonAz;
0200 clear A;
0201 clear B;
0202 
0203 % First the horizon flagging.  This is not quite correct...
0204 % MAS: Added check to make sure we have at least one horizon source.
0205 if nHorz > 0
0206     horzFlag = max(acosd( ...
0207         cosd(repmat(el,1,nHorz)) .* ...
0208         cosd(repmat(horzAz,dLength,1) - repmat(az,1,nHorz))) < ...
0209         repmat(horzRad,dLength,1),[],2);
0210 else
0211     horzFlag = logical(zeros(dLength,1));
0212 end
0213 
0214 % SJCM:  clear up some memory before the next step.
0215 clear nHorz
0216 
0217 
0218 % Next the satellite flagging.
0219 
0220 % SJCM: again, do a loop instead of the following code.
0221 % satFlag = max(acosd( ...
0222 %     sind(repmat(el,1,nSat)) .* sind(satElCorr) + ...
0223 %     cosd(repmat(el,1,nSat)) .* cosd(satElCorr) .* ...
0224 %     cosd(repmat(satAz,dLength,1) - repmat(az,1,nSat))) ...
0225 %     < satRad,[],2);
0226 % MAS: Fixed bug in looped code.  Added check to make sure we have at least
0227 % one satellite.
0228 satFlag = logical(zeros(size(el)));
0229 if nSat > 0
0230     for m=1:nSat
0231         thisSatParam = acosd( sind(el) .* sind(satElCorr(:,m)) + ...
0232             cosd(el) .* cosd(satElCorr(:,m)).*cosd(satAz(m) - az)) < satRad;
0233         
0234         satFlag = satFlag | thisSatParam;
0235     end
0236 end
0237 
0238 clear satElCorr;
0239 clear el;
0240 clear az;
0241 
0242 end