calc_loadcorrect_final

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 function [K_best,shift,dorig] = calc_loadcorrect_final(d)

   To calculate the correction factor to remove the effect
   of fluctuations in the cold load

   Do this on each of the 6 channels after alpha corrs and before r-factor corrections
   Returns arrays of best fit amplitude and shift (offset of thermal data) for each channel

   Use this function on a nice bit of static,clean data --> the factors should then be
   useable on any other data. 

   Should check this works over various timescales/ conditions before
   releasing as final

   act
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [K_best,shift,dorig] = calc_loadcorrect_final(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % function [K_best,shift,dorig] = calc_loadcorrect_final(d)
0006 %
0007 %   To calculate the correction factor to remove the effect
0008 %   of fluctuations in the cold load
0009 %
0010 %   Do this on each of the 6 channels after alpha corrs and before r-factor corrections
0011 %   Returns arrays of best fit amplitude and shift (offset of thermal data) for each channel
0012 %
0013 %   Use this function on a nice bit of static,clean data --> the factors should then be
0014 %   useable on any other data.
0015 %
0016 %   Should check this works over various timescales/ conditions before
0017 %   releasing as final
0018 %
0019 %   act
0020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 
0022 % This will work on ony the selected 'static' dataT (r-factored and
0023 % corrected)
0024 
0025 % Do this on blank, static data
0026 dcut = cutObs(d, 'blank', 'only');
0027 
0028 % change the flags so that they match the data
0029 flags = convertFlagToDataSize(dcut.flags, ...
0030     size(dcut.antenna0.receiver.data,2));
0031 
0032 numDims = size(dcut.antenna0.receiver.data,2);
0033 
0034 
0035 for chan=1:numDims
0036   
0037   if(isfield(dcut.antenna0.receiver, 'dataT'))
0038     isT = 1;
0039     temp = [dcut.antenna0.receiver.dataT(:,chan), dcut.antenna0.thermal.coldLoad(:)];
0040   else
0041     isT = 0;
0042     temp = [dcut.antenna0.receiver.data(:,chan),  dcut.antenna0.thermal.coldLoad(:)];
0043   end
0044 
0045   % apply flags
0046   %ind = logical(flags.fast(:,chan));
0047   %temp(ind,:) = nan;
0048   
0049   % Now we should remove the mean of the load temperature
0050   load_res = temp(:,2)-nanmean(temp(:,2));
0051   
0052   % Already r-factored so do not need to remove mean from the main data set
0053   sky_res = temp(:,1) - nanmean(temp(:,1));
0054 
0055   % we do not want the nan values to contribute, so set them to zero.
0056   load_res(isnan(load_res)) = 0;
0057   sky_res(isnan(sky_res)) = 0;
0058   
0059   % At this point have de-trended , mean removed data in sky_res and load_res
0060   % all flagged data masked (not in these arrays)
0061   
0062   % Now let us try just finding the factor by which we need to multiply the
0063   % load to subtract it from the sky signal
0064   % C = sum(sky.load)
0065   % Minimise via factor K such that C_min = sum((sky-Kload).load)
0066   
0067   % First we need to shift load_res around in time to match the sky_res phase
0068   % - thermal lag??
0069 
0070   % want over a full period that is not flagged
0071    
0072   % 17-Feb-2011 (MAS):  This has the potential to be really, really slow.
0073 %   for mm=1:length(ind)-100
0074 %    if(ind(mm)==1)
0075 %      consecutiveLength(mm) = 0;
0076 %    else
0077 %      ff = find(ind(mm:mm+100)==1,1);
0078 %      if(isempty(ff))
0079 %     consecutiveLength(mm) = 100;
0080 %      else
0081 %     consecutiveLength(mm) = ff;
0082 %      end
0083 %    end
0084 %  end
0085 %  fstart = find(consecutiveLength == max(consecutiveLength),1);
0086 %  fend   = fstart + max(consecutiveLength)-1;
0087 
0088 ind=0;
0089 %   % This is a lot faster.
0090   if(all(ind)==0)
0091     startInd = 1;
0092     stopInd  = length(ind);
0093   else
0094     diffInd = cat(1,0,diff(ind));
0095     startInd = find(diffInd < 0);
0096     stopInd = find(diffInd > 0);
0097   end
0098   
0099     
0100   length(startInd);
0101   length(stopInd);
0102 
0103   % check they are the same length.
0104   if(length(startInd)~=length(stopInd))
0105     if(length(startInd)<length(stopInd))
0106       if(startInd(1)>stopInd(1))
0107     stopInd(1) = [];
0108       end
0109     else
0110       if(startInd(2)<stopInd(1))
0111     startInd(1) = [];
0112       end
0113     end
0114   end
0115   
0116   if (length(startInd)>0)
0117   
0118     % First - if all the data is ok i.e no unflagged data
0119     
0120     if (length(stopInd)==0)
0121       fstart = startInd(1);
0122       fend = fstart + 99;
0123     else
0124       % This bit assumes that there is both flagged and unflagged data
0125       if startInd(1) > stopInd(1)
0126     %  keyboard
0127     startInd = cat(1,1,startInd);
0128     stopInd = cat(1,stopInd,length(ind));
0129     if max(stopInd - startInd) >= 100    % I have just added this in to prevent crashes
0130        firstLong = find((stopInd - startInd) >= 100,1); % when max(stopInd - startInd)
0131     else                                                % is less than 100
0132         firstLong = find((stopInd - startInd) >= 50,1); % MI 01/08/11
0133     end                                                  
0134     fstart = startInd(firstLong);
0135     fend = fstart + 99;
0136       else
0137 if startInd(end) > stopInd(end)
0138       stopInd = cat(1,startInd(1)-1,stopInd,length(ind));
0139       startInd = cat(1,1,startInd);
0140     if max(stopInd - startInd) >= 100
0141        firstLong = find((stopInd - startInd) >= 100,1);
0142     else 
0143         firstLong = find((stopInd - startInd) >= 50,1);
0144     end 
0145       fstart = startInd(firstLong);
0146       fend = fstart + 99;
0147 else
0148     if max(stopInd - startInd) >= 100
0149        firstLong = find((stopInd - startInd) >= 100,1);
0150     else 
0151         firstLong = find((stopInd - startInd) >= 50,1);
0152     end 
0153       fstart = startInd(firstLong);
0154       fend = fstart+99;
0155     end
0156       end
0157     end
0158     
0159     shiftIndex = 0;
0160      
0161    % for mm=fstart:fend
0162       load_shift = circshift(load_res, shiftIndex);
0163       value(shiftIndex+1) = (dot((sky_res - load_shift), load_shift));
0164       shiftIndex = shiftIndex+1;
0165    % end
0166     
0167   [H,I]=max(value);
0168   
0169   required_shift = I-1;
0170   shift(chan) = (required_shift);
0171   
0172   load_res = circshift(load_res,required_shift);
0173   
0174   test_function = @(x) (std((sky_res - x*load_res)));
0175   
0176   [x,f_val] = fminsearch(test_function, 1);
0177   
0178   K_best(chan) = x;
0179   
0180   
0181   %Put thermal load data back to how it was originally
0182   
0183   load_res = circshift(load_res,-required_shift);
0184   
0185   else
0186   % And what if all the data is flagged
0187 
0188     if (chan ==1)
0189          disp('calc_loadcorrect_final:: All your blank data is flagged') % put warning in once on first channel
0190          disp('calc_loadcorrect_final:: Will not apply a correction')
0191     end
0192   K_best(chan) = 0;   %% Once we know what the stability of this
0193   shift(chan)  = 0;   %% correction is we can use logged values
0194   end
0195   end
0196 
0197  
0198   
0199 % Apply correction to dcut data and look at power spectrum to check it
0200 % worked
0201 for chan = 1:numDims
0202     shifted_load = circshift(load_res,shift(chan));
0203     correction = K_best(chan)*(shifted_load - mean(shifted_load));
0204 
0205     if(isT)
0206       dcut.antenna0.receiver.dataT(:,chan) = ...
0207       (dcut.antenna0.receiver.dataT(:,chan))- correction;
0208     else
0209       dcut.antenna0.receiver.data(:,chan) = ...
0210       (dcut.antenna0.receiver.data(:,chan))- correction;
0211     end
0212 end
0213 
0214 
0215 % For checking - keep an array which is the original data for comparison
0216 
0217 dorig = d;
0218 
0219 
0220 return;
0221