calculateAlpha

 This function calculates the alpha values. It calls either the classic
 version (_DD), the PolOnly version or the filtered version, depending on 
 the selection flag.

 OGK, 12 March 2012

0001 function [t,A,G,varA,varG,T,r,horiz,equa,offStartPosRet,onEndPosRet,offEndPosRet,onStartPosRet] = calculateAlpha(d,selection)
0002 % This function calculates the alpha values. It calls either the classic
0003 % version (_DD), the PolOnly version or the filtered version, depending on
0004 % the selection flag.
0005 %
0006 % OGK, 12 March 2012
0007 %
0008 
0009 debug = 0;
0010 
0011 if(~isfield(d.antenna0.servo, 'equa')) % added in here for when you need to calculate alphas using data in pipeline 01/09/11
0012   display('CalculateAlpha:: Calculating RA/DEC');
0013   long=-118.2822;
0014   lat=37.2339;
0015   
0016   az = d.antenna0.servo.apparent(:,1);
0017   el = d.antenna0.servo.apparent(:,2);
0018   jd=mjd2jd(d.antenna0.receiver.utc);
0019   [equa] = horiz_coo([pi/180*(az) pi/180*(el)],jd,[pi/180*(long) ...
0020     pi/180*(lat)],'e');
0021   d.antenna0.servo.equa=equa;
0022   clear equa;
0023   clear az;
0024   clear el;
0025 end
0026 
0027 % Find the noise diode events:
0028 [off1StartPos, onEndPos, off2EndPos, onStartPos, off1EndPos, off2StartPos, N] = calcNoiseIndices(d);
0029 
0030 % If no noise diode events:
0031 if N == 0
0032     t = [];
0033     A = [];
0034     G = [];
0035     varA = [];
0036     varG = [];
0037     T = [];
0038     r = [];
0039     horiz=[];
0040     equa=[];
0041     offStartPosRet=[];
0042     onEndPosRet=[];
0043     offEndPosRet=[];
0044     onStartPosRet=[];
0045     return
0046 end
0047 
0048 onStartPosRet = onStartPos;
0049 onEndPosRet = onEndPos;
0050 if d.antenna0.receiver.utc(onStartPos(1)) < tstr2mjd('08-SEP-2010:19:00:00')
0051     offStartPosRet = off1StartPos;
0052     offEndPosRet = off1EndPos;
0053 else
0054     offStartPosRet = [off1StartPos off2StartPos];
0055     offEndPosRet = [off1EndPos off2EndPos];
0056 end
0057 
0058 % now calculate the Ion and Ioff vectors for each noise diode event.
0059 Ion = {};
0060 Ioff = {};
0061 t = zeros(N,1);
0062 horiz = zeros(N,2);
0063 equa = zeros(N,2);
0064 for k=1:N
0065     Ons = onStartPos(k);
0066     One = onEndPos(k);
0067 
0068     Ion{k} = (Ons:One)';
0069 
0070     if d.antenna0.receiver.utc(Ons) < tstr2mjd('08-SEP-2010:19:00:00')
0071         Offs = off1StartPos(k);
0072         Offe = off1EndPos(k);
0073 
0074         Ioff{k} = (Offs:Offe)';
0075     else
0076         Offs1 = off1StartPos(k);
0077         Offe1 = off1EndPos(k);
0078         Offs2 = off2StartPos(k);
0079         Offe2 = off2EndPos(k);
0080 
0081         Ioff{k} = cat(1,(Offs1:Offe1)',(Offs2:Offe2)');
0082     end
0083     
0084     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0085     % The MJD time to assign to each point:
0086     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0087     t(k) = nanmean(d.antenna0.receiver.utc(Ion{k}));
0088 
0089     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0090     % The coordinates to assign to each point:
0091     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0092     horiz(k,1:2) = nanmean(d.antenna0.servo.apparent(Ion{k},1:2),1);
0093     horiz(k,1) = nanmean(d.antenna0.servo.az(Ion{k}),1);
0094     horiz(k,2) = nanmean(d.antenna0.servo.el(Ion{k}),1);
0095     equa(k,1:2) = nanmean(d.antenna0.servo.equa(Ion{k},1:2),1);
0096 end
0097 
0098 %  Apply the flags we already have
0099 if(isfield(d, 'flags'))
0100   if(isfield(d.flags, 'switch'))
0101     d.antenna0.receiver.data(d.flags.switch,:) = nan;
0102   end
0103 end
0104 
0105 disp(['CalculateAlpha:: Selection is ',selection]);
0106 
0107 % Call the classic alpha version:
0108 if strcmp(selection,'CLASSIC')
0109     [A,G,varA,varG,T,r] = DD(d,Ion,Ioff);
0110     disp('CalculateAlpha:: IS CLASSIC');
0111 end
0112     
0113 if strcmp(selection,'FILTERED')
0114     if d.antenna0.receiver.utc(1) < tstr2mjd('01-oct-2011:00:00:00')
0115         disp('CalculateAlpha:: WARN: Cannot apply FILTERED option to data before October 2011')
0116     else
0117         % call the filtered version:
0118         [A,G,varA,varG] = filtered(d,Ion,Ioff);
0119     end
0120  disp('CalculateAlpha:: IS FILTERED');
0121 end
0122 
0123 if strcmp(selection,'POLONLY')
0124     [A,G,varA,varG] = PolOnly(d,Ion,Ioff);
0125     disp('CalculateAlpha:: IS POLONLY');
0126 end
0127 if ~exist('T','var')
0128     T = zeros(length(t),2);
0129     r = zeros(length(t),2);
0130 end
0131 
0132 if (debug)
0133  disp('CalculateAlpha:: Alpha values...');
0134  disp(A);
0135  disp('CalculateAlpha:: #################################################################');
0136  disp('CalculateAlpha:: Gain values...');
0137  disp(G);   
0138 end
0139 
0140 end
0141 
0142 function [A,G,var_A,var_G,T,r] = DD(d,Ion,Ioff)
0143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0144 % [A,G,var_A,var_G,T,r] = DD(d,Ion,Ioff)
0145 % This function calculates the alpha leakage terms for the data in d. It
0146 % finds all the noise diode off/on times and performs the necessary
0147 % calculations. It returns:
0148 %   t - N by 1 time vector in MJD
0149 %   A - N by 5 array of angle values
0150 %   G - N by 5 array of Gain factors: these convert digital units to
0151 %       Kelvin in antenna temperature. It assumes that T_noise_diode = 1 K
0152 %   T - N by 2 array of estimate of offset
0153 %   r - N by 2 array of measured r factors
0154 %   var_A = N by 5 array of variances of A values
0155 %   var_G = N by 5 array of variances of G values
0156 %
0157 % OGK, 12 Mar 2012
0158 %
0159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0160 
0161 if size(d.antenna0.receiver.data,2) ~= 10
0162     error('CalculateAlpha::DD:: Run assembleAlphaStreams first with "CLASSIC" command!');
0163 end
0164 
0165 N = length(Ion);
0166 A = zeros(N,5);
0167 G = zeros(N,5);
0168 var_A = zeros(N,5);
0169 var_G = zeros(N,5);
0170 T = zeros(N,2);
0171 r = zeros(N,2);
0172 
0173 swd = d.antenna0.receiver.data;
0174 
0175 % Work out the parameters for each noise diode on/off event
0176 for k=1:N
0177     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0178     % Work out the correction factors, classic case
0179     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0180     % Total intensity 1:
0181     [Ak,var_Ak,Gk,var_Gk,rk,Tk] = getIparameters_classic(swd(:,1),swd(:,2),Ion{k},Ioff{k});
0182     A(k,1) = Ak; var_A(k,1) = var_Ak;
0183     G(k,1) = Gk; var_G(k,1) = var_Gk;
0184     r(k,1) = rk; T(k,1) = Tk;
0185     % Polarization 1:
0186     [Ak,var_Ak,Gk,var_Gk] = getPparameters_classic(swd(:,3),swd(:,4),Ion{k},Ioff{k});
0187     A(k,2) = Ak; var_A(k,2) = var_Ak;
0188     G(k,2) = Gk; var_G(k,2) = var_Gk;
0189     % Polarization 2:
0190     [Ak,var_Ak,Gk,var_Gk] = getPparameters_classic(swd(:,5),swd(:,6),Ion{k},Ioff{k});
0191     A(k,3) = Ak; var_A(k,3) = var_Ak;
0192     G(k,3) = Gk; var_G(k,3) = var_Gk;
0193     % Polarization 3:
0194     [Ak,var_Ak,Gk,var_Gk] = getPparameters_classic(swd(:,7),swd(:,8),Ion{k},Ioff{k});
0195     A(k,4) = Ak; var_A(k,4) = var_Ak;
0196     G(k,4) = Gk; var_G(k,4) = var_Gk;
0197     % Total intensity 2:
0198     [Ak,var_Ak,Gk,var_Gk,rk,Tk] = getIparameters_classic(swd(:,9),swd(:,10),Ion{k},Ioff{k});
0199     A(k,5) = Ak; var_A(k,5) = var_Ak;
0200     G(k,5) = Gk; var_G(k,5) = var_Gk;
0201     r(k,2) = rk; T(k,2) = Tk;  
0202 end
0203 
0204 end
0205 
0206 function [A,G,varA,varG] = filtered(d,Ion,Ioff)
0207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0208 % [A,G,varA,varG] = filtered(d,Ion,Ioff)
0209 % This function calculates the alpha terms for the data in d. It
0210 % finds all the noise diode off/on times and performs the necessary
0211 % calculations. It returns:
0212 %   t - N by 1 time vector in MJD
0213 %   A - N by 1 array of angle values
0214 %   G - N by 3 array of Gain factors: these convert digital units to
0215 %       Kelvin in antenna temperature. It assumes that T_noise_diode = 1 K
0216 %
0217 % OGK, 13 Mar 2012
0218 %
0219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0220 
0221 if size(d.antenna0.receiver.data,2) ~= 8
0222     error('CalculateAlpha::filtered:: Run assembleAlphaStreams first with "FILTERED" command!');
0223 end
0224 
0225 N = length(Ion);
0226 A = zeros(N,3);
0227 varA = zeros(N,3);
0228 G = zeros(N,5);
0229 varG = zeros(N,5);
0230 data = d.antenna0.receiver.data;
0231 
0232 % Work out the parameters for each noise diode on/off event
0233 for k=1:N
0234     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0235     % Work out the correction factors
0236     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0237     % total intensity 1:
0238     [Gk,varGk] = getIparameters_filtered(data(:,1),Ion{k},Ioff{k});
0239     G(k,1) = Gk;
0240     varG(k,1) = varGk;
0241     % polarization set 1:
0242     [Ak,varAk,Gk,varGk] = getPparameters_filtered(data(:,2),data(:,3),Ion{k},Ioff{k});
0243     G(k,2) = Gk;
0244     varG(k,2) = varGk;
0245     A(k,1) = Ak;
0246     varA(k,1) = varAk;
0247     % polarization set 2:
0248     [Ak,varAk,Gk,varGk] = getPparameters_filtered(data(:,4),data(:,5),Ion{k},Ioff{k});
0249     G(k,3) = Gk;
0250     varG(k,3) = varGk;
0251     A(k,2) = Ak;
0252     varA(k,2) = varAk;
0253     % polarization set 3:
0254     [Ak,varAk,Gk,varGk] = getPparameters_filtered(data(:,6),data(:,7),Ion{k},Ioff{k});
0255     G(k,4) = Gk;
0256     varG(k,4) = varGk;
0257     A(k,3) = Ak;
0258     varA(k,3) = varAk;
0259     % total intensity 2:
0260     [Gk,varGk] = getIparameters_filtered(data(:,8),Ion{k},Ioff{k});
0261     G(k,5) = Gk;
0262     varG(k,5) = varGk;
0263 end
0264 end
0265 
0266 function [A,G,varA,varG] = PolOnly(d,Ion,Ioff)
0267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0268 % [A,G,varA,varG] = PolOnly(d,Ion,Ioff)
0269 % This function calculates the alpha leakage terms for the data in d. It
0270 % finds all the noise diode off/on times and performs the necessary
0271 % calculations. It returns:
0272 %   A - N by 5 array of angle values
0273 %   G - N by 5 array of Gain factors: these convert digital units to
0274 %       Kelvin in antenna temperature. It assumes that T_noise_diode = 1 K
0275 %
0276 % This method differs from all the others in that it works on the
0277 % difference of the diode outputs. This is the correct way of doing this:
0278 % it gets rid of common-mode signals, like 60 Hz pickup. However, this
0279 % makes it impossible to apply an r-factor correction.
0280 %
0281 % OGK, 13 Mar 2012
0282 %
0283 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0284 
0285 if size(d.antenna0.receiver.data,2) ~= 8
0286     error('CalculateAlpha::PolOnly:: Run assembleAlphaStreams first with "POLONLY" flag!');
0287 end
0288 N = length(Ion);
0289 A = zeros(N,3);
0290 G = zeros(N,5);
0291 varA = zeros(N,3);
0292 varG = zeros(N,5);
0293 data = d.antenna0.receiver.data;
0294 % Work out the parameters for each noise diode on/off event
0295 for k=1:N
0296     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0297     % Work out the correction factors
0298     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0299     % total intensity behavior is the same as that of the filtered
0300     % channels.
0301     % Total intensity 1:
0302     [Gk,varGk] = getIparameters_filtered(data(:,1),Ion{k},Ioff{k});
0303     G(k,1) = Gk;
0304     varG(k,1) = varGk;
0305 
0306     % The polarization behavior is the same as the classic data stream.
0307     % Polarization set 1:
0308     [Ak,var_Ak,Gk,var_Gk] = getPparameters_classic(data(:,2),data(:,3),Ion{k},Ioff{k});
0309     A(k,1) = Ak;
0310     varA(k,1) = var_Ak;
0311     G(k,2) = Gk;
0312     varG(k,2) = var_Gk;
0313     
0314     % Polarization set 2:
0315     [Ak,var_Ak,Gk,var_Gk] = getPparameters_classic(data(:,4),data(:,5),Ion{k},Ioff{k});
0316     A(k,2) = Ak;
0317     varA(k,2) = var_Ak;
0318     G(k,3) = Gk;
0319     varG(k,3) = var_Gk;
0320     
0321     % Polarization set 3:
0322     [Ak,var_Ak,Gk,var_Gk] = getPparameters_classic(data(:,6),data(:,7),Ion{k},Ioff{k});
0323     A(k,3) = Ak;
0324     varA(k,3) = var_Ak;
0325     G(k,4) = Gk;
0326     varG(k,4) = var_Gk;
0327     
0328     % Total intensity 2:
0329     [Gk,varGk] = getIparameters_filtered(data(:,8),Ion{k},Ioff{k});
0330     G(k,5) = Gk;
0331     varG(k,5) = varGk;
0332 end
0333 
0334 end
0335 
0336 
0337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0338 % The workhorse functions which calculate the alpha parameters and their
0339 % variances.
0340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0341 
0342 function [A,var_A,G,var_G,r,T] = getIparameters_classic(P1,P2,Ion,Ioff)
0343     % implements the noise diode memo for the classic data stream,
0344     % equations related to total intensity
0345 
0346     dP1 = nanmean(P1(Ion))-nanmean(P1(Ioff));
0347     dP2 = nanmean(P2(Ion))-nanmean(P2(Ioff));
0348     var_dP1 = nanvar(P1(Ion))+nanvar(P1(Ioff));
0349     var_dP2 = nanvar(P2(Ion))+nanvar(P2(Ioff));
0350 
0351     A = (dP1/dP2 - 1)/(dP1/dP2 + 1);
0352     var_A = 4*(dP2^2*var_dP1^2+dP1^2*var_dP2^2)/(dP1+dP2)^4;
0353     C = 1/2*((P1+P2) + (P1-P2)/A); % propto Tsky
0354     D = 1/2*((P1+P2) - (P1-P2)/A); % propto Tload
0355     G = 1/(2*(nanmean(C(Ion))-nanmean(C(Ioff))));
0356     var_G = (nanvar(P1(Ion)+P2(Ion))+nanvar(P1(Ioff)+P2(Ioff)))/...
0357       (4*(nanmean(P1(Ion)+P2(Ion))-nanmean(P1(Ioff)+P2(Ioff)))^4);
0358     r = nanmean(C(Ioff))/nanmean(D(Ioff));
0359     T = 1/(nanmean(C(Ion))/nanmean(C(Ioff)) - 1);
0360 end
0361 
0362 function [G,varG] = getIparameters_filtered(P1,Ion,Ioff)
0363     % implements the noise diode memo for the filtered data stream,
0364     % equations related to total intensity
0365     dP1 = nanmean(P1(Ion))-nanmean(P1(Ioff));
0366     G = 1/(2*dP1);
0367     varG = (nanvar(P1(Ion))+var(P1(Ioff)))/(4*dP1^4);
0368 end
0369 
0370 function [A,var_A,G,var_G] = getPparameters_classic(P1,P2,Ion,Ioff)
0371     % implements the noise diode memo for the classic data stream,
0372     % equations related to polarization
0373     % Polarization set:
0374     dP1 = nanmean(P1(Ion))-nanmean(P1(Ioff)); % cos(th)Tnd
0375     dP2 = nanmean(P2(Ion))-nanmean(P2(Ioff)); % sin(th)Tnd
0376     A = atan2(dP1,dP2);
0377     var_1 = nanvar(P1(Ion))+nanvar(P1(Ioff));
0378     var_2 = nanvar(P2(Ion))+nanvar(P2(Ioff));
0379     var_A = (dP2^2*var_1+dP1^2*var_2)/(dP1^2+dP2^2)^2;
0380     G = 1/sqrt(dP1^2+dP2^2);
0381     var_G = (dP1^2*var_1+dP2^2*var_2)/(dP1^2+dP2^2)^3;
0382 end
0383 
0384 function [A,var_A,G,var_G] = getPparameters_filtered(P1,P2,Ion,Ioff)
0385     % implements the noise diode memo for the filtered data stream,
0386     % equations related to polarization
0387     % Polarization set:
0388     dP1 = nanmean(P1(Ion))-nanmean(P1(Ioff)); % cos(th)Tnd
0389     dP2 = nanmean(P2(Ion))-nanmean(P2(Ioff)); % sin(th)Tnd
0390     
0391     A = -atan2(dP1,dP2);
0392     var_1 = nanvar(P1(Ion))+nanvar(P1(Ioff));
0393     var_2 = nanvar(P2(Ion))+nanvar(P2(Ioff));
0394     var_A = (dP2^2*var_1+dP1^2*var_2)/(dP1^2+dP2^2)^2;
0395     G = 1/sqrt(dP1^2+dP2^2);
0396     var_G = (dP1^2*var_1+dP2^2*var_2)/(dP1^2+dP2^2)^3;
0397 end
0398