Home > matutils > analyze_measureRfactor.m

analyze_measureRfactor

PURPOSE ^

[elevation,r] = analyze_measureRfactor(start,finish)

SYNOPSIS ^

function [elevation,r] = analyze_measureRfactor(start,finish)

DESCRIPTION ^

 [elevation,r] = analyze_measureRfactor(start,finish)

 This script takes in the start and finish times (as strings) of the
 measureRfactor.sch schedule, calculates the elevation-dependent
 r-factors, plots the data along with some best-fit curves, and returns
 the calculated r factors.

 OGK, Nov 2010

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [elevation,r] = analyze_measureRfactor(start,finish)
0002 % [elevation,r] = analyze_measureRfactor(start,finish)
0003 %
0004 % This script takes in the start and finish times (as strings) of the
0005 % measureRfactor.sch schedule, calculates the elevation-dependent
0006 % r-factors, plots the data along with some best-fit curves, and returns
0007 % the calculated r factors.
0008 %
0009 % OGK, Nov 2010
0010 %
0011 
0012 d = pipe_read(start,finish);
0013 
0014 % Work with the noise data:
0015 [tA,A,G,T,horiz,equa, offStartPos, onEndPos, offEndPos, onStartPos] = calculateAlpha(d);
0016 
0017 % The noise events last for 66 seconds
0018 
0019 % Apply alpha
0020 t = d.antenna0.receiver.utc;
0021 swd = d.antenna0.receiver.switchData;
0022 
0023 % Matrices of alpha and gain values:
0024 Am = zeros(length(t),8);
0025 Gm = zeros(length(t),8);
0026 
0027 % Get an elevation for each measurement:
0028 el = d.antenna0.servo.apparent(:,2);
0029 elevation = zeros(length(tA),1);
0030 for m=1:length(tA)    
0031     Isubs = abs(t-tA(m))*24*60*60 <= 35;
0032     elevation(m) = mean(el(Isubs));
0033     for k=1:8
0034         Am(Isubs,k) = A(m,k);
0035         Gm(Isubs,k) = G(m,k);
0036     end
0037 end
0038 
0039 % And now form the corrected switchData products for total intensity:
0040 
0041 % Total intensity channel 1
0042 Cs = 1/2*(swd(:,2)+swd(:,1) + 1./(1-2*Am(:,1)).*(swd(:,2)-swd(:,1)));
0043 Dl = 1/2*(swd(:,2)+swd(:,1) - 1./(1-2*Am(:,1)).*(swd(:,2)-swd(:,1)));
0044 d.antenna0.receiver.switchData(:,1) = Dl.*Gm(:,1);
0045 d.antenna0.receiver.switchData(:,2) = Cs.*Gm(:,1);
0046 Cs = 1/2*(swd(:,3)+swd(:,4) + 1./(1-2*Am(:,2)).*(swd(:,3)-swd(:,4)));
0047 Dl = 1/2*(swd(:,3)+swd(:,4) - 1./(1-2*Am(:,2)).*(swd(:,3)-swd(:,4)));
0048 d.antenna0.receiver.switchData(:,3) = Cs.*Gm(:,2);
0049 d.antenna0.receiver.switchData(:,4) = Dl.*Gm(:,2);
0050 
0051 % Total intensity channel 2
0052 Cs = 1/2*(swd(:,22)+swd(:,21) + 1./(1-2*Am(:,7)).*(swd(:,22)-swd(:,21)));
0053 Dl = 1/2*(swd(:,22)+swd(:,21) - 1./(1-2*Am(:,7)).*(swd(:,22)-swd(:,21)));
0054 d.antenna0.receiver.switchData(:,21) = Dl.*Gm(:,7);
0055 d.antenna0.receiver.switchData(:,22) = Cs.*Gm(:,7);
0056 Cs = 1/2*(swd(:,23)+swd(:,24) + 1./(1-2*Am(:,8)).*(swd(:,23)-swd(:,24)));
0057 Dl = 1/2*(swd(:,23)+swd(:,24) - 1./(1-2*Am(:,8)).*(swd(:,23)-swd(:,24)));
0058 d.antenna0.receiver.switchData(:,23) = Cs.*Gm(:,8);
0059 d.antenna0.receiver.switchData(:,24) = Dl.*Gm(:,8);
0060 
0061 
0062 % Form the switchData products for the polarization channels:
0063 % Q11, U11 combination:
0064 c5 = swd(:,6)-swd(:,5);
0065 c6 = swd(:,10)-swd(:,9);
0066 d.antenna0.receiver.switchData(:,5) = -Gm(:,3).*cos(Am(:,3)).*c5;
0067 d.antenna0.receiver.switchData(:,6) = -Gm(:,3).*sin(Am(:,3)).*c6;
0068 d.antenna0.receiver.switchData(:,9) = -Gm(:,3).*sin(Am(:,3)).*c5;
0069 d.antenna0.receiver.switchData(:,10) = Gm(:,3).*cos(Am(:,3)).*c6;
0070 
0071 % Form the switchData products for the polarization channels:
0072 % Q12, U12 combination:
0073 c5 = swd(:,8)-swd(:,7);
0074 c6 = swd(:,12)-swd(:,11);
0075 d.antenna0.receiver.switchData(:,8) = -Gm(:,4).*cos(Am(:,4)).*c5;
0076 d.antenna0.receiver.switchData(:,7) = -Gm(:,4).*sin(Am(:,4)).*c6;
0077 d.antenna0.receiver.switchData(:,12) =-Gm(:,4).*sin(Am(:,4)).*c5;
0078 d.antenna0.receiver.switchData(:,11) = Gm(:,4).*cos(Am(:,4)).*c6;
0079 
0080 % Form the switchData products for the polarization channels:
0081 % Q21, U21 combination:
0082 c5 = swd(:,14)-swd(:,13);
0083 c6 = swd(:,18)-swd(:,17);
0084 d.antenna0.receiver.switchData(:,13) = -Gm(:,5).*cos(Am(:,5)).*c5;
0085 d.antenna0.receiver.switchData(:,14) = -Gm(:,5).*sin(Am(:,5)).*c6;
0086 d.antenna0.receiver.switchData(:,17) = -Gm(:,5).*sin(Am(:,5)).*c5;
0087 d.antenna0.receiver.switchData(:,18) =  Gm(:,5).*cos(Am(:,5)).*c6;
0088 
0089 % Form the switchData products for the polarization channels:
0090 % Q22, U22 combination:
0091 c5 = swd(:,16)-swd(:,15);
0092 c6 = swd(:,20)-swd(:,19);
0093 d.antenna0.receiver.switchData(:,16) = -Gm(:,6).*cos(Am(:,6)).*c5;
0094 d.antenna0.receiver.switchData(:,15) = -Gm(:,6).*sin(Am(:,6)).*c6;
0095 d.antenna0.receiver.switchData(:,20) = -Gm(:,6).*sin(Am(:,6)).*c5;
0096 d.antenna0.receiver.switchData(:,19) =  Gm(:,6).*cos(Am(:,6)).*c6;
0097 
0098 % Now calculate the r factor at each event:
0099 N = length(tA);
0100 r = zeros(N,4);
0101 swd = d.antenna0.receiver.switchData;
0102 
0103 % Work out the parameters for each noise diode on/off event
0104 for k=1:N
0105   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0106   % The relevant start and end positions for this noise diode event
0107   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0108 
0109   if d.antenna0.receiver.utc(offStartPos(k)) < tstr2mjd('08-SEP-2010:19:00:00')
0110       
0111     Offs = offStartPos(k);
0112     Offe = offEndPos(k);
0113 
0114     swdoff = swd(Offs:Offe,:);
0115   else
0116     Offs1 = offStartPos(k);
0117     Offe1 = onStartPos(k)-1;
0118     Offs2 = onEndPos(k)+1;
0119     Offe2 = offEndPos(k);
0120 
0121     swdoff = cat(1,swd(Offs1:Offe1,:),swd(Offs2:Offe2,:));
0122   end
0123   
0124   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0125   % Work out the r-factor corrections
0126   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0127     r(k,1) = nanmean(swdoff(:,2))/nanmean(swdoff(:,1));
0128     r(k,2) = nanmean(swdoff(:,3))/nanmean(swdoff(:,4));
0129     r(k,3) = nanmean(swdoff(:,22))/nanmean(swdoff(:,21));
0130     r(k,4) = nanmean(swdoff(:,23))/nanmean(swdoff(:,24));
0131 end
0132 
0133 atmosphere = @(x)(x(1)*(1./cos((90-elevation)*pi/180)-1)+x(2));
0134 errsq1 = @(x)(sum( ( r(:,1) - atmosphere(x) ).^2 ));
0135 errsq2 = @(x)(sum( ( r(:,2) - atmosphere(x) ).^2 ));
0136 errsq3 = @(x)(sum( ( r(:,3) - atmosphere(x) ).^2 ));
0137 errsq4 = @(x)(sum( ( r(:,4) - atmosphere(x) ).^2 ));
0138 
0139 % Use fminsearch to find the best fits for all four functions
0140 [x1,fval1] = fminsearch(errsq1,[0.1, 1.5]);
0141 [x2,fval2] = fminsearch(errsq2,[0.1, 1.5]);
0142 [x3,fval3] = fminsearch(errsq3,[0.1, 1.5]);
0143 [x4,fval4] = fminsearch(errsq4,[0.1, 1.5]);
0144 
0145 % The mean of the elevations above 37 degrees:
0146 mn = mean(r(elevation >= 37,:),1);
0147 
0148 figure
0149 subplot(2,2,1)
0150 plot(elevation,r(:,1),'k.',elevation,atmosphere(x1),'r-',elevation(elevation>=37),mn(1)*ones(size(elevation(elevation>=37))),'b-')
0151 ylabel('Diode 1 r-Factor')
0152 legend('Data',sprintf('r = %4.3f [sec(za) - 1] + %4.3f',x1(1),x1(2)),sprintf('r = %4.3f, mean above 37 deg',mn(1)))
0153 subplot(2,2,2)
0154 plot(elevation,r(:,2),'k.',elevation,atmosphere(x2),'r-',elevation(elevation>=37),mn(2)*ones(size(elevation(elevation>=37))),'b-')
0155 ylabel('Diode 2 r-Factor')
0156 legend('Data',sprintf('r = %4.3f [sec(za) - 1] + %4.3f',x2(1),x2(2)),sprintf('r = %4.3f, mean above 37 deg',mn(2)))
0157 subplot(2,2,3)
0158 plot(elevation,r(:,3),'k.',elevation,atmosphere(x3),'r-',elevation(elevation>=37),mn(3)*ones(size(elevation(elevation>=37))),'b-')
0159 xlabel('Elevation [deg]')
0160 legend('Data',sprintf('r = %4.3f [sec(za) - 1] + %4.3f',x3(1),x3(2)),sprintf('r = %4.3f, mean above 37 deg',mn(3)))
0161 ylabel('Diode 3 r-Factor')
0162 subplot(2,2,4)
0163 plot(elevation,r(:,4),'k.',elevation,atmosphere(x4),'r-',elevation(elevation>=37),mn(4)*ones(size(elevation(elevation>=37))),'b-')
0164 legend('Data',sprintf('r = %4.3f [sec(za) - 1] + %4.3f',x4(1),x4(2)),sprintf('r = %4.3f, mean above 37 deg',mn(4)))
0165 xlabel('Elevation [deg]')
0166 ylabel('Diode 4 r-Factor')
0167 
0168 end
Generated on Sun 14-Jun-2015 17:12:45 by m2html © 2005