interpRegisters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function d = interpRegisters(d)

  only the cryocon and the az/el need to be interpolated

  Modified:
   27-Jan-2011 (MAS) -- Made the interpolation of corrupted data more
   efficient.
   3-DEC-2012 (MAS) -- Placed the indCryoBad definition inside an 'if'
   statment so that this doesn't crash if the user doesn't have the
   thermal info.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function d = interpRegisters(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %  function d = interpRegisters(d)
0005 %
0006 %  only the cryocon and the az/el need to be interpolated
0007 %
0008 %  Modified:
0009 %   27-Jan-2011 (MAS) -- Made the interpolation of corrupted data more
0010 %   efficient.
0011 %   3-DEC-2012 (MAS) -- Placed the indCryoBad definition inside an 'if'
0012 %   statment so that this doesn't crash if the user doesn't have the
0013 %   thermal info.
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 
0016 hasRx   = 0;
0017 hasCryo = 0;
0018 hasServo= 0;
0019 %display('interpRegisters:: interpRegisters');
0020 if(isfield(d, 'antenna0'))
0021   if(isfield(d.antenna0, 'servo'))
0022     if(isfield(d.antenna0.servo, 'fast_az_pos'))
0023       hasServo = 1;
0024     end
0025   end
0026   
0027   if(isfield(d.antenna0, 'thermal'))
0028     if(isfield(d.antenna0.thermal, 'ccTemperatureLoad'))
0029       hasCryo = 1;
0030     end
0031   end
0032   
0033   if(isfield(d.antenna0, 'receiver'))
0034     if(isfield(d.antenna0.receiver, 'utc'))
0035       hasRx = 1;
0036     end
0037   end
0038 end
0039 
0040 %display('interpRegisters:: interpRegisters1');
0041 hasAny = hasServo + hasRx + hasCryo;
0042 if(hasAny<2)
0043   %display('interpRegisters:: Nothing to inpterpolate')
0044   return;
0045 end
0046 
0047 if(~hasRx)
0048   display('interpRegisters:: Nothing to interpolate')
0049   return;
0050 end
0051 
0052 
0053 utcLength = length(d.antenna0.receiver.utc);
0054 
0055 %display('interpRegisters:: interpRegisters2');
0056 % check to see that the backend reported back data on time for the whole
0057 % track -- at this point we know it has receiver data.
0058 indRxBad = zeros(size(d.antenna0.receiver.utc));
0059 if(utcLength~=length(unique(d.antenna0.receiver.utc)))
0060 %display('interpRegisters:: interpRegisters3');
0061   % look for where the frame was dropped
0062   a = sort(d.antenna0.receiver.utc);
0063   % Values with utc == 0 are flagged later.
0064   % Do not bother putting these through the for loop. (MAS  27-Jan-2011)
0065   f = diff(a)==0 & a(2:end) ~= 0;  
0066   
0067   % Just get the unique utc values.  The old way repeated values in the
0068   % following loop.  (MAS  27-Jan-2011)
0069   aUnique = unique(a(f));
0070 
0071   % For loop modified to handle previous two changes.  Should be much more
0072   % efficient.  (MAS  27-Jan-2011)
0073   for m=1:length(aUnique)
0074     ff = find(d.antenna0.receiver.utc == aUnique(m));
0075     indRxBad(ff(2:length(ff))) = 1;
0076     %display('interpRegisters:: interpRegisters4');
0077     %length(ff)
0078     %m
0079   end    
0080 % display('interpRegisters:: interpRegisters5');
0081   % check for zero values.
0082   indRxBad = indRxBad | d.antenna0.receiver.utc==0;
0083   
0084   % Now find the flagged regions.  If several flagged points are adjacent,
0085   % then this makes the following for loop much more efficient.
0086   % (MAS  27-Jan-2011)
0087   indRxBadStart = find(diff(indRxBad) > 0) + 1;
0088   indRxBadEnd = find(diff(indRxBad) < 0);
0089   
0090   if indRxBad(1) == 1
0091     indRxBadStart = cat(1,1,indRxBadStart);      
0092   end
0093   if indRxBad(end) == 1
0094     indRxBadEnd = cat(1,indRxBadEnd,utcLength);      
0095   end
0096   
0097 
0098   % fix the times in the servo and keep track of bad values.
0099   ff = find(indRxBad);
0100   timeDiff = median(diff(d.antenna0.receiver.utc(~indRxBad)));
0101 
0102   % New for loop.  Designed to be much more efficient in the case where
0103   % several adjacent times are flagged.  (MAS  27-Jan-2011)
0104   for m=1:length(indRxBadStart)
0105     interpLength = indRxBadEnd(m) - indRxBadStart(m) + 1;
0106       
0107     if indRxBadStart(m) == 1 && indRxBadEnd(m) ~= utcLength
0108       d.antenna0.receiver.utc(indRxBadStart(m):indRxBadEnd(m)) = ...
0109           d.antenna0.receiver.utc(indRxBadEnd(m)+1) ...
0110           - timeDiff * (interpLength:-1:1)';
0111     elseif indRxBadStart(m) == 1 && indRxBadEnd(m) == utcLength
0112       display('interpRegisters:: ALL OF YOUR BACKEND DATA ARE BAD');
0113       display('interpRegisters:: DO NOT RELY ON THESE DATA FOR MAPMAKING');
0114       d.antenna0.receiver.utc = nan(length(d.antenna0.receiver.utc), 1);
0115     else
0116       d.antenna0.receiver.utc(indRxBadStart(m):indRxBadEnd(m)) = ...
0117           d.antenna0.receiver.utc(indRxBadStart(m)-1) ...
0118           + timeDiff * (1:interpLength)';
0119     end
0120   end
0121   
0122 %  for m=1:length(ff)
0123 %    d.antenna0.receiver.utc(ff(m)) = d.antenna0.receiver.utc(ff(m)-1) + timeDiff;
0124 %  end
0125 %display('interpRegisters:: interpRegisters5');
0126   display('interpRegisters:: Number of bad backend times in track:');
0127   display(sprintf('interpRegisters:: %i seconds out of %i seconds', length(ff)/100, ...
0128       length(indRxBad)/100));
0129   
0130 end    
0131 
0132 % first thing we have to do is check that the thermal and cryo times did not
0133 % skip -- this happens in the servo's case when the 1PPS tick does not match
0134 % the internal clock, and for the cryo when the gpib doesn't report back on
0135 % time.
0136 if(hasServo)
0137   indServoBad = zeros(size(d.antenna0.servo.utc));
0138   %% in some rare cases, the servo reports a few bad data points.  we can ...
0139       %% check as follows:
0140   az5hz = interp1(d.array.frame.utc-0.5./60/60/24, d.antenna0.servo.slow_az_pos, d.antenna0.servo.utc);
0141   azdiff = abs(az5hz - d.antenna0.servo.fast_az_pos);
0142   indServoBad =  azdiff > 5 & azdiff < 355;
0143   indServoBad =  indServoBad | azdiff> 364;
0144   [s e] = findStartStop(indServoBad);
0145   %% issues with the timing are usually on 5 sample increments.
0146   for m=1:length(s)
0147     if(mod(s(m),5)~=1)
0148       aa = mod(s(m),5);
0149       s(m) = s(m) - aa + 1;
0150     end
0151     if(mod(e(m),5)~=0)
0152       aa = mod(e(m),5);
0153       e(m) = e(m) + 5 - aa;
0154     end
0155     if(s(m) < 1)
0156       s(m) = 1;
0157     end
0158     if(e(m) > length(indServoBad))
0159       e(m) = length(indServoBad);
0160     end
0161   end
0162 else
0163   indServoBad = [];
0164 end
0165 
0166 if(hasServo)
0167 
0168   utcLength = length(d.antenna0.servo.utc);
0169 
0170   if(utcLength~=length(unique(d.antenna0.servo.utc)))
0171     % this happens when the servo reports back the wrong number of samples,
0172     % typically due to a 1PPS error.
0173     a = sort(d.antenna0.servo.utc);
0174     % Values with utc == 0 are flagged later.
0175     % Don not bother putting these through the for loop. (MAS  27-Jan-2011)
0176     f = diff(a)==0 & a(2:end) ~= 0;  
0177 
0178 
0179     % Just get the unique utc values.  The old way repeated values in the
0180     % following loop.  (MAS  27-Jan-2011)
0181     aUnique = unique(a(f));
0182 
0183     % For loop modified to handle previous two changes.  Should be much more
0184     % efficient.  (MAS  27-Jan-2011)
0185     for m=1:length(aUnique)
0186       ff = find(d.antenna0.servo.utc == aUnique(m));
0187       indServoBad(ff(2:length(ff))) = 1;
0188     end    
0189 
0190     % check for zero values in the time.
0191     indServoBad = indServoBad | d.antenna0.servo.utc==0;
0192   
0193     
0194     % Now find the flagged regions.  If several flagged points are adjacent,
0195     % then this makes the following for loop much more efficient.
0196     % (MAS  27-Jan-2011)
0197     
0198     indServoBadStart = find(diff(indServoBad) > 0) + 1;
0199     indServoBadEnd = find(diff(indServoBad) < 0);
0200   
0201     if indServoBad(1) == 1
0202       indServoBadStart = cat(1,1,indServoBadStart);      
0203     end
0204     if indServoBad(end) == 1
0205       indServoBadEnd = cat(1,indServoBadEnd,utcLength);      
0206     end
0207   
0208 
0209     % fix the times in the servo and keep track of bad values.
0210     timeDiff = median(diff(d.antenna0.servo.utc(~indServoBad)));
0211 
0212     % New for loop.  Designed to be much more efficient in the case where
0213     % several adjacent times are flagged.  (MAS  27-Jan-2011)
0214     for m=1:length(indServoBadStart)
0215       interpLength = indServoBadEnd(m) - indServoBadStart(m) + 1;
0216       
0217       if indServoBadStart(m) == 1 && indServoBadEnd(m) ~= utcLength
0218         d.antenna0.servo.utc(indServoBadStart(m):indServoBadEnd(m)) = ...
0219             d.antenna0.servo.utc(indServoBadEnd(m)+1) ...
0220             - timeDiff * (interpLength:-1:1)';
0221       elseif indServoBadStart(m) == 1 && indServoBadEnd(m) == utcLength
0222         display('interpRegisters:: ALL OF YOUR SERVO DATA ARE BAD');
0223         display('interpRegisters:: DO NOT RELY ON THESE DATA FOR MAPMAKING');
0224         d.antenna0.servo.utc = nan(utcLength, 1);
0225       else
0226         d.antenna0.servo.utc(indServoBadStart(m):indServoBadEnd(m)) = ...
0227             d.antenna0.servo.utc(indServoBadStart(m)-1) ...
0228             + timeDiff * (1:interpLength)';
0229       end
0230     end
0231 
0232   end  
0233   
0234 
0235 %    for m=1:length(ff);
0236 %      d.antenna0.servo.utc(ff(m)) = d.antenna0.servo.utc(ff(m)-1) + timeDiff;
0237 %    end
0238 
0239   ff = find(indServoBad);
0240 
0241   if length(ff) ~= utcLength
0242 
0243     % next we interpolate onto the new times, and keep track of them to flag.
0244     utc    = d.antenna0.servo.utc;
0245     az_pos = d.antenna0.servo.fast_az_pos;  
0246     el_pos = d.antenna0.servo.fast_el_pos;  
0247     az_err = d.antenna0.servo.fast_az_err;
0248     el_err = d.antenna0.servo.fast_el_err;
0249     utc(ff)    = [];
0250     az_pos(ff) = [];
0251     az_err(ff) = [];
0252     el_pos(ff) = [];
0253     el_err(ff) = [];
0254   
0255     d.antenna0.servo.fast_az_pos = interp1(utc, az_pos, d.antenna0.servo.utc);
0256     d.antenna0.servo.fast_az_err = interp1(utc, az_err, d.antenna0.servo.utc);
0257     d.antenna0.servo.fast_el_pos = interp1(utc, el_pos, d.antenna0.servo.utc);  
0258     d.antenna0.servo.fast_el_err = interp1(utc, el_err, d.antenna0.servo.utc);  
0259   else
0260       
0261     % In this case, these data are borked.  Just put something here...
0262     d.antenna0.servo.fast_az_pos = nan(utcLength,1);
0263     d.antenna0.servo.fast_az_err = nan(utcLength,1);
0264     d.antenna0.servo.fast_el_pos = nan(utcLength,1);  
0265     d.antenna0.servo.fast_el_err = nan(utcLength,1);        
0266   end
0267     
0268   
0269   display('interpRegisters:: Number of bad servo times in track (due to noisy 1PPS):');
0270   display(sprintf('interpRegisters:: %i seconds out of %i seconds', length(ff)/5, length(indServoBad)/5));
0271 
0272 end
0273 
0274 if(hasCryo)
0275 
0276   %This goes inside the if statement, in case the user has not requested
0277   % the thermal data.  MAS.  3-Dec-2012.
0278   indCryoBad = zeros(size(d.antenna0.thermal.utc));
0279     
0280     
0281   if round(length(d.antenna0.receiver.utc) / length(d.antenna0.thermal.utc)) ~= 20
0282     d.antenna0.thermal.utc = d.antenna0.receiver.utc(1:20:end);
0283     indCryoBad = zeros(size(d.antenna0.thermal.utc));
0284   end
0285   if round(length(d.antenna0.receiver.utc) / length(d.antenna0.thermal.ccTemperatureLoad)) ~= 20
0286     d.antenna0.thermal.ccTemperatureLoad = ...
0287         interp1( ...
0288         d.antenna0.receiver.utc(1:round(length(d.antenna0.receiver.utc)/length(d.antenna0.thermal.ccTemperatureLoad)):end), ...
0289         d.antenna0.thermal.ccTemperatureLoad,d.antenna0.thermal.utc,'nearest','extrap');
0290   end
0291     
0292   utcLength = length(d.antenna0.thermal.utc);
0293 
0294   if(utcLength~=length(unique(d.antenna0.thermal.utc)))
0295     % this happens when the servo reports back the wrong number of samples,
0296     % typically due to a 1PPS error.
0297     a = sort(d.antenna0.thermal.utc);
0298     % Values with utc == 0 are flagged later.
0299     % Do not bother putting these through the for loop. (MAS  27-Jan-2011)
0300     f = diff(a)==0 & a(2:end) ~= 0;  
0301 
0302 
0303     % Just get the unique utc values.  The old way repeated values in the
0304     % following loop.  (MAS  27-Jan-2011)
0305     aUnique = unique(a(f));
0306 
0307     % For loop modified to handle previous two changes.  Should be much more
0308     % efficient.  (MAS  27-Jan-2011)
0309     for m=1:length(aUnique)
0310       ff = find(d.antenna0.thermal.utc == aUnique(m));
0311       indCryoBad(ff(2:length(ff))) = 1;
0312     end    
0313 
0314 
0315     % check for zero values in the time.
0316     indCryoBad = indCryoBad | d.antenna0.thermal.utc==0;
0317   
0318     
0319     % Now find the flagged regions.  If several flagged points are adjacent,
0320     % then this makes the following for loop much more efficient.
0321     % (MAS  27-Jan-2011)
0322     indCryoBadStart = find(diff(indCryoBad) > 0) + 1;
0323     indCryoBadEnd = find(diff(indCryoBad) < 0);
0324   
0325     if indCryoBad(1) == 1
0326       indCryoBadStart = cat(1,1,indCryoBadStart);      
0327     end
0328     if indCryoBad(end) == 1
0329       indCryoBadEnd = cat(1,indCryoBadEnd,utcLength);      
0330     end
0331   
0332 
0333     % fix the times in the cryo and keep track of bad values.
0334     timeDiff = median(diff(d.antenna0.thermal.utc(~indCryoBad)));
0335 
0336     % New for loop.  Designed to be much more efficient in the case where
0337     % several adjacent times are flagged.  (MAS  27-Jan-2011)
0338     for m=1:length(indCryoBadStart)
0339       interpLength = indCryoBadEnd(m) - indCryoBadStart(m) + 1;
0340       
0341       if indCryoBadStart(m) == 1 && indCryoBadEnd(m) ~= utcLength
0342         d.antenna0.thermal.utc(indCryoBadStart(m):indCryoBadEnd(m)) = ...
0343             d.antenna0.thermal.utc(indCryoBadEnd(m)+1) ...
0344             - timeDiff * (interpLength:-1:1)';
0345       elseif indCryoBadStart(m) == 1 && indCryoBadEnd(m) == utcLength
0346         display('interpRegisters:: ALL OF YOUR CRYO DATA ARE BAD');
0347         display('interpRegisters:: DO NOT RELY ON THESE DATA FOR MAPMAKING');
0348         d.antenna0.thermal.utc = nan(utcLength, 1);
0349       else
0350         d.antenna0.thermal.utc(indCryoBadStart(m):indCryoBadEnd(m)) = ...
0351             d.antenna0.thermal.utc(indCryoBadStart(m)-1) ...
0352             + timeDiff * (1:interpLength)';
0353       end
0354     end
0355 
0356   end
0357 
0358   % Same as before
0359   ff = find(indCryoBad);
0360 %  disp(max(ff));
0361   if length(ff) ~= utcLength
0362   
0363     utc    = d.antenna0.thermal.utc;
0364     cryo   = d.antenna0.thermal.ccTemperatureLoad;
0365 
0366 %     disp(size(cryo));
0367 %    disp(size(utc));
0368 %     disp(size(repmat(cryo,[1 round(length(utc)/length(cryo))])'));
0369 %
0370 %     if length(cryo) ~= length(utc)
0371 %         cryo = reshape(repmat(cryo,[1 round(length(utc)/length(cryo))])',length(utc),1);
0372 %     end
0373     
0374     utc(ff)    = [];
0375     cryo(ff)   = [];
0376 
0377     d.antenna0.thermal.ccTemperatureLoad = interp1(utc, cryo, ...
0378     d.antenna0.thermal.utc);
0379   else
0380     d.antenna0.thermal.ccTemperatureLoad = nan(utcLength,1);      
0381   end
0382   
0383   display('interpRegisters:: Number of bad cryocon times in track (due gpib communication dropout):');
0384   display(sprintf('interpRegisters:: %i samples out of %i samples', length(find(indCryoBad)), length(indCryoBad)));
0385 end
0386 
0387 % now we have to interpolate the time of the servo utc onto the time of the
0388 % thermal utc so we can match them up.
0389 if(hasCryo && hasServo)
0390   % match the missing servo times to the thermal sample, and vice-versa.
0391   % this is most easily done by checking the difference between the servo
0392   % and thermal utc (they should be roughly 1.05 seconds).
0393 
0394 %   disp(size(d.antenna0.servo.utc));
0395 %   disp(size(d.antenna0.thermal.utc));
0396   if(d.array.frame.utc(1)<55244)
0397     % offset should be 0.5seconds higher
0398     timeDiff = (d.antenna0.servo.utc - d.antenna0.thermal.utc)*24*60*60 - ...
0399     0.45;
0400   elseif(d.array.frame.utc(1)>=55244 & d.array.frame.utc(1) < 56138)
0401     timeDiff = (d.antenna0.servo.utc - d.antenna0.thermal.utc)*24*60*60 + ...
0402     0.05;
0403   else  % July 30, 2012:  new time offset in Control system
0404     timeDiff = (d.antenna0.servo.utc - d.antenna0.thermal.utc)*24*60*60 + 0.015;
0405   end
0406   
0407   indBadBoth = abs(timeDiff)>=0.01 | indServoBad | indCryoBad;
0408   
0409   indServoBad = indBadBoth;
0410   indCryoBad  = indBadBoth;
0411 end
0412   
0413 % we do not actually want to cut the data, we need to maintain it.  we just
0414 % want to flag it.
0415 dcut = d;
0416 % interpolate the servo and cryo onto the same timeframe as the receiver data
0417 
0418 if(hasServo)
0419   ff = find(indServoBad);
0420   % check if all data is bad
0421   if(length(ff) + 2 >= length(dcut.antenna0.servo.utc))
0422     display('interpRegisters:: ALL OF YOUR SERVO DATA ARE BAD');
0423     display('interpRegisters:: DO NOT RELY ON THESE DATA FOR MAPMAKING');
0424     d.antenna0.servo.az = nan(length(d.antenna0.receiver.utc), 1);
0425     d.antenna0.servo.el = nan(length(d.antenna0.receiver.utc), 1);
0426   else
0427     dcut.antenna0.servo.utc(ff) = [];
0428     dcut.antenna0.servo.fast_az_pos(ff) = [];
0429     dcut.antenna0.servo.fast_el_pos(ff) = [];
0430     dcut.antenna0.servo.fast_az_err(ff) = [];
0431     dcut.antenna0.servo.fast_el_err(ff) = [];  
0432     iAz = interp1(dcut.antenna0.servo.utc, dcut.antenna0.servo.fast_az_pos, ...
0433     dcut.antenna0.receiver.utc);
0434     iEl = interp1(dcut.antenna0.servo.utc, dcut.antenna0.servo.fast_el_pos, ...
0435     dcut.antenna0.receiver.utc);
0436     d.antenna0.servo.az = iAz;
0437     d.antenna0.servo.el = iEl;    
0438   end
0439 end
0440 if(hasCryo)
0441   ff = find(indCryoBad);
0442   if(length(ff) + 2 >= length(dcut.antenna0.thermal.utc))
0443     display('interpRegisters:: ALL OF YOUR CRYO DATA ARE BAD');
0444     display('interpRegisters:: DO NOT RELY ON THESE DATA FOR MAPMAKING');
0445     d.antenna0.thermal.coldLoad = nan(length(d.antenna0.receiver.utc), 1);
0446   else
0447     dcut.antenna0.thermal.utc(ff) = [];
0448     dcut.antenna0.thermal.ccTemperatureLoad(ff) = [];
0449     
0450     iLoad = interp1(dcut.antenna0.thermal.utc, ...
0451     dcut.antenna0.thermal.ccTemperatureLoad,dcut.antenna0.receiver.utc, ...
0452     'spline');
0453     d.antenna0.thermal.coldLoad = iLoad;
0454   end
0455 end
0456 
0457 % now we interpolate the servo and cryo con onto the same timescale.
0458 
0459 dcut = d;
0460 if(hasServo && hasCryo)
0461   ff = find(indBadBoth);
0462   utc = dcut.antenna0.thermal.utc;
0463   utc(ff) = [];
0464   load = dcut.antenna0.thermal.ccTemperatureLoad;
0465   load(ff) = [];
0466   
0467   % If all data are bad, if-then-else prevents crash.  (MAS 27-Jan-2011).
0468   if length(ff) ~= length(indBadBoth)  
0469     % interpolate the cryo onto the servo time.
0470     iLoad2 = interp1(utc, load, dcut.antenna0.servo.utc);
0471     d.antenna0.thermal.utc = d.antenna0.servo.utc;
0472     d.antenna0.thermal.ccTemperatureLoad = iLoad2;
0473   else
0474     iLoad2 = nan(length(dcut.antenna0.servo.utc),1);
0475     d.antenna0.thermal.utc = d.antenna0.servo.utc;
0476     d.antenna0.thermal.ccTemperatureLoad = iLoad2;
0477               
0478   end
0479 end
0480 
0481 if(hasServo)
0482   d.flags.interpFlags.medium = indServoBad;
0483   if(hasCryo)
0484     d.flags.interpFlags.medium = indServoBad;
0485   end
0486 elseif(hasCryo)
0487   d.flags.interpFlags.medium = indCryoBad;
0488 end
0489 
0490 if(hasCryo)
0491   d.flags.interpFlags.fast   = indRxBad;
0492 end
0493 
0494 
0495 return;