opacity

 Function opacity.m
 Inputs d = opacity(d)
 d the CBASS data struct
 outputs the data to d.antenna0.receiver.dataTcorr

0001 % Function opacity.m
0002 % Inputs d = opacity(d)
0003 % d the CBASS data struct
0004 % outputs the data to d.antenna0.receiver.dataTcorr
0005 
0006 function d = opacity(d)
0007 
0008 
0009 skydip = d.index.skydip.fast & ~d.index.source.fast & ~d.index.noise.fast;
0010 skydipStart=find((diff(d.index.skydip.fast)==1)); %get the skydip start positions
0011 skydipEnd=find((diff(d.index.skydip.fast)==-1));  %get the skydip end positions
0012 onSourceStart=find((diff(d.index.source.fast)==1)); % get the on source data
0013 onSourceEnd=find((diff(d.index.source.fast)==1)); 
0014 try
0015 d.antenna0.receiver.dataTcorr = d.antenna0.receiver.dataT; %copy the raw data to the output position
0016 catch
0017 disp('opacity:: No dataT')
0018 d.antenna0.receiver.dataTcorr = d.antenna0.receiver.data; %copy the raw data to the output position
0019 end
0020 
0021 
0022 
0023 numskydips = min(length(skydipStart),length(skydipEnd));
0024 for i=1:numskydips
0025     index=[skydipStart(i):skydipEnd(i)]; %get the indexing for the first skydip
0026     maxEl = find(d.antenna0.servo.el(index)==max(d.antenna0.servo.el(index))); %find the maximum elevation in the skydip
0027     minEl = find(d.antenna0.servo.el(index)==min(d.antenna0.servo.el(index))); %find the minimum elevation in the skydip
0028     start_scan = min(index(maxEl),index(minEl)); %not sure whether skydip is a dip or a rise so do this to check
0029     end_scan = max(index(maxEl),index(minEl));
0030     scan=[start_scan:end_scan]; %this skydip indexing
0031     %normdat=data./normVal;
0032     cscEl=csc(deg2rad(d.antenna0.servo.el)); %get the cosecant angles
0033     px1=polyfit(cscEl(scan),d.antenna0.receiver.dataTcorr(scan,1),1); %fit a straight line to the skydip vs cosecant
0034     px2=polyfit(cscEl(scan),d.antenna0.receiver.dataTcorr(scan,6),1);
0035     
0036     zero1=px1(1).*1 + px1(2); %get the estimated minimum atmospheric contribution (i.e at el=90 and cosecant(90)=1
0037     zero2=px2(1).*1 + px2(2);
0038     opacity1 = px1(1); %get the opacity
0039     opacity2 = px2(1);
0040     opacity{i} = [skydipStart(i) skydipEnd(i) opacity1 zero1 opacity2 zero2]; %store the various values to the structures
0041     
0042 end
0043 
0044 for i=1:length(opacity)
0045     opstruct.x(i)=mean(opacity{i}(2),opacity{i}(1)); %get a estimate of the index value when the opacity was obtained
0046     opstruct.op1(i) = opacity{i}(3); %opacity for channel 1 at this index value
0047     opstruct.op2(i) = opacity{i}(5); %opacity for channel 2 at this index value
0048     opstruct.zero1(i) = opacity{i}(4); %zero atmosphere value
0049     opstruct.zero2(i) = opacity{i}(6);
0050 end
0051 
0052 ind1 = find(~isnan(opstruct.op1));
0053 ind2 = find(~isnan(opstruct.op2));
0054 ind3 = find(~isnan(opstruct.zero1));
0055 ind4 = find(~isnan(opstruct.zero2));
0056 
0057 try
0058     opstruct.popacity1=polyfit(opstruct.x(ind1),opstruct.op1(ind1),1); %fit straight line to estimate the opacity change through the night
0059     opstruct.popacity2=polyfit(opstruct.x(ind2),opstruct.op2(ind2),1); %same again
0060     opstruct.pzero1=polyfit(opstruct.x(ind3),opstruct.zero1(ind3),1); %same again
0061     opstruct.pzero2=polyfit(opstruct.x(ind4),opstruct.zero2(ind4),1); %same again
0062 
0063 
0064 opacity1I = polyval(opstruct.popacity1,[1:length(d.antenna0.receiver.dataTcorr)])'; %get the opacity values interpolated across the night
0065 opacity2I = polyval(opstruct.popacity2,[1:length(d.antenna0.receiver.dataTcorr)])'; %get the opacity values interpolated across the night
0066 zero1I = polyval(opstruct.pzero1,[1:length(d.antenna0.receiver.dataTcorr)])';%get the zero values interpolated across the night
0067 zero2I = polyval(opstruct.pzero2,[1:length(d.antenna0.receiver.dataTcorr)])';%get the zero values interpolated across the night
0068 catch
0069 opacity1I = opstruct.op1(1);
0070 opacity2I = opstruct.op2(1);
0071 zero1I = opstruct.zero1(1);
0072 zero2I=opstruct.zero2(1);
0073 end
0074 
0075 d.antenna0.receiver.dataTcorr(:,1) = d.antenna0.receiver.dataTcorr(:,1) - zero1I - opacity1I.*cscEl; %correct for opacity
0076 d.antenna0.receiver.dataTcorr(:,6) = d.antenna0.receiver.dataTcorr(:,6) - zero2I - opacity2I.*cscEl;%correct for opacity
0077 p1 = polyfit([1:length(d.antenna0.receiver.dataTcorr)]',d.antenna0.receiver.dataTcorr(:,1),1);
0078 p2 = polyfit([1:length(d.antenna0.receiver.dataTcorr)]',d.antenna0.receiver.dataTcorr(:,6),1);
0079 %d.antenna0.receiver.dataTcorr(:,1) = d.antenna0.receiver.dataTcorr(:,1) - median(d.antenna0.receiver.dataTcorr(:,1));
0080 %d.antenna0.receiver.dataTcorr(:,6) = d.antenna0.receiver.dataTcorr(:,6) - median(d.antenna0.receiver.dataTcorr(:,6));
0081 %if ~(isnan(p1(1)) && isnan(p1(2)) && isnan(p2(1)) && isnan(p2(2))    )
0082     d.antenna0.receiver.dataTcorr(:,1) = d.antenna0.receiver.dataTcorr(:,1)-polyval(p1,[1:length(d.antenna0.receiver.dataTcorr)]')+5; %subtract a baseline
0083     d.antenna0.receiver.dataTcorr(:,6) = d.antenna0.receiver.dataTcorr(:,6)-polyval(p1,[1:length(d.antenna0.receiver.dataTcorr)]')+5;
0084 %else
0085 %    disp('opacity:: Cant fit a baseline in opacity.m- not removing one')
0086 %end
0087 
0088 
0089 end