Home > reduc > calcs > calib_beam_fits.m

calib_beam_fits

PURPOSE ^

find calibrator and fit a Gaussian over it WHERE the width of the

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 find calibrator and fit a Gaussian over it WHERE the width of the
 Gaussian is the width of the beam.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % find calibrator and fit a Gaussian over it WHERE the width of the
0002 % Gaussian is the width of the beam.
0003 
0004 % first apply alpha acorrectiona and correct for non linearity:
0005 
0006 d = determineIndices(d);
0007 d = apparentAzEl(d);
0008 
0009 
0010 if(~isfield(d.antenna0.servo, 'equa')) 
0011   display('Calculating RA/DEC');
0012   long=-118.2822;
0013   lat=37.2339;
0014   
0015   az = d.antenna0.servo.apparent(:,1);
0016   el = d.antenna0.servo.apparent(:,2);
0017   jd=mjd2jd(d.antenna0.receiver.utc);
0018   [equa] = horiz_coo([pi/180*(az) pi/180*(el)],jd,[pi/180*(long) ...
0019         pi/180*(lat)],'e');
0020   d.antenna0.servo.equa=equa;
0021   clear equa;
0022   clear az;
0023   clear el;
0024 end
0025 
0026 s = 'FILTERED';
0027 dO = assembleAlphaStreams(d,s); 
0028 [tA,Ac,Gc,Af,Gf] = calculateAlpha(dO,s); 
0029 dO = applyAlpha(dO,tA,Ac,Gc,Af,Gf,s); 
0030 d_out = calculateStokes(dO,s);
0031 
0032 
0033 d = calibrate_linearity(d);
0034 
0035 % find scan data
0036 
0037 d = framecut(d, bitsearch(d.array.frame.features, 7, 'any'));
0038 
0039 C1 = d.antenna0.receiver.data(:,1);
0040 C6 = d.antenna0.receiver.data(:,6);
0041 
0042 name = 'CasA'; % input your calibrator
0043 so = d.antenna0.tracker.source;
0044 H = find(ismember(so,name));
0045 H = [H(1): H(length(H)-1)];
0046 
0047 time = d.antenna0.receiver.utc;
0048 Time = (time - time(1)).*24*60;
0049 Az = d.antenna0.servo.apparent(:,1);
0050 El = d.antenna0.servo.apparent(:,2);
0051 InCo = [Az.*(pi/180) El.*(pi/180)];
0052 
0053 % need to get from Az/El to RA/Dec
0054 
0055 ra = zeros(length(time),1);
0056 dec = zeros(length(time),1);
0057 
0058 for i = 1: length(time)
0059 
0060     lat = 37.2313;
0061     lon = -118.283; 
0062     utc = time(i); 
0063     JD=mjd2jd(utc);
0064     outCo = horiz_coo(InCo(i,1:2), JD, [lon*(pi/180) lat*(pi/180)], 'e'); 
0065     ra(i) = outCo(1,1).*180/pi;
0066     dec(i) = outCo(1,2).*180/pi;
0067     
0068 end
0069 
0070 % cut data to only include named calibrator
0071 
0072 LD = length(C1);
0073 
0074 Pp = [1:99:LD];
0075 Gg = [2:LD/100];
0076 
0077 
0078 Cc = Pp(Gg);
0079 Ee = Pp(Gg-1);
0080 
0081 KK = (LD/100) - 1;
0082 
0083 for kk = 1 : floor(KK)
0084     
0085     AC1(1,kk) = mean(C1(Ee(kk): Cc(kk)));
0086     AC6(1,kk) = mean(C6(Ee(kk): Cc(kk)));
0087     A2(1,kk) = mean(ra(Ee(kk) : Cc(kk)));
0088     A3(1,kk) = mean(dec(Ee(kk) : Cc(kk)));
0089     A4(1,kk) = mean(Time(Ee(kk) : Cc(kk)));
0090     A5(1,kk) = mean(El(Ee(kk) : Cc(kk)));
0091     A6(1,kk) = mean(Az(Ee(kk) : Cc(kk)));
0092            
0093 end 
0094 
0095 D1 = AC1';
0096 D6 = AC6';
0097 ra = A2';
0098 dec = A3';
0099 Time = A4'; 
0100 El = A5';
0101 Az = A6';
0102 
0103     newTime = Time(H);
0104     newRa = ra(H);
0105     newDec = dec(H);
0106     newD1 = D1(H);
0107     newD6 = D6(H);  
0108     newEl = El(H);
0109     newAz = Az(H);
0110           
0111 % now to get from RA/Dec to offsets
0112 
0113 if strcmp(name, 'TauA') == 1
0114    ra1 = 83.6333;
0115    dec1 = 22.2167;
0116 elseif strcmp(name, 'CygA') == 1
0117        ra1 = 299.8679;
0118        dec1 = 40.7339;
0119 elseif strcmp(name, 'CasA') == 1
0120        ra1 = 350.8500;
0121        dec1 = 58.8150;
0122 end
0123 
0124 newDec = newDec - dec1; % turn everything into offsets from ra1, dec1
0125 newRa = newRa - ra1;
0126 
0127 sky_dec = newDec;
0128 sky_ra = newRa;
0129 
0130 
0131 
0132     % make sure we haven't got a scan tail at wrong ra and dec as we move in from another source
0133     
0134     list_ra = zeros(length(sky_ra));
0135         
0136     for i = 1 : length(sky_ra)
0137         if sky_ra(i) > 5 
0138            sky_ra(i) = NaN;
0139            sky_dec(i) = NaN;
0140            newD1(i) = NaN;
0141            newD6(i) = NaN;
0142            list_ra(i) = 1;
0143         elseif sky_ra(i) < -5
0144                sky_ra(i) = NaN;
0145                sky_dec(i) = NaN;
0146                newD1(i) = NaN;
0147                newD6(i) = NaN;
0148                list_ra(i) =1;
0149         end
0150     end
0151     
0152  % try actaully fittin the gaussians
0153 
0154     xx1 = find(newD1 == nanmax(newD1));
0155     x_01_x = sky_ra(xx1); % x position
0156     x_01_y = sky_dec(xx1); % y position
0157     x_01_x = x_01_x(1);
0158     x_01_y = x_01_y(1);
0159     
0160     xx6 = find(newD6 == nanmax(newD6));
0161     x_06_x = sky_ra(xx6); % x position
0162     x_06_y = sky_dec(xx6); % y position
0163     x_06_x = x_06_x(1);
0164     x_06_y = x_06_y(1); 
0165     
0166  
0167 L = length(newD1); 
0168 
0169 % to find offset gradient by fitting a straight line to the non Gaussian
0170 % part of scan
0171 
0172 off1 = msbackadj_copied(newTime, newD1, 'WindowSize',1,'StepSize',2);
0173 straight1 = newD1-off1;
0174 ascS1 = sort(straight1);
0175 off6 = msbackadj_copied(newTime, newD6, 'WindowSize',1,'StepSize',2);
0176 straight6 = newD6-off6;
0177 
0178 if list_ra(1) == 0 && list_ra(L) == 0
0179 
0180     a = straight1(1);
0181     b = straight1(L);
0182     c = straight6(1);
0183     d = straight6(L);
0184     e = sky_ra(1);
0185     f = sky_ra(L);
0186     g = sky_dec(1);
0187     h = sky_dec(L);
0188 else
0189     a = straight1(10);
0190     b = straight1(L-15);
0191     c = straight6(10);
0192     d = straight6(L-15);
0193     e = sky_ra(10);
0194     f = sky_ra(L-15);
0195     g = sky_dec(10);
0196     h = sky_dec(L-15);
0197 end
0198 
0199 fit1 = cat(1, a, b);
0200 fit6 = cat(1,c,d);
0201 ra_cut = cat(1,e,f);
0202 dec_cut = cat(1,g,h);
0203  
0204 Px1 = polyfit(ra_cut, fit1,1);
0205 Py1 = polyfit(dec_cut, fit1,1);
0206 grad1X = Px1(1);
0207 grad1Y = Py1(1);
0208 
0209 Px6 = polyfit(ra_cut, fit6,1);
0210 Py6 = polyfit(dec_cut, fit6,1);
0211 grad6X = Px6(1);
0212 grad6Y = Py6(1);
0213 
0214 asc1 = sort(newD1);
0215 offset1 = asc1(1);
0216 asc6 = sort(newD6);
0217 offset6 = asc6(1);
0218 
0219 A1 = max(newD1) - offset1; % Gaussian peak amplitudes
0220 A6 = max(newD6) - offset6;
0221 
0222 Sigma_fit = 0.73/(2*sqrt(2*log(2)));   
0223 % This is cos we scan at 0.5 deg per second for calib scans
0224 Sigma_x_fit = Sigma_fit/0.5;
0225 Sigma_y_fit = Sigma_fit/0.5;
0226 
0227 Sigma_fit_tot = sqrt((Sigma_x_fit)^2 + (Sigma_y_fit)^2);
0228                                                                                                              
0229 BBeta1 = [A1, x_01_x, Sigma_x_fit, x_01_y, Sigma_y_fit, grad1X, grad1Y, offset1];
0230 BBeta6 = [A6, x_06_x, Sigma_x_fit, x_06_y, Sigma_y_fit, grad1Y, grad6Y, offset6];
0231 
0232 xy1 = cat(1, sky_ra, sky_dec);
0233 [betanew1,resid1, J1, Cov1]=nlinfit(xy1,newD1,@gauss2D_fit,BBeta1);
0234 [betanew6,resid6,J6,Cov6]=nlinfit(xy1,newD6,@gauss2D_fit,BBeta6);
0235 
0236 Z1 = gauss2D_draw(betanew1, sky_ra, sky_dec);
0237 Z6 = gauss2D_draw(betanew6, sky_ra, sky_dec); 
0238 
0239 figure(1)
0240 mesh(sky_ra,sky_dec,Z1)
0241 colormap(copper)
0242 hold on;
0243 plot3(sky_ra, sky_dec, newD1, 'o')
0244 hold off
0245 figure(2)
0246 mesh(sky_ra,sky_dec,Z6);
0247 colormap(copper)
0248 hold on;
0249 plot3(sky_ra, sky_dec, newD6, 'o')
0250 hold off;
0251 
0252 
0253 Channel1_Amp = betanew1(1)
0254 Channel6_Amp = betanew6(1)
0255 ra_offset1 = betanew1(2)
0256 dec_offset1 = betanew1(4)
0257 ra_offset6 = betanew6(2)
0258 dec_offset6 = betanew6(4)
0259 El1 = newEl(xx1)
0260 El6 = newEl(xx6)
0261 
0262 amp_error1 = (Cov1(1,1)/betanew1(1))*100
0263 amp_error6 = (Cov6(1,1)/betanew6(1))*100
0264
Generated on Sun 14-Jun-2015 17:12:45 by m2html © 2005