rFactorCorrection

 This function calculates the r-factor correction to the data in structure d.
 It uses the methods outlined in Mennella et al (2003), A&A, Vol 410, pp
 1089.
 Method = 1: r = <Sky>/<Load>
 Method = 2: r = rms(Sky)/rms(Load)
 Method = 3: r = optimization of r until 1/f is lowest

 It processes the data in d.antenna0.receiver.switchData, and calculates
 the r-factors for each diode. r is a 4 by 1 vector, with r(1) = diode 1,
 r(2) = diode 2, r(3) = diode 11, and r(4) = diode 12.
 

 Oliver King

0001 function [r,A] = rFactorCorrection(d,method)
0002 % This function calculates the r-factor correction to the data in structure d.
0003 % It uses the methods outlined in Mennella et al (2003), A&A, Vol 410, pp
0004 % 1089.
0005 % Method = 1: r = <Sky>/<Load>
0006 % Method = 2: r = rms(Sky)/rms(Load)
0007 % Method = 3: r = optimization of r until 1/f is lowest
0008 %
0009 % It processes the data in d.antenna0.receiver.switchData, and calculates
0010 % the r-factors for each diode. r is a 4 by 1 vector, with r(1) = diode 1,
0011 % r(2) = diode 2, r(3) = diode 11, and r(4) = diode 12.
0012 %
0013 %
0014 % Oliver King
0015 
0016 % if method == 1
0017 %     disp('Warning: method 1 is not implemented yet. Defaulting to RMS method')
0018 %     method = 2;
0019 % end
0020 
0021 r = zeros(4,1);
0022 A = zeros(4,1);
0023 
0024 % Table of what is being seen in each diode state:
0025 % Column  Diode   State   Source
0026 % 1       1       1       Load
0027 % 2       1       2       Sky
0028 % 3       2       1       Sky
0029 % 4       2       2       Load
0030 % 21      11      1       Load
0031 % 22      11      2       Sky
0032 % 23      12      1       Sky
0033 % 24      12      2       Load
0034 %
0035 
0036 d1L = d.antenna0.receiver.switchData(:,1);
0037 d1S = d.antenna0.receiver.switchData(:,2);
0038 d2S = d.antenna0.receiver.switchData(:,3);
0039 d2L = d.antenna0.receiver.switchData(:,4);
0040 d11L = d.antenna0.receiver.switchData(:,21);
0041 d11S = d.antenna0.receiver.switchData(:,22);
0042 d12S = d.antenna0.receiver.switchData(:,23);
0043 d12L = d.antenna0.receiver.switchData(:,24);
0044         
0045 switch method
0046     case 1
0047         % r = mean(sky)/mean(load)
0048         r(1) = mean(d1S)/mean(d1L);
0049         r(2) = mean(d2S)/mean(d2L);
0050         r(3) = mean(d11S)/mean(d11L);
0051         r(4) = mean(d12S)/mean(d12L);
0052         
0053     case 2
0054         % r = rms(sky)/rms(load)
0055         r(1) = std(d1S-mean(d1S))/std(d1L-mean(d1L));
0056         r(2) = std(d2S-mean(d2S))/std(d2L-mean(d2L));
0057         r(3) = std(d11S-mean(d11S))/std(d11L-mean(d11L));
0058         r(4) = std(d12S-mean(d12S))/std(d12L-mean(d12L));
0059         
0060     case 3
0061         % r is obtained by minimizing 1/f knee frequency of sky - r*load
0062         % Do for each diode in turn
0063         fk = zeros(4,1);
0064         disp('r   =>  fknee:')
0065         [a,b] = findr(d1L,d1S);
0066         r(1) = a; fk(1) = b;
0067         fprintf('%3.2f  =>  %4.3f\n',r(1),fk(1))
0068         [a,b] = findr(d2L,d2S);
0069         r(2) = a; fk(2) = b;
0070         fprintf('%3.2f  =>  %4.3f\n',r(2),fk(2))
0071         [a,b] = findr(d11L,d11S);
0072         r(3) = a; fk(3) = b;
0073         fprintf('%3.2f  =>  %4.3f\n',r(3),fk(3))
0074         [a,b] = findr(d12L,d12S);
0075         r(4) = a; fk(4) = b;
0076         fprintf('%3.2f  =>  %4.3f\n',r(4),fk(4))
0077      case 4
0078         % r is obtained by minimizing 1/f knee frequency of sky - r*load
0079         % Do for each diode in turn
0080         fk = zeros(4,1);
0081         disp('r   =>  fknee:')
0082         [a,b] = findRAndA(d1L,d1S);
0083         r(1) = a(1); fk(1) = b; A(1) = a(2);
0084         fprintf('%3.2f  =>  %4.3f\n',r(1),fk(1))
0085         [a,b] = findRAndA(d2L,d2S);
0086         r(2) = a(1); fk(2) = b; A(2) = a(2);
0087         fprintf('%3.2f  =>  %4.3f\n',r(2),fk(2))
0088         [a,b] = findRAndA(d11L,d11S);
0089         r(3) = a(1); fk(3) = b; A(3) = a(2);
0090         fprintf('%3.2f  =>  %4.3f\n',r(3),fk(3))
0091         [a,b] = findRAndA(d12L,d12S);
0092         r(4) = a(1); fk(4) = b; A(4) = a(2);
0093         fprintf('%3.2f  =>  %4.3f\n',r(4),fk(4))
0094 end
0095 
0096 end
0097 
0098 function [r,fknee] = findr(L,S)
0099 % Subfunction which finds the value of r which minimizes the knee frequency
0100 % of the spectrum of S-r*L
0101 c = {S L 100};
0102 options = optimset('MaxFunEvals',100000,'MaxIter',100000,'Algorithm','active-set');%,...
0103 %    'PlotFcns',{@optimplotx,@optimplotfval});
0104 % [x,fval] = fminsearch(@(x) errFunc(x,c),1,options); % initial value for r is 1
0105 [r,fknee] = fmincon(@(r) errFunc(r,c),1,[],[],[],[],0.5,2,[],options); % 0.1 < r < 2
0106 end
0107 
0108 function [x,fknee] = findRAndA(L,S)
0109 % Subfunction which finds the value of r and offset A which minimizes the
0110 % knee frequency of the spectrum of S-r*L
0111 c = {S L 100};
0112 options = optimset('MaxFunEvals',100000,'MaxIter',100000,'Algorithm','active-set',...
0113     'PlotFcns',{@optimplotx,@optimplotfval});
0114 % [x,fval] = fminsearch(@(x) errFunc(x,c),1,options); % initial value for r is 1
0115 [x,fknee] = fmincon(@(x) errFuncRAndA(x,c),[1;10],[],[],[],[],[0.1;-10000],...
0116     [2; 10000],[],options); % 0.1 < r < 2
0117 end
0118 
0119 function fknee = errFunc(r,c)
0120 % error function to minimize
0121 S = c{1};
0122 L = c{2};
0123 V = S-r*L;
0124 fs = c{3};
0125 [a,A,B,fknee] = fit1overf(V,fs,0);
0126 end
0127 
0128 function fknee = errFuncRAndA(x,c)
0129 % error function to minimize
0130 S = c{1};
0131 L = c{2};
0132 V = (S+x(2))-x(1)*(L+x(2));
0133 fs = c{3};
0134 [a,A,B,fknee] = fit1overf(V,fs,0);
0135 end