Description of calibratorScan

0001 function [dataSaved sourceName dates] =  calibratorScan(startDate,endDate,fileLocation)
0002 %This function produces calibrator scan plots (defined in the obs_log.html)
0003 %and saves them
0004 %Input Requirements
0005 %startDate : StartDate in C-BASS format
0006 %endDate : End Date in C-BASS format
0007 %fileLocation : Location to store the files
0008 % example usage
0009 % calibratorScan('25-Oct-2010:03:45:10','10-Nov-2010:04:10:45',[home,'/angela/calibrationScans/'])
0010 % file output in this format
0011 % Off Source Noise Diode contrib
0012 % On Source ND contrib
0013 % sky signal off source
0014 % Sky signal+ND off source
0015 % Sky+source
0016 %Sky+source+ND
0017 %compression ration off Source
0018 %compression raio on Source
0019 % (3columns) az peak, azhpbw, azoffset
0020 % (3 columns) same for elevation
0021 %azimuth during scan
0022 %elevation during scan
0023 
0024 %CJC 11 Nov 2010
0025 %ACT 16 Nov 2010 - added outputs of source name, dates & temperatures
0026 %clear
0027 
0028 
0029 %%%YOU SHOULD AMEND THE FOLLOWING TO YOUR LOCAL INSTALLATION
0030 
0031 [home,installeddir] = where_am_i();
0032 logfile = ([home, '/',installeddir,'/log/obs_log.html']);
0033 pipeline_file = @pipelinedData %Pick your pipeline file to use
0034 %%%%%%%%%%%%%%%%%%%%
0035 
0036 
0037 runningCounter=0;
0038 close all
0039 
0040 %t3=extract_obslog(logfile,'scan_calibrators_sjcm.sch',startDate,endDate);
0041 t3=extract_obslog(logfile,'scan_calibrators',startDate,endDate);
0042 %disp('calibratorScan:: HERE -1')
0043 t=t3;%[t1;t2;t3];
0044 numScans=size(t3);
0045 dataSaved=zeros(1,8)
0046 figure;
0047 for j=1:numScans(1)
0048 %    disp('calibratorScan:: HERE -2')
0049 %j=2;
0050 %i=1;
0051 
0052 %figure
0053 try
0054 d=pipe_read(t3{j,1},t3{j,2});
0055 
0056 
0057 % Interpolate the ND temperatures onto the receiver.utc timescale for later
0058 % use - ACT
0059 d.antenna0.thermal.fasttemps = interp1(d.antenna0.tracker.utc,d.antenna0.thermal.dlpTemperatureSensors, d.antenna0.receiver.utc);
0060 
0061 %
0062 
0063 % %% Do the alpha and r-factor correction
0064 % % 1. Calculate the alpha corrections
0065 %     disp('calibratorScan:: Applying the alpha and r-factor corrections')
0066 %
0067 %     try
0068 %     [t,A,G,T,horiz,equa, offStartPos, onEndPos, offEndPos, onStartPos] = calculateAlpha(d);
0069 %     catch
0070 %         disp('calibratorScan:: (Failed to calculate Alphas')
0071 %     end
0072 %
0073 % % 2. Applies the alpha corrections to the data and overwrites the switchData register
0074 %     try
0075 %     d = applyAlpha2(d,t,A,G);
0076 %     catch
0077 %         disp('calibratorScan:: Failed to apply the alpha corrections')
0078 %     end
0079 %
0080 % % 3. Some pre-processing before before r-factor correction:
0081 %     try
0082 %     thisAlpha = [t A G T horiz equa];
0083 %     d.correction.alpha.indices = [offStartPos' onEndPos' offEndPos' onStartPos'];
0084 %     d.correction.alpha.values = thisAlpha;
0085 %     catch
0086 %         disp('calibratorScan:: Failed to pre-process data ready for r-factors')
0087 %     end
0088 %
0089 % % 4. Calculate the r-factor during the noise diode off events (this is why the above steps need to have been done)
0090 %     try
0091 %     r = calculateRfactor(d);
0092 %     meanr = mean(r)
0093 %     catch
0094 %         disp('calibratorScan:: Failted to calculate r-factors')
0095 %     end
0096 %
0097 % % 5. Apply the r-factor, by subtracting a smoothed load signal (scaled by the r factor) from the sky signal.
0098 % %    smoothLength is the length of the smoothing window function in seconds. This overwrites the data register.
0099 %     try
0100 %     smoothLength = 1
0101 %     rval =  [1.9,1.8,1.8,1.7];
0102 %     d = calculateStokes2_samer(d,r,smoothLength,rval);
0103 %     catch
0104 %         disp('calibratorScan:: Failed to apply r-factors')
0105 %     end
0106 %%
0107 
0108 
0109 
0110 date = t3(j,1)
0111 indnoise=find(d.index.noise.fast); %get the noise diode indices
0112 indskydip=find(d.index.skydip.fast); %get the skydip data
0113 indcalib = find(d.index.calibrator.fast); %get the calibrator scan features
0114 offSourceFast=interp1(d.antenna0.tracker.utc,double(d.antenna0.tracker.offSource),d.antenna0.receiver.utc,'nearest'); %here we just interpolate from 5Hz to 100Hz
0115 indoffSource=find(offSourceFast); %and here we find the off Source data
0116 disp('calibratorScan:: HERE 0')
0117 sources = unique(d.antenna0.tracker.source); %find the uniques sources defined in this scan period
0118 for i=1:length(sources)
0119     
0120    disp('calibratorScan:: HERE 1')
0121     sourceInd = strcmp(sources(i),d.antenna0.tracker.source); %get data relating to one of the sources (might only be one!!)
0122     testData = framecut(d,sourceInd); %cut to only this source data
0123     feat=interp1(testData.array.frame.utc,double(testData.array.frame.features),testData.antenna0.receiver.utc,'nearest');
0124     noiseInd=find(testData.index.noise.fast==1); %get the noise diode events
0125     noiseIndDiff=diff(testData.index.noise.fast); %being very sneaky I find where the diode turns on and off
0126     changeNoiseEventOn=find(noiseIndDiff==1);%find the places where the noise diode turns on
0127     changeNoiseEventOff=find(noiseIndDiff==-1);%find the places where the noise diode turns off
0128     disp('calibratorScan:: HERE 2')
0129     numberNoiseEvents=size(changeNoiseEventOn); %now we need to know how many calibrator scans through this source occurred? Sometime more than one- I use the number of ND events to calculate- there are 3 noise diode events per calibrator scan in scan_calibrators_sjcm.sch
0130     try
0131     [az,el]=calcAzEl(testData); %calculate the source azimuth and elevation
0132     azErr = (testData.antenna0.servo.az - az)./sec(deg2rad(testData.antenna0.servo.el)) ; %use the source azimuth calculated before to calculate offset
0133     elErr = testData.antenna0.servo.el - el ; %use the source elevation calculated before to calculate offset
0134     catch
0135         disp('calibratorScan:: Failed to calculate the az/el offsets');
0136     end
0137     
0138     for l=1:numberNoiseEvents(1)/3
0139         
0140         runningCounter=runningCounter+1; %have a running counter of the number of scans processed- mainly for labelling purposes
0141         azimuthScan= find((abs(diff(testData.antenna0.servo.az))>0.2/100) & (abs(azErr(1:length(azErr)-1))<1.5));%get the az scans (very sneaky!)- basically find areas where we are moving >0.2deg/s in azimuth and are also nearby the source
0142         elevationScan= find(abs(diff(testData.antenna0.servo.el))>0.2/100 & (abs(elErr(1:length(elErr)-1))<1.5)); %get the el scans (ditto)
0143         
0144         
0145         k=3*(l-1)+1; %there are three ND events so need to index appropriately
0146         NoiseEvent{k}=[changeNoiseEventOn(k)+10:changeNoiseEventOff(k)-10]; %first noise diode on Source
0147         NoiseEvent{k+1}=[changeNoiseEventOn(k+1)+10:changeNoiseEventOff(k+1)-10]; %second Noise diode on Source
0148         NoiseEvent{k+2}=[changeNoiseEventOn(k+2)+10:changeNoiseEventOff(k+2)-10]; %third Noise diode off Source
0149     
0150         %here we need to exclude any azimuth scans outside of the area
0151         %defined by the first two noise diode events- see the plots
0152         temp=find(azimuthScan>min(NoiseEvent{k+1}));
0153         temp2 = find(azimuthScan<max(NoiseEvent{k}) );
0154         azimuthScan([temp2;temp])=[]; %only use the scans in between the two noise diode events
0155          
0156         %and do the same as above for the elevation
0157         temp=find(elevationScan>min(NoiseEvent{k+1}) );
0158         temp2 = find(elevationScan<max(NoiseEvent{k}) );
0159         elevationScan([temp2;temp])=[];%only use the scans in between the two noise diode events
0160         
0161         OffSourceNoNoiseStart{l} = [min(NoiseEvent{k+2})-400:min(NoiseEvent{k+2})-50];%get offSource index values before ND and off Source
0162         OffSourceNoNoiseEnd{l} = [max(NoiseEvent{k+2})+50:max(NoiseEvent{k+2})+400]; %get offSource index values after ND and offSource
0163         OnSourceNoNoiseStart{l} = [min(NoiseEvent{k+1})-400:min(NoiseEvent{k+1})-50];%get onSource index values before ND and onSource
0164         OnSourceNoNoiseEnd{l} = [max(NoiseEvent{k+1})+50:max(NoiseEvent{k+1})+400];%get onSource index values after ND and onSource
0165     
0166        % Angela Added a line to grab the temperature in the cone near the
0167        % noise diode when (sky+ND)
0168         NDTemp{l} = mean(mean(testData.antenna0.thermal.fasttemps(NoiseEvent{k+2},2:4)));
0169     AirTemp{l} = mean(testData.array.weather.airTemperature);
0170         %
0171         NoiseDiodeOffSource{l}=mean(testData.antenna0.receiver.data(NoiseEvent{k+2},1)); %get the Noise diode signal (sky+ND) we see when off the source
0172         NoiseDiodeOnSource{l}=mean(testData.antenna0.receiver.data(NoiseEvent{k+1},1)); %get the Noise diode signal (sky+source+ND) we see when on the source
0173         SkyOffSource{l} = mean(testData.antenna0.receiver.data([OffSourceNoNoiseStart{l} ],1)); %get the sky signal (sky) we see when off the source
0174         SkyOnSource{l} = mean(testData.antenna0.receiver.data([OnSourceNoNoiseStart{l}],1));%get the sky signal (sky+source) we see when off the source
0175     
0176         OutputIntensityValues{l}=[SkyOffSource{l};NoiseDiodeOffSource{l};SkyOnSource{l};NoiseDiodeOnSource{l}];
0177     
0178         NDContrib{l}(1)=OutputIntensityValues{l}(2)-OutputIntensityValues{l}(1); %offSource Noise diode signal
0179         NDContrib{l}(2)=OutputIntensityValues{l}(4)-OutputIntensityValues{l}(3); %On Source Noise Diode signal
0180         
0181     
0182         compressionRatio{l}(1) = (NDContrib{l}(1)/SkyOffSource{l}); %compression ratio off source
0183         compressionRatio{l}(2) =(NDContrib{l}(2)/SkyOnSource{l}); %compression ration on source
0184       %this is just so the plots don not drive me mad
0185             clf 
0186             h1=subplot(3,1,1)
0187             hold
0188             atemp=[min(NoiseEvent{k}):max(NoiseEvent{k+2})];
0189             plot(atemp,testData.antenna0.receiver.data(atemp,1));
0190   
0191        %and now we generate the plots-
0192         plot(NoiseEvent{k+2},testData.antenna0.receiver.data(NoiseEvent{k+2},1),'r');
0193         plot([OffSourceNoNoiseStart{l} OffSourceNoNoiseEnd{l}],testData.antenna0.receiver.data([OffSourceNoNoiseStart{l} OffSourceNoNoiseEnd{l}],1),'y');
0194         plot(NoiseEvent{k+1},testData.antenna0.receiver.data(NoiseEvent{k+1},1),'g');
0195         plot([OnSourceNoNoiseStart{l} OnSourceNoNoiseEnd{l}],testData.antenna0.receiver.data([OnSourceNoNoiseStart{l} OnSourceNoNoiseEnd{l}],1),'c');
0196         plot(azimuthScan,testData.antenna0.receiver.data(azimuthScan,1),'m.');
0197         plot((elevationScan),testData.antenna0.receiver.data(elevationScan,1),'.k');
0198         
0199         NDString{l}=['ND on = ',num2str(NDContrib{l}(2)),10,'ND off = ',num2str(NDContrib{l}(1)),10,'Sky = ',num2str(SkyOffSource{l}),10,'Sky+ND = ',num2str(NoiseDiodeOffSource{l}),10,'Sky + Source = ',num2str(SkyOnSource{l}),10,'Sky+Source+ND = ',num2str(NoiseDiodeOnSource{l}),10,'Diode Slope off Source %= ',num2str(100*compressionRatio{l}(1)),10,'Diode Slope on Source %= ',num2str(100*compressionRatio{l}(2))] %the 10 is for a newline!!
0200         text(max(NoiseEvent{k}),max(testData.antenna0.receiver.data(:,1)-1.5), NDString{l});
0201         
0202         %now to plot the cross scans
0203         subplot(3,1,2)
0204         x=azErr(azimuthScan);
0205         y=testData.antenna0.receiver.data(azimuthScan,1);
0206         plot(x,y,'.')
0207         %%here we get estimates for the azimuth gaussian fit
0208         A=max(y)-min(y);
0209         x0 = x(find(y==max(y)));
0210         HPBWInd=find((y>min(y)+0.47.*A) & (y<min(y)+0.53.*A));
0211         t1=x(HPBWInd);
0212         t2=t1<x0;
0213         lowSide=mean(t1(t2));
0214         t1=x(HPBWInd);
0215         t2=t1>x0;
0216         highSide=mean(t1(t2));
0217         width = highSide-lowSide;
0218         sigma = width./2.3548; %relation between sigma and HPBW
0219         m=0; %first guess!
0220         c= min(y); %first guess
0221         beta=[A x0 sigma m c];
0222         beta0=nlinfit(x,y,@gaussfit,beta);
0223         hold
0224         plot(x,gaussfit(beta0,x),'.r');
0225         HPBWend=x0+sigma*2.3548/2;
0226         HPBWstart=x0-sigma*2.3548/2;
0227         HPBWDeg=HPBWend-HPBWstart;
0228         StringForPlot=['A = ',num2str(A),10,'HPBW = ',num2str(HPBWDeg),10,'Offset = ',num2str(x0)];
0229         text(x0,min(y)+1, StringForPlot);
0230         aAZ=A;
0231         HPBWaz=HPBWDeg;
0232         offsetAz=x0;
0233         %xlim([-3 +3]);
0234         grid on
0235         set(gca,'XMinorGrid','on');
0236         
0237         
0238         title('Azimuth Cross Scans');
0239         xlabel('Degrees');
0240         ylabel('ADC units');
0241         xlim([-2 2]);
0242      
0243         %%here we get estimates for the elevation gaussian fit
0244         subplot(3,1,3)
0245         x=elErr(elevationScan);
0246         y=testData.antenna0.receiver.data(elevationScan,1);
0247         %plot(x,y,'.')
0248         A=max(y)-min(y);
0249         x0 = x(find(y==max(y)));
0250         HPBWInd=find((y>min(y)+0.47.*A) & (y<min(y)+0.53.*A));
0251         t1=x(HPBWInd);
0252         t2=t1<x0;
0253         lowSide=mean(t1(t2));
0254         t1=x(HPBWInd);
0255         t2=t1>x0;
0256         highSide=mean(t1(t2));
0257         width = highSide-lowSide;
0258         sigma = width./2.3548; %relation between sigma and HPBW
0259         m=0; %first guess!
0260         c= min(y); %first guess
0261         beta=[A x0 sigma m c];
0262         beta0=nlinfit(x,y,@gaussfit,beta);
0263         %hold
0264         
0265         hold
0266         plot(x,y,'.');
0267         plot(x,gaussfit(beta0,x),'.r');
0268         HPBWend=x0+sigma*2.3548/2;
0269         HPBWstart=x0-sigma*2.3548/2;
0270         HPBWDeg=HPBWend-HPBWstart;
0271         StringForPlot=['A = ',num2str(A),10,'HPBW = ',num2str(HPBWDeg),10,'Offset = ',num2str(x0)];
0272         aEL=A;
0273         HPBWel=HPBWDeg;
0274         offsetEl=x0;
0275         text(x0,min(y)+1, StringForPlot);
0276         title('Elevation Cross Scans');
0277         xlabel('Degrees');
0278         ylabel('ADC units');
0279          xlim([-2 2])
0280    
0281          grid on
0282         set(gca,'XMinorGrid','on');
0283         
0284         %here we calculate the average of the azimuth and elevation of the
0285         %source during the cross scans
0286         azimuths=testData.antenna0.servo.az(azimuthScan);
0287         elevations = testData.antenna0.servo.el(elevationScan);
0288         azerrs=azErr(azimuthScan);
0289         elerrs=elErr(elevationScan);
0290         azimuthIndex=find((azerrs<(offsetAz+10/1000) & azerrs>(offsetAz-10/1000)));
0291         elevationIndex=find((elerrs<(offsetEl+10/1000) & elerrs>(offsetEl-10/1000)));
0292         azimuthDuringScans=mean(azimuths(azimuthIndex));
0293         elevationDuringScans=mean(elevations(elevationIndex));
0294         %%%%%%%%%
0295         
0296         
0297         titleText=['Calibrator Scan through Source ',char(sources(i)),10,mjd2string(testData.antenna0.receiver.utc(1)),' ',mjd2string(testData.antenna0.receiver.utc(length(testData.antenna0.receiver.utc)))]
0298         title(h1,titleText);
0299         xlabel(h1,'Sample Number [100Hz]');
0300         ylabel(h1,'Channel 1 ADC');
0301         legend(h1,'Data','Noise Event Off Source','Off Source Sky','Noise Event On Source','On Source Sky','Azimuth Cross Scans','Elevation Cross scans');
0302         fullscreen = get(0,'ScreenSize');
0303         set(gcf,'Position',[0 -50 fullscreen(3) fullscreen(4)])
0304    
0305         set(gcf,'PaperPositionMode','auto')
0306         %
0307         size(NDTemp)
0308         NDTemp{l}
0309         print('-dpng', [fileLocation,num2str(runningCounter),'CalibScan','--',char(sources(i)),num2str(j),'-',num2str(i),'-',num2str(l),'.png'],'-r300');
0310         dataSaved(runningCounter,:)=[NDContrib{l}(2) NDContrib{l}(1) ...
0311           SkyOffSource{l} NoiseDiodeOffSource{l} SkyOnSource{l} ...
0312           NoiseDiodeOnSource{l} 100*compressionRatio{l}(1) 100*compressionRatio{l}(2) aAZ HPBWaz offsetAz aEL HPBWel offsetEl azimuthDuringScans elevationDuringScans NDTemp{l} AirTemp{l}];
0313         sourceName{runningCounter}=sources(i);
0314         dates{runningCounter}=date;
0315     
0316     end
0317 end
0318 catch
0319     disp('calibratorScan:: Failed the entire batch process')
0320 end
0321 end
0322    figure 
0323    plot(dataSaved(:,5),dataSaved(:,1))
0324    hold
0325    plot(dataSaved(:,3),dataSaved(:,2),'r')
0326    title('Noise Diode Linearity')
0327    xlabel('Sky Temperature [ADC]')
0328    ylabel('Noise Diode contribution')
0329    legend('On Source ND','Off Source ND')
0330    dlmwrite([fileLocation,'CalibratorDate',startDate,'-',endDate,'.csv'],dataSaved)
0331    close all
0332    
0333     
0334     
0335
calibratorScan

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
