cutObs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
   function [d indRx] = cutObs(d, type, logical)

   cuts the corresponding to a given observation type, given by the
   logical.

   valid types:  'source','elscan', 'calibrator', 'abscal', 'blank', 'skydip',
   'radio_point_cross', 'radio_point_scan', 'opt_point', 'beammap', 
    'noise_event', 'noise','ncp'

   valid logical: 'only', 'and', 'or', 'not'

   for example:
   for data that is labelled 'skydip', regardless of what other features
   are set:
   d = cutObs(d, 'skydip', 'or');
   if you want data that is only labelled as 'skydip', call
   as:  (keep in mind that some data is taken with many different labels.
   d = cutObs(d, 'skydip', 'only');

0001 function [d indRx] = cutObs(d, type, logical)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %   function [d indRx] = cutObs(d, type, logical)
0006 %
0007 %   cuts the corresponding to a given observation type, given by the
0008 %   logical.
0009 %
0010 %   valid types:  'source','elscan', 'calibrator', 'abscal', 'blank', 'skydip',
0011 %   'radio_point_cross', 'radio_point_scan', 'opt_point', 'beammap',
0012 %    'noise_event', 'noise','ncp'
0013 %
0014 %   valid logical: 'only', 'and', 'or', 'not'
0015 %
0016 %   for example:
0017 %   for data that is labelled 'skydip', regardless of what other features
0018 %   are set:
0019 %   d = cutObs(d, 'skydip', 'or');
0020 %   if you want data that is only labelled as 'skydip', call
0021 %   as:  (keep in mind that some data is taken with many different labels.
0022 %   d = cutObs(d, 'skydip', 'only');
0023 
0024 %   for data on 'calibrator' or 'noise', cutObs(d, {'calibrator', 'noise'},
0025 %   'or');
0026 %   for noise observations on a source only (not on calibrator):
0027 %     cutObs(d, {'noise', 'source'}, 'and');
0028 %
0029 %  for all observations that are not the noise source:
0030 %     cutObs(d, 'noise', 'not');
0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 
0033 %check the inputs
0034 if(nargin<2)
0035   error('Wrong number of inputs');
0036   help cutObs
0037 elseif(nargin==2)
0038   logical = 'or';
0039 end
0040 
0041 % change type to a cell array if only one input
0042 if(ischar(type))
0043   newType{1} = type;
0044   type = newType;
0045 end
0046 
0047 % check that all types are valid.
0048 featureFields;
0049 validTypes = fieldNames;
0050 
0051 typeValid = zeros(1,length(type));
0052 for m=1:length(type)
0053   if(~isempty(strmatch(type(m), validTypes, 'exact')))
0054     typeValid(m) = 1;
0055   end
0056 end
0057 if(~all(typeValid))
0058   error('Invalid type');
0059 end
0060 
0061 % 'only' option only works for a single type input
0062 if(strmatch(logical, 'only', 'exact'))
0063   if(length(type)>1)
0064     error('Current function supports the ''only'' option for single type');
0065   end
0066 end
0067 if(strmatch(logical, 'not', 'exact'))
0068   if(length(type)>1)
0069     error('Current function supports the ''only'' option for single type');
0070   end
0071 end
0072 
0073 
0074 
0075 % first determine the arrays on which to cut.
0076 % determine the corresponding inds for each data size
0077 [indReg, indServ, indRx] = determineInd(d, type, logical);
0078 
0079 % remove fields that we can't cut
0080 [d buf] = separate_fields(d);
0081 
0082 % cut them
0083 d = framecutSub(d, indRx, indServ, indReg);
0084 
0085 % add ones we couldn't cut back in
0086 d = mergestruct(d, buf);
0087 
0088 return;
0089 
0090 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0091 % old framecut program, works when things are properly sorted.
0092 
0093 function d = framecutSub(d, indRx, indServ, indReg)
0094 
0095 names=fieldnames(d);
0096 for i=1:length(names)
0097   if(eval(sprintf('isstruct(d.%s)',names{i})))
0098     eval(sprintf('d.%s=framecutSub(d.%s,indRx, indServ, indReg);',names{i},names{i}));
0099   else
0100     eval(sprintf('thisSize = size(d.%s,1);', names{i}));
0101     if(thisSize == length(indRx))
0102       eval(sprintf('d.%s=regcut(d.%s,indRx);',names{i},names{i}));      
0103     elseif(thisSize == length(indReg))
0104       eval(sprintf('d.%s=regcut(d.%s,indReg);',names{i},names{i}));            
0105     elseif(thisSize == length(indServ))
0106       eval(sprintf('d.%s=regcut(d.%s,indServ);',names{i},names{i}));            
0107     else
0108       error('Do not know which dimensino to cut over');
0109     end
0110   end
0111 end
0112 
0113 return
0114 
0115 %%%%%%%%%%%%%%%%%%%%%%%%
0116 function r=regcut(r,ind)
0117 
0118 switch(ndims(r))
0119   case 1
0120     r=r(ind);
0121   case 2
0122     r=r(ind,:);
0123   case 3
0124     r=r(ind,:,:);
0125   case 4
0126     r=r(ind,:,:,:);
0127 end
0128 
0129 return
0130 
0131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0132 function field = whichField(d, ind)
0133 
0134 lind  = length(ind);
0135 lserv = length(d.antenna0.servo.utc);
0136 lrec  = length(d.antenna0.receiver.utc);
0137 lnorm = length(d.array.frame.utc);
0138 
0139 if(lind == lserv)
0140   field = 'servo';
0141 elseif(lind == lrec)
0142   field = 'receiver';
0143 elseif(lind == lnorm)
0144   field = 'regular';
0145 else
0146   error('Vector over which to cut does not match dimensions in data');
0147 end
0148 
0149 return;
0150 
0151 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0152 function [ind_slow, ind_medium, ind_fast] = determineInd(d, type, logical)
0153 
0154 timeTypes = {'slow', 'medium', 'fast'};
0155 
0156 featureFields;
0157 validTypes = fieldNames;
0158 
0159 % first we determine the
0160 switch(logical)
0161   case 'or'
0162     % easiest case, we just combine everything
0163     ind_slow = zeros(size(d.index.source.slow));
0164     ind_medium = zeros(size(d.index.source.medium));
0165     ind_fast = zeros(size(d.index.source.fast));
0166 
0167     for m=1:length(type)
0168       for k=1:length(timeTypes)
0169     eval(sprintf('thisInd_%s = d.index.%s.%s;', timeTypes{k}, type{m}, ...
0170         timeTypes{k}));
0171     eval(sprintf('ind_%s = ind_%s | thisInd_%s;', timeTypes{k}, ...
0172         timeTypes{k}, timeTypes{k}));
0173       end
0174     end
0175 
0176   case 'and'
0177     % second easiest case, just do the same as above, but with ands instead
0178     % of ors
0179     ind_slow = ones(size(d.index.source.slow));
0180     ind_medium = ones(size(d.index.source.medium));
0181     ind_fast = ones(size(d.index.source.fast));
0182     for m=1:length(type)
0183       for k=1:length(timeTypes)
0184     eval(sprintf('thisInd_%s = d.index.%s.%s;', timeTypes{k}, type{m}, ...
0185         timeTypes{k}));
0186     eval(sprintf('ind_%s = ind_%s & thisInd_%s;', timeTypes{k}, ...
0187         timeTypes{k}, timeTypes{k}));
0188       end
0189     end    
0190 
0191   case 'not'
0192     % easiest case, we just combine everything
0193     ind_slow = zeros(size(d.index.source.slow));
0194     ind_medium = zeros(size(d.index.source.medium));
0195     ind_fast = zeros(size(d.index.source.fast));
0196 
0197     for m=1:length(type)
0198       for k=1:length(timeTypes)
0199     eval(sprintf('thisInd_%s = d.index.%s.%s;', timeTypes{k}, type{m}, ...
0200         timeTypes{k}));
0201     eval(sprintf('ind_%s = ind_%s | thisInd_%s;', timeTypes{k}, ...
0202         timeTypes{k}, timeTypes{k}));
0203       end
0204     end    
0205     
0206     ind_slow = ~ind_slow;
0207     ind_medium = ~ind_medium;
0208     ind_fast = ~ind_fast;
0209     
0210     
0211   case 'only'
0212     % in this case, you want the one where you have only the source you want
0213     % and nothing else
0214     ind_slow = zeros(size(d.index.source.slow));
0215     ind_medium = zeros(size(d.index.source.medium));
0216     ind_fast = zeros(size(d.index.source.fast));
0217     for m=1:length(validTypes)
0218       if(strmatch(validTypes{m}, type{1}, 'exact'))
0219     % this is the one we want to keep
0220     for k=1:length(timeTypes)
0221       eval(sprintf('thisInd_good_%s = d.index.%s.%s;', timeTypes{k}, validTypes{m}, ...
0222           timeTypes{k}));
0223     end
0224       else
0225     for k=1:length(timeTypes)
0226       eval(sprintf('thisInd_%s = d.index.%s.%s;', timeTypes{k}, validTypes{m}, ...
0227           timeTypes{k}));
0228       eval(sprintf('ind_%s = ind_%s | thisInd_%s;', timeTypes{k}, ...
0229           timeTypes{k}, timeTypes{k}));    
0230     end    
0231       end
0232     end
0233     ind_slow = thisInd_good_slow & ~ind_slow;
0234     ind_medium = thisInd_good_medium & ~ind_medium;
0235     ind_fast = thisInd_good_fast & ~ind_fast;
0236 end
0237 
0238 return;
0239  
0240      
0241 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0242 function [indReg, indServ, indRx] = determineIndOld(d, ind, field)
0243 
0244 halfSec  = 0.5/60/60/24;
0245 tenthSec = 0.1/60/60/24; 
0246 
0247 switch(field)
0248   case 'regular'
0249      indReg = ind;
0250      indRx = zeros(size(d.antenna0.receiver.utc));
0251      indServ = zeros(size(d.antenna0.servo.utc));
0252 
0253      times = d.array.frame.utc(ind);
0254      for m=1:length(times)
0255      f = find(d.antenna0.receiver.utc>=(times(m)-halfSec) & d.antenna0.receiver.utc<=(times(m)+halfSec));
0256          indRx(f) = 1;
0257      f = find(d.antenna0.servo.utc>=times(m)-halfSec & d.antenna0.servo.utc<=times(m)+halfSec);
0258          indServ(f) = 1;
0259      end
0260 
0261   case 'receiver'
0262      indRx = ind;
0263      indServ = zeros(size(d.antenna0.servo.utc));
0264      indReg  = zeros(size(d.array.frame.utc));
0265 
0266      times = d.antenna0.receiver.utc(ind);
0267      % all we have to do is see if the times are within 0.5s for regular and tenthSec for servo.
0268      for m=1:length(times)
0269      f = find(d.array.frame.utc>=times(m)-halfSec & d.array.frame.utc<=times(m)+halfSec);
0270          indReg(f) = 1;
0271      f = find(d.antenna0.servo.utc>=times(m)-tenthSec & d.antenna0.servo.utc<=times(m)+tenthSec);
0272          indServ(f) = 1;
0273      end
0274 
0275   case 'servo'
0276      indServ = ind;
0277      indRx   = zeros(size(d.antenna0.receiver.utc));
0278      indReg  = zeros(size(d.array.frame.utc));
0279 
0280      times = d.antenna0.servo.utc(ind);
0281      % all we have to do is see if the times are within 0.5s for regular and 0.1 for receiver
0282      for m=1:length(times)
0283      f = find(d.array.frame.utc>=times(m)-halfSec & d.array.frame.utc<=times(m)+halfSec);
0284          indReg(f) = 1;
0285      f = find(d.antenna0.servo.utc>=times(m)-tenthSec & d.antenna0.servo.utc<=times(m)+tenthSec);
0286          indRx(f) = 1;
0287      end
0288 end
0289 
0290 indRx   = logical(indRx);
0291 indServ = logical(indServ);
0292 indReg  = logical(indReg);
0293 
0294 return;