load_correct

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

 function [C] = correct_load(d)

   To apply a correction factor to remove the effect
   of fluctuations in the cold load
   Do this on each of the 6 channels
    act

  modified sjcm to fit with pipeline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [K_best,sky_res,load_res,d] = load_correct(d)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % function [C] = correct_load(d)
0006 %
0007 %   To apply a correction factor to remove the effect
0008 %   of fluctuations in the cold load
0009 %   Do this on each of the 6 channels
0010 %    act
0011 %
0012 %  modified sjcm to fit with pipeline
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 
0015 for chan=1:6
0016   % First let's select the flagged data points only
0017   temp = [d.antenna0.receiver.data(:,chan) d.antenna0.thermal.coldLoad];
0018   temp(d.flags.fast(:,chan)) = [];
0019 
0020   % Now we should remove the mean of the load temperature
0021   load_res = temp(:,2)-mean(temp(:,2));
0022   
0023   % Now let's detrend and remove the mean from the sky data
0024   sky_res = temp(:,1);
0025   
0026   % At this point have de-trended , mean removed data in sky_res and load_res
0027   % all flagged data masked (not in these arrays)
0028   
0029   % Now let's try just finding the factor by which we need to multiply the load to subtract it from the sky signal
0030   % C = sum(sky.load)
0031   % Minimise via factor K such that C_min = sum((sky-Kload).load)
0032   
0033   % First we need to shift load_res around in time to match the sky_res phase
0034   % - thermal lag??
0035   
0036   for i=0:100
0037     load_shift = circshift(load_res,i);
0038     index = i+101;
0039     value(index) = (dot((sky_res - load_shift), load_shift));
0040     ivalue(index) = i;
0041   end
0042   
0043   
0044   [H,I]=max(value);
0045   required_shift = ivalue(I);
0046   
0047   d.antenna0.thermal.coldLoad = circshift(d.antenna0.thermal.coldLoad,required_shift);
0048   load_res = circshift(load_res,required_shift);
0049   test_function = @(x) (std((sky_res - x*load_res)));
0050   
0051   [x,f_val] = fminsearch(test_function, 1);
0052   
0053   K_best = x;
0054   
0055   % Now subtract (K*load) from the sky signal and plot
0056   blah = K_best*(d.antenna0.thermal.coldLoad - mean(d.antenna0.thermal.coldLoad(~d.flags.fast(:,chan))));
0057   
0058   corrected  = d.antenna0.receiver.data(:,chan) - blah;
0059   original = d.antenna0.receiver.data(:,chan);
0060   d.antenna0.receiver.data(:,chan)=corrected;
0061   
0062   % Now look at power spectrum - temporarily set flagged datapoints to zero
0063   corrected(d.flags.fast(:,chan)) = 0;
0064   original(d.flags.fast(:,chan)) = 0;
0065   
0066   %Put thermal load data back to how it was originally
0067   d.antenna0.thermal.coldLoad = circshift(d.antenna0.thermal.coldLoad,-required_shift);
0068 end
0069 
0070 return;
0071