Description of load

0001 function [A B C] = load_measure(d,fitstep,fitlen)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % function [A B C] = load_measure(d,fitstep,fitlen)
0006 %
0007 %   Calculates the A, B, and C parameters of the ccTemperatureLoad or ccColdPlate vector.
0008 %   These are used for applying the cold load templates to the data
0009 %   vectors.
0010 %
0011 %   I/O:
0012 %    - 'd' is the data structure to which we wish to apply the template.
0013 %    - 'fitstep' is the stepsize (in minutes) at which A, B, and C will be
0014 %    calculated.  I recommend 0.5.
0015 %    - 'fitlen' is the amount of data (in minutes) which will be used in
0016 %    the calculation of a single set of A, B, and C values.  I recommend 5.
0017 %
0018 %   MAS -- 29-March-2012
0019 %
0020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 
0022 
0023 % Convert the fitting parameters to number of 5Hz samples.
0024 fitStepS = fitstep * 60 * 5;
0025 fitLenS = fitlen * 60 * 5;
0026 
0027 
0028 % Get the cold load vector.
0029 if d.antenna0.thermal.utc(end) < tstr2mjd('30-Jul-2012:00:00:00')
0030     notNAN = ~isnan(d.antenna0.thermal.ccTemperatureLoad);
0031     coldLoad = d.antenna0.thermal.ccTemperatureLoad(notNAN);
0032     coldLoadUtc = d.antenna0.thermal.utc(notNAN);
0033 else
0034     notNAN = ~isnan(d.antenna0.thermal.ccColdPlate);
0035     coldLoad = d.antenna0.thermal.ccColdPlate(notNAN);
0036     coldLoadUtc = d.antenna0.thermal.utc(notNAN);
0037 end
0038 
0039 % Filter the cold load vector to remove sub-Hz variations.
0040 [home,installeddir] = where_am_i();
0041 f = open([home '/' installeddir '/reduc/load/filters/load_filter_1.2.mat']);
0042 filtLoad = filter(f.Num, 1, coldLoad);
0043 fLen = length(f.Num);
0044 
0045 coldLoad  = filtLoad(fLen:end);
0046 coldLoadUtc  = coldLoadUtc(ceil(fLen/2):end-floor(fLen/2));
0047 
0048 loadLen = length(coldLoad);
0049 
0050 coldLoadX = 2* pi * 86400 * (coldLoadUtc - min(d.antenna0.receiver.utc));
0051 
0052 
0053 % Now we step over the filtered cold load vector, calculating the A, B, and
0054 % C values.
0055 fitLoop = 1;
0056 k = 1;
0057 startPoint = [1e-3 1.2 0];
0058 while (fitLoop == 1)
0059     
0060     startCut = 1 + (k-1) * fitStepS;
0061     endCut = fitLenS + (k-1) * fitStepS;
0062 
0063     if endCut >= loadLen
0064         startCut = loadLen - fitLenS + 1;
0065         endCut = loadLen;
0066 
0067         fitLoop = 0;
0068     end
0069 
0070     if (startCut > (endCut - 10) || startCut < 1)
0071         A(k) = 0;
0072         B(k) = 0;
0073         C(k) = 0;
0074         U(k) = 0;
0075         fFlag(k) = 1;
0076         continue
0077     end
0078     
0079     coldCut = detrend(coldLoad(startCut:endCut));
0080     xCut = coldLoadX(startCut:endCut);
0081     utcCut = coldLoadUtc(startCut:endCut);
0082     
0083     fitLoad = fit(xCut, coldCut, 'sin1', 'Startpoint', startPoint);
0084 
0085     
0086     if fitLoad.a1 < 0
0087         A(k) = abs(fitLoad.a1);
0088         C(k) = mod(fitLoad.c1 + pi, 2*pi);
0089     else
0090         A(k) = fitLoad.a1;
0091         C(k) = mod(fitLoad.c1,2*pi);
0092     end
0093 
0094     B(k) = fitLoad.b1;
0095     U(k) = mean(utcCut);
0096 
0097 
0098     confLoad = confint(fitLoad);
0099     fErr = range(confLoad);
0100 
0101     if (max(fErr(1)/A(k),fErr(2)/B(k)) < 0.1)
0102         startPoint = [A(k) B(k) C(k)];
0103         fFlag(k) = 0;
0104     else
0105         fFlag(k) = 1;
0106     end
0107         
0108     
0109     k = k + 1;
0110 end
0111 
0112 % Unwrap the C (phase) term.
0113 wrLen = find(diff(C) > pi);
0114 for k=1:length(wrLen)
0115     C(wrLen(k)+1:end) = C(wrLen(k)+1:end) - 2*pi;
0116 end
0117 
0118 wrLen = find(diff(C) < -pi);
0119 for k=1:length(wrLen)
0120     C(wrLen(k)+1:end) = C(wrLen(k)+1:end) + 2*pi;
0121 end
0122 
0123 
0124 if sum(fFlag == 0) > 1
0125     % Interpolate and smooth over the full 100Hz time vector.
0126     A = smooth(interp1(U(fFlag == 0), A(fFlag == 0), ...
0127         d.antenna0.receiver.utc, 'nearest', 'extrap'),1.5*fitStepS*20);
0128     B = smooth(interp1(U(fFlag == 0), B(fFlag == 0), ...
0129         d.antenna0.receiver.utc, 'nearest', 'extrap'),1.5*fitStepS*20);
0130     C = smooth(interp1(U(fFlag == 0), C(fFlag == 0), ...
0131         d.antenna0.receiver.utc, 'nearest', 'extrap'),1.5*fitStepS*20);
0132 
0133 %     Aint = smooth(interp1(U(fFlag == 0), A(fFlag == 0), ...
0134 %         d.antenna0.receiver.utc, 'nearest', 'extrap'),1.5*fitStepS*20);
0135 %     Bint = smooth(interp1(U(fFlag == 0), B(fFlag == 0), ...
0136 %         d.antenna0.receiver.utc, 'nearest', 'extrap'),1.5*fitStepS*20);
0137 %     Cint = smooth(interp1(U(fFlag == 0), C(fFlag == 0), ...
0138 %         d.antenna0.receiver.utc, 'nearest', 'extrap'),1.5*fitStepS*20);
0139 %
0140 %     figure;
0141 %     plot(24*60*(d.antenna0.receiver.utc-d.antenna0.receiver.utc(1)),Aint);
0142 %     hold on;
0143 %     plot(24*60*(U-d.antenna0.receiver.utc(1)),A,'r.');
0144 %     hold off;
0145 %     figure;
0146 %     plot(24*60*(d.antenna0.receiver.utc-d.antenna0.receiver.utc(1)),Bint);
0147 %     hold on;
0148 %     plot(24*60*(U-d.antenna0.receiver.utc(1)),B,'r.');
0149 %     hold off;
0150 %     figure;
0151 %     plot(24*60*(d.antenna0.receiver.utc-d.antenna0.receiver.utc(1)),Cint);
0152 %     hold on;
0153 %     plot(24*60*(U-d.antenna0.receiver.utc(1)),C,'r.');
0154 %     hold off;
0155 %
0156 %     A = Aint;
0157 %     B = Bint;
0158 %     C = Cint;
0159 else
0160     A = 1e-3 * ones(length(d.antenna0.receiver.utc),1);
0161     B = 1.2 * ones(length(d.antenna0.receiver.utc),1);
0162     C = zeros(length(d.antenna0.receiver.utc),1);    
0163 end
0164 
0165 end
load_measure

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
