calibrate_IQUV2

 By default, save the calibration matrix

0001 function [calmat data P fval hessian]=calibrate_IQUV2(cal_source,save_path,start_date,source_name,IntensityType,save_data,new_bad_rasters)
0002 
0003 % By default, save the calibration matrix
0004 if(~exist('save_data'))
0005        save_data = 1
0006 end
0007 if(~exist('new_bad_rasters'))
0008        new_bad_rasters = []
0009 end
0010 %save_path ='/test_fits_June2013/';
0011 
0012 %load('./test_fits/07-Nov-2012:07:22:02_current_PA_IQUV.mat')
0013 [home,installeddir] = where_am_i();
0014 maindir= [home,'/',installeddir,'/',save_path,'/'];
0015 
0016 load([maindir,start_date,'_',source_name,'_',IntensityType,'_PA_IQUV.mat'])
0017 % Need to form amplitudes of I,Q,U,and corresponding PA
0018 
0019 %for(i=1:length(gfit))
0020 %NOTE - ORDER AS IS REQUIRED FOR DESCART - raster code will have to
0021 %disentangle in TauARaster_quick2.m
0022 
0023 % Remove any bad_rasters - e.g. if you have already calibrated and seen
0024 % dodgy data and want to re-do without it
0025 for b=1:length(new_bad_rasters)
0026     good_rasters=good_rasters(good_rasters~=new_bad_rasters(b))
0027 end
0028 
0029 
0030 good_rasters=good_rasters
0031 %good_rasters=good_rasters([1:8, 10:end])%:5 7:10 12:end])%10:end]); % don't use raster 9
0032 for i=1:length(good_rasters);
0033     data(i,1)=gfit{good_rasters(i),1}.a1; %I1
0034     data(i,2)=gfit{good_rasters(i),8}.a1; %I2
0035     data(i,3)=gfit{good_rasters(i),2}.a1; %Q1
0036     data(i,4)=gfit{good_rasters(i),4}.a1; %Q2
0037     data(i,5)=gfit{good_rasters(i),3}.a1; %U1
0038     data(i,6)=gfit{good_rasters(i),5}.a1; %U2
0039     data(i,7)= PA_D(good_rasters(i)); %Parallactic angle
0040 end
0041 %% Generate model data
0042 
0043 % Old - hardwired
0044 % Model data for TauA I=578Jy (3.4K), pol frac = 7%, Pol angle=38deg
0045 % Need to get a good reference for this.
0046 % if (strmatch(cal_source,'taua'))
0047 %     model = [578 0.07 38]; %
0048 % end
0049 %
0050 % if (strmatch(cal_source,'m42'))
0051 %     model = [415 0.0 0]; %
0052 % end
0053 
0054 %New - use same calculation as is already used for astro cal in pipeline
0055 
0056 mjd=tstr2mjd(start_date);
0057 [sourceNum sourceFlux fluxErr polang polfrac] = ...
0058     sourceCorrespondance({cal_source}, mjd);
0059 
0060 
0061 % Want to keep things in K (rather than Jy) so convert this sourceFlux:
0062 % read in telescope parameters from /constants/telescop_constants.m
0063 telescope_constants;
0064 
0065 Tsrc = sourceFlux.*(dishArea*apEff)./(2*k_b*1e26);
0066 
0067 model = [Tsrc polfrac polang]
0068     
0069 % Try finding out what the calibration matrix is for this data
0070 
0071 [calmat fval hessian]= polcal2(data, model);
0072 %[calmat resnorm residual jacobian]=polcal2(data,model)
0073 %calmat = [  0.9,  0.9,  0.0,  0.0, 0.0, 0.0;...
0074 %           -0.2, -0.3,  0.8,  0.7, 0.5, 0.6;...
0075 %           -0.3, -0.2, -0.4, -0.5, 0.9, 0.8];
0076 
0077 % % Now try applying this calibration data
0078 % for i=1:length(good_rasters)
0079 %     caldata(i,1:3) = data(i,1:3)*calmat; %corrected I,Q,U
0080 % end
0081 
0082 % Plot the results
0083 %close all
0084 
0085 % Plot the raw data
0086 %figure(1)
0087 figure (30)
0088 plot_labels = {'I1','I2','Q1','Q2','U1','U2'};
0089 legend_label = {'data','modelled data'};
0090 calmat
0091 for i=1:6 % loop over channels
0092 subplot(2,3,i);
0093 plot(data(:,7),data(:,i));
0094 hold all
0095 end
0096 
0097 
0098 %figure(2)
0099 % Plot the model*calmat --> observed data
0100 
0101 calmodel=zeros(size(data,1),6);
0102 for i=1:6;
0103     subplot(2,3,i);
0104     for npoints = 1:size(data,1); % loop over all parallactic angles and create 'observed' data from model at each
0105     phi = -data(npoints,7);
0106     modelvec = [model(1), model(1)*model(2)*cosd(2*(model(3)+phi)), model(1)*model(2)*sind(2*(model(3)+phi))];
0107     calmodel(npoints,:)=modelvec*calmat;
0108     end
0109 
0110 plot(data(:,7),calmodel(:,i))
0111 xlabel(plot_labels(i));
0112 legend(legend_label,'Location','SouthEast');
0113 if(i==1 | i==2)
0114     ylim([0 3.6])
0115 end
0116 end
0117 
0118 
0119 
0120 fits_name = strcat(start_date,'_',source_name,'_',IntensityType);
0121             
0122 save_file = strcat([maindir,'/',fits_name,'_caldata.mat']);
0123 plot_file1 =  strcat([maindir,'/',fits_name,'_polcal_plot.png']) ;
0124 plot_file2 =  strcat([maindir,'/',fits_name,'_polcal_plot2.png']) ;
0125 if(save_data==1)
0126  save(save_file,'calmat')
0127  print('-f30','-dpng',plot_file2)
0128 print('-f1','-dpng',plot_file1)
0129 end
0130     
0131 %close all