lukesTrigPointGaussianFitter

0001 addpath('./lukes_code')
0002 addpath('./Angelas_Raster_Code')
0003 %% Fit gauss to trig points
0004 %% Select which chunks of data you want to use
0005 % 0dbm - 60db, osc freq of 10Hz, centre freq of 5.104GHz
0006 start_dates = ['10-May-2014:09:28:32'; '10-May-2014:09:56:44'; '10-May-2014:10:15:04'; '10-May-2014:10:34:45'];
0007 end_dates = ['10-May-2014:09:44:02'; '10-May-2014:10:12:14'; '10-May-2014:10:30:34'; '10-May-2014:10:50:15'];
0008 angles = [0; +45; -45; 90 ];
0009 angles_string = {'0 degrees', '+45 degrees', '+90 degrees', '-45 degrees' };
0010 
0011 %% Parameter Choices
0012 % Put data into CBASS-N style
0013 chan=78;
0014 nchan=4;
0015 plot_labels = {'LL','RR','Q','U'};
0016 LL_index = 1;     % cbassS2N spits out d.data(:,1) = LL ; 6 = RR ; 2 = QQ ; 3 = UU
0017 RR_index = 6;
0018 QQ_index = 2;
0019 UU_index = 3;
0020 % Calculate the az and El offset form the source
0021 az_trig=16;
0022 el_trig=6.5;
0023 % Fitting the Gaussians
0024 plot_beamMaps_ = 1;
0025 fit_gaussian_ = 1;
0026 overlay_gaussian_ = 1;
0027 plot_cuts_ = 1;
0028 save_plots_ = 1;
0029 pixel_size = 0.1; % degree
0030 colLR = 1.5e5;
0031 colQU = 1e4;
0032 % FWHM calculation
0033 FWHMaz = zeros(4);
0034 FWHMel = zeros(4);
0035 
0036 %% Run the code
0037 for j=1:4
0038     d = read_arcSouth( char(start_dates(j,:)), char(end_dates(j,:)) )
0039     
0040     % Put data into Cbass-N style
0041     disp( 'Puttiing data into CBASS-N style' )
0042     d = cbassS2N(d,chan);
0043      
0044     % Calculate the az and El offset form the source
0045     disp( 'Calculating az and el offsets from the source' )
0046     d.azOffSave = (d.antenna0.servo.apparent(:,1)-az_trig).*cosd(d.antenna0.servo.apparent(:,2));
0047     d.elOffSave = d.antenna0.servo.apparent(:,2)-el_trig;
0048     d.angle = sqrt(d.elOffSave.^2 + d.azOffSave.^2);
0049     
0050     % Separate data in to raster, elscan and azscan
0051     disp( 'Separating data into raster, elscan and azscan' )
0052     draster = framecut(d,d.index.beammap.fast);
0053     dazscan = framecut(d,d.index.radio_point_scan.fast);
0054     delscan = framecut(d,d.index.elscan.fast);
0055     
0056     % NaN Flagging
0057     disp( 'Flagging NaNs' )
0058     nanflag = isnan(draster.antenna0.receiver.data(:,1));
0059     flag = ~nanflag;
0060     
0061     % Demodulate the relevant parameters
0062     disp( 'Demodulating the data' )
0063     
0064     ichoice1 = demodStartIndexFinder( draster.antenna0.receiver.data(flag,LL_index) );
0065     ichoice2 = demodStartIndexFinder( draster.antenna0.receiver.data(flag,RR_index) );
0066     ichoice3 = demodStartIndexFinder( draster.antenna0.receiver.data(flag,QQ_index) );
0067     ichoice4 = demodStartIndexFinder( draster.antenna0.receiver.data(flag,UU_index) );
0068     
0069     ichoice = round( mean( [ichoice1 ichoice2 ] ) );
0070     
0071     LL_demod = demodBG( ichoice, draster.antenna0.receiver.data(flag,LL_index) );
0072     RR_demod = demodBG( ichoice, draster.antenna0.receiver.data(flag,RR_index) );
0073     QQ_demod = demodBG( ichoice, draster.antenna0.receiver.data(flag,QQ_index) );
0074     UU_demod = demodBG( ichoice, draster.antenna0.receiver.data(flag,UU_index) );
0075     dataRaster_demod = [ LL_demod', RR_demod', QQ_demod', UU_demod' ];
0076     
0077     azOffRaster_demod = demodNoBG( ichoice, draster.azOffSave(flag) );
0078     elOffRaster_demod = demodNoBG( ichoice, draster.elOffSave(flag) );
0079     
0080     % Fitting the Gaussians and plotting the rasters
0081     azmin = min(azOffRaster_demod);
0082     azmax = max(azOffRaster_demod);
0083     elmin = min(elOffRaster_demod);
0084     elmax = max(elOffRaster_demod);
0085     
0086     for i=1:4
0087         [TauAmap{i}, xidx, yidx, TauAmapIdx{i}] = bin_quickly2d(azmin, azmax, elmin, elmax, pixel_size,...
0088             azOffRaster_demod, elOffRaster_demod,dataRaster_demod(:,i));
0089         [XI, YI] = meshgrid( xidx, yidx );
0090         
0091     if plot_beamMaps_ == 1
0092         disp( 'Plotting Beam Maps' )
0093         figure( j*10 )
0094             suptitle( angles_string{j} )
0095             subplot(4,1,i)
0096             if i<3
0097                 clim=[0 colLR];
0098             else
0099                     clim=[-colQU colQU];
0100             end
0101             imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i},clim)
0102         colormap( 'default' )
0103             title(plot_labels{i})
0104             axis equal
0105             ylim( [elmin elmax] );
0106             xlim( [azmin azmax] );
0107             colorbar
0108     end
0109 
0110         if fit_gaussian_ == 1
0111             % Gaussian Fitting
0112             disp( 'Fiting Gaussian' )
0113             % Initial starting guess
0114             %[fitMaxVal1, fitMaxInd1] = max( abs(TauAmap{i}) );
0115             %[fitMaxVal2, fitMaxInd2] = max( fitMaxVal1 );
0116             %[sf gof t] = Egauss2D_fit_act(XI,YI,TauAmap{i});
0117             % Want to fit anegative gaussian to +90 Q and -45 U
0118             if i==j & i>2
0119                 [sf gof t] = Egauss2D_fit_negative_act(XI,YI,TauAmap{i});
0120             else
0121                 [sf gof t] = Egauss2D_fit_act(XI,YI,TauAmap{i});
0122             end
0123             gfit{i, j} = sf;
0124             gfit_gof{i, j} = gof;
0125             a1 = gfit{i, j}.a1;
0126             theta = gfit{i, j}.theta;
0127             sigma_x = gfit{i, j}.sigma_x;
0128             sigma_y = gfit{i, j}.sigma_y;
0129             x0 = gfit{i, j}.x0;
0130             y0 = gfit{i, j}.y0;
0131             z = a1*exp( - ( (cos(theta)^2/2/sigma_x^2 + sin(theta)^2/2/sigma_y^2)*(XI-x0).^2 + 2*(-sin(2*theta)/4/sigma_x^2 + sin(2*theta)/4/sigma_y^2)*(XI-x0).*(YI-y0) + (sin(theta)^2/2/sigma_x^2 + cos(theta)^2/2/sigma_y^2)*(YI-y0).^2  ));
0132             FWHMaz( i, j ) = 2.35482*sigma_x;
0133             FWHMel( i, j ) = 2.35482*sigma_y;
0134         end
0135         
0136         if overlay_gaussian_ ==1
0137             disp( 'Plotting Fit' )
0138             figure(j*10+1)
0139             suptitle( angles_string{j} )
0140             subplot(4,1,i)
0141             if i<3
0142                 clim=[0 colLR];
0143             else
0144                 clim=[-colQU colQU];
0145             end
0146             imagesc(fliplr(xidx),fliplr(yidx),TauAmap{i},clim)
0147             hold;
0148             contour(fliplr(xidx),fliplr(yidx),z,'LineColor', [0 0 0])
0149             colormap('default')
0150             title(plot_labels{i})
0151             axis equal
0152             ylim( [elmin elmax] );
0153             xlim( [azmin azmax] );
0154             colorbar
0155         end
0156         
0157         if plot_cuts_ == 1
0158             disp( 'Plotting cuts through the fit(blue) and data(red)' )
0159             if i==j & i>2
0160                 [peakvak, peakindex] = min( min(z) );
0161             else
0162                 [peakval, peakindex] = max( max(z) );
0163             end
0164             figure( j*10+2)
0165             suptitle( angles_string{j} )
0166             subplot(4, 1, i)
0167             plot( fliplr(yidx), TauAmap{i}(:,peakindex), '.r' )
0168             hold
0169             plot( fliplr(yidx), z(:,peakindex),'b' )
0170             ylabel( plot_labels{i} )
0171             xlabel( 'Elevation Offset' )
0172         end
0173     end
0174     
0175     if save_plots_ ==1
0176         if j==1
0177             saveas(10, 'plus0degBeamMap','epsc')
0178             saveas(11, 'plus0degFitCountour','epsc')
0179             saveas(12, 'plus0degFitCuts','epsc')
0180         end
0181         if j==2
0182                saveas(20, 'plus45degBeamMap','epsc')
0183             saveas(21, 'plus45degFitContour','epsc')
0184             saveas(22, 'plus45degFitCuts','epsc')
0185         end
0186         if j==3
0187                saveas(30, 'plus90degBeamMap','epsc')
0188             saveas(31, 'plus90degFitContour','epsc')
0189             saveas(32, 'plus90degFitCuts','epsc')
0190         end
0191         if j==4
0192                saveas(40, 'minus45degBeamMap','epsc')
0193             saveas(41, 'minus45degFitContour','epsc')
0194             saveas(42, 'minus45degFitCuts','epsc')
0195         end
0196         close all
0197     end
0198 end
0199     
0200 save( 'FWHM.mat', 'FWHMaz','FWHMel' )