fit_IQUV2_radec

0001 function fit_IQUV2_radec(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=1;
0033 
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 % Calculate RA and DEC in degrees if they don't already exist
0067 % Use new function to get J2000 coords
0068   
0069 %if(~isfield(dgood.antenna0.servo, 'equa'))
0070       disp('Getting Equatorial co-ordinates in J2000');
0071       dgood = calcRADECJ2000(dgood);
0072 %  end
0073 %
0074 %   disp('Getting Equatorial co-ordinates');
0075 %   % we do not need to calculate the ra/dec until we get ready to write things
0076 %   % out.
0077 %   display('Calculating RA/DEC');
0078 %   long=-118.2822;
0079 %   lat=37.2339;
0080 %
0081 %   az = dgood.antenna0.servo.apparent(:,1);
0082 %   el = dgood.antenna0.servo.apparent(:,2);
0083 %   jd=mjd2jd(dgood.antenna0.receiver.utc);
0084 %   [equa] = horiz_coo([pi/180*(az) pi/180*(el)],jd,[pi/180*(long) ...
0085 %     pi/180*(lat)],'e');
0086 %   dgood.antenna0.servo.equa=equa;
0087 %
0088 %   clear az;
0089 %   clear el;
0090   
0091 for m= 1:length(s)
0092     
0093     n = n+1;
0094     ind = zeros(size(dgood.index.radio_point_scan.fast));
0095     ind(s(m):e(m)) = 1;
0096     ind = logical(ind);
0097     dcut  = framecut(dgood, ind);
0098     azApp = interp1(dcut.antenna0.tracker.utc, ...
0099         dcut.antenna0.tracker.horiz_topo(:,1), dcut.antenna0.receiver.utc);
0100     azOffSave = (azApp) - dcut.antenna0.servo.apparent(:,1);
0101     azOffSave = -azOffSave;
0102     azOffSave = azOffSave.*cos(dcut.antenna0.servo.el*pi/180);% Scale azimuth by cos(elevation)???CHECK
0103     
0104    
0105     elApp = interp1(dcut.antenna0.tracker.utc, ...
0106         dcut.antenna0.tracker.horiz_topo(:,2), dcut.antenna0.receiver.utc);
0107     elOffSave = elApp - dcut.antenna0.servo.apparent(:,2);
0108     elOffSave = -elOffSave;
0109     angle = sqrt(elOffSave.^2 + azOffSave.^2);
0110     
0111     % Make sure we ignore NaNs
0112     a = (isnan(azOffSave) | isnan(elOffSave)| isnan(dcut.antenna0.receiver.data(:,1))  );
0113     
0114     % Grab latitude and longitude of antenna
0115     
0116     [antenna_no antenna_name] = identifyTelescope(dcut);
0117     telescope_constants
0118     % Calculate the mean parallactic angle
0119        
0120     AZ = deg2rad(mean(azApp(~a)));
0121     EL = deg2rad(mean(elApp(~a)));
0122     
0123     
0124         
0125     % Calculate mean PA uses DEC at current epoch
0126     % Also need JD
0127     
0128     LAT_D = antenna_latitude(antenna_no)  % degrees
0129     LON_D = antenna_longitude(antenna_no) % degrees E
0130     
0131     LAT = deg2rad(LAT_D); % radians
0132     LON = deg2rad(LON_D); % radians
0133     
0134     JD = mjd2jd(dcut.antenna0.receiver.utc);
0135     LST=lst(JD,deg2rad(LON_D),'m'); % Uses East longitude in radians
0136                                     % Gives LST in fraction of days
0137     
0138     % Temporarily calculate current RA DEC
0139     dcut = calcRADECcurrent(dcut);
0140     RA_current = dcut.antenna0.servo.equa(:,1);  %radians
0141     DEC_current = dcut.antenna0.servo.equa(:,2); %radians
0142     
0143     [PA]=parallactic_angle((mean(RA_current)),(mean(DEC_current)),mean(LST),(LAT)); % Everything must be in radians, output is in radians
0144     PA_D(m) = rad2deg(PA); % gets definition in same way as real data. %% CHeck sign here - 13/8/2014 - ACT
0145     
0146     % Make sure that we store J2000 in the data structure
0147     dcut = calcRADECJ2000(dcut);
0148     
0149     
0150     
0151     % Form mean I,Q,U,V,P for this raster - Why?? Not sure this get used
0152     % anywhere any more
0153     
0154     I1(m)= mean(dcut.antenna0.receiver.data(~a,1));
0155     Q1(m) = mean(dcut.antenna0.receiver.data(~a,2));
0156     U1(m) = mean(dcut.antenna0.receiver.data(~a,3));
0157     Q2(m) = mean(dcut.antenna0.receiver.data(~a,4));
0158     U2(m) = mean(dcut.antenna0.receiver.data(~a,5));
0159     Q(m) = mean(dcut.antenna0.receiver.data(~a,6));
0160     U(m) = mean(dcut.antenna0.receiver.data(~a,7));
0161     I2(m)= mean(dcut.antenna0.receiver.data(~a,8));
0162     V(m)= mean(dcut.antenna0.receiver.data(~a,9));
0163     P(m) = sqrt(Q(m).^2 + U(m).^2);
0164 
0165     
0166     
0167   
0168    
0169 %  end
0170   
0171   % Calculate offset from nominal RA and DEC in degrees - need equa in
0172   % J2000 since TauA pos stored as J2000
0173 
0174 [RA_tau DEC_tau] = sourcePos(source_name); % deg
0175 RA_off = rad2deg(dcut.antenna0.servo.equa(:,1))-RA_tau; % degrees
0176 DEC_off = rad2deg(dcut.antenna0.servo.equa(:,2))-DEC_tau;
0177 RA_off = RA_off.*cosd(rad2deg(dcut.antenna0.servo.equa(:,2))); % scale RA to declination
0178 %pixel_size = 0.1;
0179 % RA = rad2deg(RA);
0180 % DEC= rad2deg(DEC);
0181 % ramin = RA_tau - 3;
0182 % ramax = RA_tau + 3;
0183 % decmin = DEC_tau - 3;
0184 % decmax = DEC_tau + 3;
0185 % ramin = - 3;
0186 % ramax = + 3;
0187 % decmin = - 3;
0188 % decmax = + 3;
0189     %Define a grid in degrees
0190 
0191 if(strmatch('N',cbass))
0192     ti=-3:cell_size:3;
0193  else
0194    ti =-10:cell_size:10;  % Kludge for South until pointing better
0195 end
0196     
0197     % Make sure we ignore NaNs
0198     a = (isnan(RA_off) | isnan(DEC_off) | isnan(dcut.antenna0.receiver.data(:,1)) );
0199     
0200     if(nchan==9)
0201         chantofit=[1,2,3,4,5,8]; % Do fits for I1,Q1,U1,I2,Q2,U2
0202         %chantofit=[2];
0203     else
0204             chan2fit=[1,2]
0205     end
0206     
0207     for (chan=chantofit); % only plot the I channel data for speed
0208         disp(['Fitting Gaussian ',num2str(m),' to channel: ',num2str(chan)])
0209         if(size(dcut.antenna0.receiver.data(~a,chan),1)>0)
0210             % Grid the data
0211             [XI,YI]=meshgrid(ti,ti);
0212             ZI = griddata(RA_off(~a),DEC_off(~a),dcut.antenna0.receiver.data(~a,chan),XI,YI);
0213             
0214             
0215             %Fit a gaussian to each one if you want
0216             
0217             if(fit_gauss)
0218             %disp('hello');
0219                 [sf gof t] = Egauss2D_fit_act(XI,YI,ZI);
0220                  gfit{m,chan} = sf;
0221                 % If you want to plot the fit
0222                 %plot(gfit{5},[XI(:),YI(:)],GoodI_c(:))
0223                 %clear sf gof t
0224             end
0225             
0226             
0227             if(plot_contours)
0228                 % Plot mini contour maps on 3x3 grid on each page
0229                 fig_no = ceil(m/16);
0230                 subplot_no = m-(fig_no-1)*16;
0231                 
0232                 %subplot(nfigs_x,nfigs_y,n);
0233                 %Make contours based on the max_value in the data set for I,Q,U,V
0234                 VI = [1:-0.1:-1]*contour_max(chan);
0235                 
0236                 figure(fig_no+2);
0237                 subplot(4,4,subplot_no)
0238                 contourf(XI,YI,ZI,VI);
0239                 axis square;
0240                 grid on;
0241                 xlabel('RA Offset,deg');
0242                 ylabel('DEC Offset,deg');
0243                 title(num2str(m));
0244             end
0245         end
0246     end
0247 end
0248 
0249 % Make sure we know where the good rasters are  - ie ones which have a
0250 % succesful gfit
0251 % This gets saved in ***_PA_IQUV.mat later
0252 
0253 empty=cellfun(@isempty,gfit);
0254 good_rasters=find(~empty(:,1))';
0255 
0256 % Save first page of plots as an example file
0257 if(savefileon)
0258     if(exist('fits_name'))
0259         map_name = strcat([maindir,'/',fits_name,'_Contour_RAdec.png'])
0260         %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig'])
0261     else
0262         fits_name = strcat(start_date,'_',source_name)
0263         map_name = strcat([maindir,'/',fits_name,'_Contour_RAdec.png'])
0264         %fig_name = strcat([maindir,'/',fits_name,'_Pointing.fig'])
0265     end
0266     figure(1)
0267     set(gcf, 'PaperPositionMode', 'auto');
0268     %set(gcf,'PaperType','A4');
0269     print('-f1','-dpng',map_name)
0270 end
0271 % Save the gfit, PA, IQUV and dgood data for re-use later without having to load and
0272 % calibrate again
0273 
0274 
0275 fits_name = strcat(start_date,'_',source_name)
0276 if(exist('calmat'))
0277     save_file = strcat([maindir,'/',fits_name,'_',IntensityType,'_PA_IQUV_cal_RAdec.mat'])
0278 else
0279     save_file = strcat([maindir,'/',fits_name,'_',IntensityType,'_PA_IQUV_RAdec.mat'])
0280 end
0281 
0282 save(save_file,'gfit','PA_D','I1','I2','Q','U','V','P','good_rasters')
0283 %close all
0284 end