fit_IQUV2

0001 function fit_IQUV2(nchan, plot_labels, cbass,IntensityType,save_path,start_date,source_name,dgood,calmat)
0002 
0003 if(~exist('dgood'))
0004     
0005     %save_path ='/test_fits_June2013/' % Adds this on to path of cbass_analysis
0006     [home,installeddir] = where_am_i();
0007     addpath([home,'/',installeddir,'/Angelas_Raster_Code/']);
0008     
0009 %     addpath([home,'/',installeddir,'/angela_cbass/angela_beam/sim_beams/']);
0010 %     addpath([home,'/',installeddir,'/angela_cbass/angela_beam/Beam_plotting/'])
0011 %     addpath([home,'/',installeddir,'/angela_cbass/angela_beam/rasters/'])
0012 %     addpath([home,'/',installeddir,'/angela_cbass/angela_beam/Tsys_stuff/'])
0013 %     addpath([home,'/',installeddir,'/angela_cbass/angela_beam/reducs/'])
0014 
0015 maindir= [home,'/',installeddir,'/',save_path,'/'];
0016     
0017     load([maindir,start_date,'_',source_name,'_',IntensityType,'_dgood.mat'])
0018 end
0019 
0020 %%%% NEED TO CHANGE SIZE OF MATRIX HERE %%%
0021 if(exist('calmat'))
0022     for i=1:length(dgood.antenna0.receiver.data(:,1))
0023         dgood.antenna0.receiver.data(i,1:3)=dgood.antenna0.receiver.data(i,1:3)*calmat;  % Mulitply I,Q,U by the calibration matrix
0024     end
0025 end
0026 %%%%%
0027 
0028 cell_size= 0.1;
0029 fit_gauss= 1;
0030 n=0;
0031 plot_contours=0;
0032 savefileon=0;
0033 gfit=[];
0034 [home,installeddir]=where_am_i();
0035 %save_path ='/test_fits/';
0036 maindir= [home,'/',installeddir,'/',save_path,'/'];
0037 
0038 % Now make contour maps if you wish
0039 
0040 %Identify the raster data - changed flag after 1 Dec
0041 % if dgood.antenna0.receiver.utc(1) > tstr2mjd('01-Dec-2012:07:46:43')
0042 %    flagname = 'beammap';
0043 %    dgood.index.radio_point_scan.fast = dgood.index.beammap.fast;
0044 %    dgood.index.radio_point_scan.slow = dgood.index.beammap.slow;
0045 % else
0046 flagname = 'radio_point_scan'
0047 
0048 % end
0049 
0050 % Find start and end points of each raster
0051 [s e] = findStartStop(dgood.index.radio_point_scan.fast);
0052 
0053 %Identify the maximum of I,Q,U,V to use as a normalization for each raster
0054 
0055 points = (dgood.index.radio_point_scan.fast==1);
0056 for i=1:nchan
0057     contour_max(i) = max(dgood.antenna0.receiver.data(points,i));
0058 end
0059 
0060 %disp(['Will now fit gaussians to each I,Q,U,V raster. ', num2str(length(s)),' to do...'])
0061 
0062 % ACT 16/1/2014 - only fit gaussians to I1 and I2  (or LL and RR) - use them for pointing
0063 % correction - shoudl be same in both
0064 disp(['Will now fit gaussians to each I1 and I2 raster. ', num2str(length(s)),' to do...'])
0065 
0066 for m= 1:length(s)
0067     
0068     n = n+1;
0069     ind = zeros(size(dgood.index.radio_point_scan.fast));
0070     ind(s(m):e(m)) = 1;
0071     ind = logical(ind);
0072     dcut  = framecut(dgood, ind);
0073     azApp = interp1(dcut.antenna0.tracker.utc, ...
0074         dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc);
0075     azOffSave = (azApp) - dcut.antenna0.servo.apparent(:,1);
0076     azOffSave = -azOffSave;
0077     azOffSave = azOffSave.*cos(dcut.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK
0078     
0079     
0080     elApp = interp1(dcut.antenna0.tracker.utc, ...
0081         dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc);
0082     elOffSave = elApp - dcut.antenna0.servo.apparent(:,2);
0083     elOffSave = -elOffSave;
0084     angle = sqrt(elOffSave.^2 + azOffSave.^2);
0085     
0086     % Make sure we ignore NaNs
0087     a = (isnan(azOffSave) | isnan(elOffSave)| isnan(dcut.antenna0.receiver.data(:,1))  );
0088     
0089     % Calculate the mean parallactic angle
0090            
0091     % Grab latitude and longitude of antenna
0092     
0093     [antenna_no antenna_name] = identifyTelescope(dcut);
0094     telescope_constants
0095     
0096         
0097     % Calculate mean PA uses DEC at current epoch
0098     % Also need JD
0099     
0100     LAT_D = antenna_latitude(antenna_no)  % degrees
0101     LON_D = antenna_longitude(antenna_no) % degrees E
0102     
0103     LAT = deg2rad(LAT_D); % radians
0104     LON = deg2rad(LON_D); % radians
0105     
0106     JD = mjd2jd(dcut.antenna0.receiver.utc);
0107     LST=lst(JD,deg2rad(LON_D),'m'); % Uses East longitude in radians
0108                                     % Gives LST in fraction of days
0109     
0110     % Temporarily calculate current RA DEC
0111     dcut = calcRADECcurrent(dcut);
0112     RA_current = dcut.antenna0.servo.equa(:,1);  %radians
0113     DEC_current = dcut.antenna0.servo.equa(:,2); %radians
0114     
0115     [PA]=parallactic_angle(((RA_current)),((DEC_current)),(LST),(LAT)); % Everything must be in radians, output is in radians
0116     PA_D(m) = rad2deg(mean(PA)); % gets definition in same way as real data. %% CHeck sign here - 13/8/2014 - ACT
0117     
0118     % Make sure that we store J2000 in the data structure
0119     dcut = calcRADECJ2000(dcut);
0120     
0121     
0122     
0123 
0124 %%%%
0125 %     % OLD - now use parallactic_angle function same as in FakeRaster.m -
0126 %     % ACt 13/8/2014
0127 %    DEC = deg2rad(mean((dcut.antenna0.tracker.equat_geoc(:,2))));
0128 %    sinPA = -((sin(AZ).*cos(LAT))./cos(DEC)); % radians
0129 %    sinPA(abs(sinPA)>=1)=1;
0130 %    PA = asin(sinPA);
0131 %    PA_D(m)= rad2deg(PA); % deg
0132 %%%%%
0133 
0134     % Form mean I,Q,U,V,P for this raster - Why?? Not sure this get used
0135     % anywhere any more
0136     
0137     I1(m)= mean(dcut.antenna0.receiver.data(~a,1));
0138     Q1(m) = mean(dcut.antenna0.receiver.data(~a,2));
0139     U1(m) = mean(dcut.antenna0.receiver.data(~a,3));
0140     Q2(m) = mean(dcut.antenna0.receiver.data(~a,4));
0141     U2(m) = mean(dcut.antenna0.receiver.data(~a,5));
0142     Q(m) = mean(dcut.antenna0.receiver.data(~a,6));
0143     U(m) = mean(dcut.antenna0.receiver.data(~a,7));
0144     I2(m)= mean(dcut.antenna0.receiver.data(~a,8));
0145     V(m)= mean(dcut.antenna0.receiver.data(~a,9));
0146     P(m) = sqrt(Q(m).^2 + U(m).^2);
0147     
0148     %Define a grid in degrees
0149 
0150 if(strmatch('N',cbass))
0151     ti=-3:cell_size:3;
0152  else
0153    ti =-10:cell_size:10;  % Kludge for South until pointing better
0154 end
0155     
0156     % Make sure we ignore NaNs
0157     a = (isnan(azOffSave) | isnan(elOffSave) | isnan(dcut.antenna0.receiver.data(:,1)) );
0158     
0159     if(nchan==9)
0160         chantofit=[1,2,3,4,5,8]; % Do fits for I1,Q1,U1,I2,Q2,U2
0161     else
0162             chantofit=[1,2]
0163     end
0164     
0165     for (chan=chantofit); % only plot the I channel data for speed
0166         disp(['Fitting Gaussian ',num2str(m),' to channel: ',num2str(chan)])
0167         if(size(dcut.antenna0.receiver.data(~a,chan),1)>0)
0168             % Grid the data
0169             [XI,YI]=meshgrid(ti,ti);
0170             ZI = griddata(azOffSave(~a),elOffSave(~a),dcut.antenna0.receiver.data(~a,chan),XI,YI);
0171             imagesc(ZI)
0172             
0173             %Fit a gaussian to each one if you want
0174             
0175             if(fit_gauss)
0176             %disp('hello');
0177                 fit_params = [];
0178                 [sf gof t] = Egauss2D_fit_act(XI,YI,ZI,fit_params,chan,gfit,m);
0179                 gfit{m,chan} = sf;
0180                 % If you want to plot the fit
0181                 %plot(gfit{5},[XI(:),YI(:)],GoodI_c(:))
0182                 %clear sf gof t
0183             end
0184             
0185             
0186             if(plot_contours)
0187                 % Plot mini contour maps on 3x3 grid on each page
0188                 fig_no = ceil(m/16);
0189                 subplot_no = m-(fig_no-1)*16;
0190                 
0191                 %subplot(nfigs_x,nfigs_y,n);
0192                 %Make contours based on the max_value in the data set for I,Q,U,V
0193                 VI = [1:-0.1:0]*contour_max(chan);
0194                 
0195                 figure(fig_no+2);
0196                 subplot(4,4,subplot_no)
0197                 contour(XI,YI,ZI,VI);
0198                 axis square;
0199                 grid on;
0200                 xlabel('Az Offset,deg');
0201                 ylabel('El Offset,deg');
0202                 title(num2str(m));
0203             end
0204         end
0205     end
0206 end
0207 
0208 % Make sure we know where the good rasters are  - ie ones which have a
0209 % succesful gfit
0210 % This gets saved in ***_PA_IQUV.mat later
0211 
0212 empty=cellfun(@isempty,gfit);
0213 good_rasters=find(~empty(:,1))';
0214 
0215 % Save first page of plots as an example file
0216 if(savefileon)
0217     if(exist('fits_name'))
0218         map_name = strcat([maindir,'/',fits_name,'_Contour.png'])
0219         %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig'])
0220     else
0221         fits_name = strcat(start_date,'_',source_name)
0222         map_name = strcat([maindir,'/',fits_name,'_Contour.png'])
0223         %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig'])
0224     end
0225     figure(1)
0226     set(gcf, 'PaperPositionMode', 'auto');
0227     %set(gcf,'PaperType','A4');
0228     print('-f1','-dpng',map_name)
0229 end
0230 % Save the gfit, PA, IQUV and dgood data for re-use later without having to load and
0231 % calibrate again
0232 
0233 
0234 fits_name = strcat(start_date,'_',source_name)
0235 if(exist('calmat'))
0236     save_file = strcat([maindir,'/',fits_name,'_',IntensityType,'_PA_IQUV_cal.mat'])
0237 else
0238     save_file = strcat([maindir,'/',fits_name,'_',IntensityType,'_PA_IQUV.mat'])
0239 end
0240 
0241 save(save_file,'gfit','PA_D','I1','I2','Q','U','V','P','good_rasters')
0242 close all
0243 end