absCBF_AZEL

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function calibrator vals = calib_beam_fitsAZEL(d)

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

0001   function [utc, source_num, Amp, newAmp, el_off, az_off, tot_off, SigAz, ...
0002             SigEl, EL, chi, error, SkyDip, Theo, Pipe, pang, AZ]  = absCBF_AZEL(d, popt)
0003 
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 %  function calibrator vals = calib_beam_fitsAZEL(d)
0007 %
0008 % find calibrator and fit a Gaussian over it WHERE the width of the
0009 % Gaussian is the width of the beam.
0010 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0011 
0012 % first run pipeline
0013 
0014 d = reduceData_absCal(d, 'redScript.red', 0, 'testtest2');
0015 
0016 % get opacity data
0017 
0018 skydip = NaN;
0019 SkyDipMethod2v2
0020 skydip = OpacityOne;
0021 TheoOp = CalcTheoOp(d);
0022 pipe = calculateTau(d);
0023 pipeTheo = last(pipe(:,4));
0024 
0025 % find scan data
0026 
0027 d = framecut(d, bitsearch(d.array.frame.features, 7, 'any'));
0028 so = d.antenna0.tracker.source;
0029 calibs = {'CasA', 'CygA', 'TauA'}; 
0030 tally = ismember(calibs, so);
0031 
0032 if tally(1) == 1
0033    name = {'CasA'};
0034 end
0035 if tally(2) == 1
0036    name = {'CygA'};
0037 end
0038 if tally(3) ==1
0039    name = {'TauA'};
0040 end
0041 if tally(2) && tally(3) ==1
0042    name = {'CygA', 'TauA'};
0043 end
0044 if tally(1) && tally(3) ==1
0045    name = {'TauA','CasA'};
0046 end
0047 if tally(1) && tally(2) ==1
0048    name = {'CygA', 'CasA'};  
0049 end
0050 
0051 if length(name) == 2
0052    repeat = 2;
0053 else repeat = 1;
0054 end
0055 
0056 for hhh = 1:repeat
0057  
0058     name1 = name(hhh);
0059 
0060 if strcmp(name1, 'CasA')
0061    name_num = 1; 
0062 elseif strcmp(name1, 'CygA')
0063        name_num = 2; 
0064 elseif strcmp(name1, 'TauA')
0065        name_num = 3;
0066 end
0067 
0068 test = strcmp(name1, so);
0069 testData = framecut(d, test);
0070 d1 = testData;
0071 D1 = d1.antenna0.receiver.data(:,1);
0072 time = d1.antenna0.receiver.utc;
0073 Time = (time - time(1)).*24*60;
0074 Az = d1.antenna0.servo.apparent(:,1);
0075 El = d1.antenna0.servo.apparent(:,2);
0076 LD = length(D1);
0077 
0078 AzSource = interp1(d1.antenna0.tracker.utc, d1.antenna0.tracker.horiz_topo(:,1), d1.antenna0.receiver.utc);
0079 ElSource = interp1(d1.antenna0.tracker.utc, d1.antenna0.tracker.horiz_topo(:,2), d1.antenna0.receiver.utc);
0080 programmed_Az_offset = interp1(d1.antenna0.tracker.utc, d1.antenna0.tracker.scan_off(:,1), d1.antenna0.receiver.utc);
0081 programmed_El_offset = interp1(d1.antenna0.tracker.utc, d1.antenna0.tracker.scan_off(:,2), d1.antenna0.receiver.utc);
0082 
0083 %replace NaN values with their following values
0084 AzN = isnan(AzSource);
0085 LAzN = length(find(AzN));
0086 AzOffN = isnan(programmed_Az_offset);
0087 LAzOffN = length(find(AzOffN));
0088 
0089 AzSource(1:LAzN) = repmat(AzSource(LAzN+1), 1, LAzN);
0090 ElSource(1:LAzN) = repmat(ElSource(LAzN+1), 1, LAzN);
0091 programmed_Az_offset(1:LAzOffN) = repmat(programmed_Az_offset(LAzOffN+1), 1, LAzOffN);
0092 programmed_El_offset(1:LAzOffN) = repmat(programmed_El_offset(LAzOffN+1), 1, LAzOffN);
0093 
0094 
0095     newTime = Time;
0096     newD1 = D1;
0097     newEl = El;
0098     newAz = Az;
0099     new_time = time;  
0100     new_AzSource = AzSource;
0101     new_ElSource = ElSource;
0102     new_prog_El_offset = programmed_El_offset;
0103     new_prog_Az_offset = programmed_Az_offset;
0104     
0105     % make sure we don't have scan tail as telescope slews from another
0106     % source
0107    
0108     pop = last(newAz);
0109     inc = pop * 1.15;
0110     content = find(newAz < inc);
0111    
0112     newTime = newTime(content);
0113     newD1 = newD1(content);
0114     newEl = newEl(content);
0115     newAz = newAz(content);
0116     new_time =   new_time(content);
0117     new_AzSource = new_AzSource(content);
0118     new_ElSource = new_ElSource(content);
0119     new_prog_El_offset = new_prog_El_offset(content);
0120     new_prog_Az_offset = new_prog_Az_offset(content);
0121 
0122    
0123     % there's more than one Az/El scan for a source so we need to split them up
0124     % and get their offsets
0125     
0126     [maxtab] = peakdet(newAz, 4, newTime);
0127     Div = length(maxtab) -1; 
0128     newL = floor(length(newAz)/Div);
0129     
0130     for i = 1 : Div
0131         sky_El(i+(i-1)*newL-(i-1)*1:newL*i) = newEl(i+(i-1)*newL-(i-1)*1:newL*i) - new_ElSource(i+(i-1)*newL-(i-1)*1:newL*i);
0132         sky_Az(i+(i-1)*newL-(i-1)*1:newL*i) = (newAz(i+(i-1)*newL-(i-1)*1:newL*i) - new_AzSource(i+(i-1)*newL-(i-1)*1:newL*i)).*cosd(newEl(i+(i-1)*newL-(i-1)*1:newL*i));
0133     end
0134     
0135     
0136 % now lets fit all the scans
0137 
0138 for j = 1 : Div
0139     
0140     sky_El_a = sky_El(j+(j-1)*newL-(j-1)*1:newL*j);
0141     sky_Az_a = sky_Az(j+(j-1)*newL-(j-1)*1:newL*j);
0142     newD1_a = newD1(j+(j-1)*newL-(j-1)*1:newL*j);
0143     newTime_a = newTime(j+(j-1)*newL-(j-1)*1:newL*j);
0144     new_AzSource_a = new_AzSource(j+(j-1)*newL-(j-1)*1:newL*j);
0145     new_ElSource_a = new_AzSource(j+(j-1)*newL-(j-1)*1:newL*j);
0146     newAz_a = newAz(j+(j-1)*newL-(j-1)*1:newL*j);
0147     newEl_a = newEl(j+(j-1)*newL-(j-1)*1:newL*j);
0148     new_prog_El_offset_a = new_prog_El_offset(j+(j-1)*newL-(j-1)*1:newL*j);
0149     new_prog_Az_offset_a = new_prog_Az_offset(j+(j-1)*newL-(j-1)*1:newL*j);
0150     
0151   % find where we are actually on source according to the telescope
0152     ON = find(new_prog_El_offset_a ==0 & new_prog_El_offset_a==0,1);
0153     on_Az = newAz_a(ON); 
0154     on_El = newEl_a(ON);
0155        
0156  % start gettting Gaussian parameters
0157    % xx1 = find(newD1_a == nanmax(newD1_a));
0158     xx1 = ON;
0159     x_01_x = sky_Az_a(xx1); % central x position
0160     x_01_y = sky_El_a(xx1); % central y position
0161     x_01_x = x_01_x(1);
0162     x_01_y = x_01_y(1);
0163     
0164 L = length(newD1_a); 
0165 
0166 % to find offset gradient by fitting a straight line to the non Gaussian
0167 % part of scan
0168 
0169 off1 = remove_background(newTime_a, newD1_a, 'WindowSize',1,'StepSize',0.2);
0170 straight1 = newD1_a-off1;
0171 ascS1 = sort(straight1);
0172 
0173 %if list_az(1) == 0 && list_az(L) == 0
0174 
0175     a = straight1(1);
0176     b = straight1(L);
0177     e = sky_Az_a(1);
0178     f = sky_Az_a(L);
0179     g = sky_El_a(1);
0180     h = sky_El_a(L);
0181  
0182 fit1 = cat(1, a, b);
0183 az_cut = cat(1,e,f);
0184 el_cut = cat(1,g,h);
0185  
0186 Px1 = polyfit(az_cut, fit1,1);
0187 Py1 = polyfit(el_cut, fit1,1);
0188 grad1X = Px1(1);
0189 grad1Y = Py1(1);
0190 
0191 asc1 = sort(newD1_a);
0192 offset1 = asc1(1);
0193 
0194 A1 = newD1_a(ON) - straight1(ON); 
0195 
0196 Sigma_fit = 0.73/(2*sqrt(2*log(2)));   
0197 % This is cos we scan at 0.5 deg per second for calib scans NO NOT USING
0198 Sigma_x_fit = Sigma_fit;%/0.5;
0199 Sigma_y_fit = Sigma_fit;%/0.5;
0200 
0201 BBeta1 = [A1, x_01_x, Sigma_x_fit, x_01_y, Sigma_y_fit, grad1X, grad1Y, offset1];
0202 %BBeta1 = [A1, x_01_x, x_01_y, grad1X, grad1Y, offset1];
0203 
0204 for I = 1 : length(newD1_a)
0205     if sky_Az_a(I) < 0.5 || sky_Az_a(I) > 1.7 && sky_Az_a(I) < 3.0
0206        within1(I) = I;
0207     else within1(I) = 0;
0208     end
0209 end
0210 
0211 for I = 1 : length(newD1_a)
0212     if sky_Az_a(I) > -0.5 || sky_Az_a(I) < -1.7 && sky_Az_a(I) > -3.0
0213        within2(I) = I;
0214     else within2(I) = 0;
0215     end
0216 end
0217 
0218 for I = 1 : length(newD1_a)
0219     if sky_El_a(I) < 0.5 || sky_El_a(I) > 1.7 && sky_Az_a(I) < 3.0
0220        within3(I) = I;
0221     else within3(I) = 0;
0222     end
0223 end
0224 
0225 for I = 1 : length(newD1_a)
0226     if sky_El_a(I) > -0.5 || sky_El_a(I) < -1.7  && sky_Az_a(I) > -3.0
0227        within4(I) = I;
0228     else within4(I) = 0;
0229     end
0230 end
0231 
0232 if strcmp(name1, 'CygA') == 1
0233    g = find(newD1_a  > mode(newD1_a)+0.03);
0234    gg = last(g);
0235    tracker = sky_El_a(gg);
0236   if tracker > 0   
0237      for I = 1 : length(newD1_a)
0238         if sky_El_a(I) < 0.5
0239             justCyg1(I) = I;
0240         else justCyg1(I) = 0;
0241         end
0242      end
0243   elseif tracker < 0
0244         for I = 1 : length(newD1_a)
0245         if sky_El_a(I) > -0.5
0246             justCyg1(I) = I;
0247         else justCyg1(I) = 0;
0248         end
0249      end
0250   end
0251 end  
0252 
0253 within1 = within1(find(within1));  
0254 within2 = within2(find(within2)); 
0255 within3 = within3(find(within3));  
0256 within4 = within4(find(within4)); 
0257 
0258 tog1 = intersect(within1, within3);
0259 tog2 = intersect(within2, within4);
0260 tog = intersect(tog1, tog2);
0261 tog = unique(tog);
0262 
0263 if strcmp(name1, 'CygA') == 1
0264 cygTog1 = justCyg1(find(justCyg1));
0265 cygTog = intersect(cygTog1, tog);
0266 end
0267 
0268 if strcmp(name1, 'CygA') == 1
0269    sky_Az_a2 = sky_Az_a(cygTog);
0270    sky_El_a2 = sky_El_a(cygTog);
0271    fitData = newD1_a(cygTog);
0272 else
0273    sky_Az_a2 = sky_Az_a(tog);
0274    sky_El_a2 = sky_El_a(tog);
0275    fitData = newD1_a(tog);
0276 end
0277 
0278 lenAz = length(sky_Az_a2);
0279 xy1 = cat(2, sky_Az_a2, sky_El_a2, lenAz);
0280 
0281 [betanew1,resid1, J1, Cov1]=nlinfit(xy1',fitData,@gauss2D_fit,BBeta1);
0282 Z1 = gauss2D_draw(betanew1, sky_Az_a, sky_El_a);
0283 
0284 if popt == 1
0285 figure(j)
0286 h = plot3(sky_Az_a2, sky_El_a2, fitData,'o');
0287 set(h(1), 'MarkerEdgeColor', 'none', 'MarkerFaceColor', 'k')
0288 hold on 
0289 plot3(sky_Az_a, sky_El_a, newD1_a, 'co',sky_Az_a, sky_El_a, Z1,'r')
0290 xlabel('SkyAz')
0291 ylabel('SkyEl')
0292 hold off
0293 figure(j+2)
0294 if strcmp(name1, 'CygA') == 1
0295 plot(newTime_a(cygTog), fitData,'ko')
0296 else plot(newTime_a(tog), fitData,'ko')
0297 end
0298 hold on 
0299 plot(newTime_a, newD1_a, 'c.',newTime_a, Z1,'r');
0300 xlabel('SkyAz')
0301 ylabel('SkyEl')
0302 hold off
0303 end
0304 
0305 Channel1_Amp(j) = betanew1(1);
0306 az_offset1(j) = sky_Az_a2(ON);
0307 el_offset1(j) = sky_El_a2(ON);
0308 El1(j) = newEl_a(xx1);
0309 Az1(j) = newAz_a(xx1);
0310 utcTime(j) = new_time(xx1);
0311 
0312 amp_error1(j) = (Cov1(1,1)/Channel1_Amp(j))*100;
0313 chi1(j) = chi2_beamfits(betanew1, xy1', fitData);
0314 
0315     %find out where we are actually looking when telescope thinks offsets
0316     %are 0
0317     real_Az = on_Az + az_offset1(j); %sky_Az_a(ON);
0318     real_El = on_El + el_offset1(j); %sky_El_a(ON);
0319     
0320 total_off1(j) = spaceangle(real_Az, real_El, on_Az, on_El, 'deg');
0321 powerOff1(j) = exp(-(total_off1(j)^2)/0.1922);
0322 new_Channel1_Amp(j) = Channel1_Amp(j)/powerOff1(j);
0323 
0324 SigmaAz(j) = betanew1(3);
0325 SigmaEl(j) = betanew1(5);
0326 
0327 % paralactic angle
0328 Pang(j) = parangle(real_Az, real_El, 'N');
0329 
0330 end 
0331 
0332 utc(1:j,hhh) = utcTime;
0333 source_num(1:j,hhh) = name_num;
0334 Amp(1:j,hhh) = Channel1_Amp;
0335 newAmp(1:j,hhh) = new_Channel1_Amp;
0336 el_off(1:j,hhh) = el_offset1;
0337 az_off(1:j,hhh) = az_offset1;
0338 tot_off(1:j,hhh) = total_off1;
0339 SigAz(1:j,hhh) = SigmaAz;
0340 SigEl(1:j,hhh) = SigmaEl;
0341 EL(1:j,hhh) = El1;
0342 chi(1:j,hhh) = chi1;
0343 error(1:j,hhh) = amp_error1;
0344 SkyDip(1:j,hhh) = skydip;
0345 Theo(1:j,hhh) = TheoOp;
0346 Pipe(1:j,hhh) = pipeTheo;
0347 pang(1:j,hhh) = Pang;
0348 AZ(1:j,hhh) = Az1;
0349 
0350 end
0351 
0352 if repeat ==2
0353 utc = cat(1,utc(:,1), utc(:,2));
0354 source_num = cat(1,source_num(:,1), source_num(:,2));
0355 Amp = cat(1,Amp(:,1), Amp(:,2));
0356 newAmp = cat(1,newAmp(:,1), newAmp(:,2));
0357 el_off = cat(1,el_off(:,1), el_off(:,2));
0358 az_off = cat(1,az_off(:,1), az_off(:,2));
0359 tot_off = cat(1,tot_off(:,1), tot_off(:,2));
0360 SigAz = cat(1,SigAz(:,1), SigAz(:,2));
0361 SigEl = cat(1,SigEl(:,1), SigEl(:,2));
0362 EL = cat(1,EL(:,1), EL(:,2));
0363 chi = cat(1, chi(:,1), chi(:,2));
0364 error = cat(1,error(:,1), error(:,2));
0365 SkyDip = cat(1,SkyDip(:,1), SkyDip(:,2));
0366 Theo = cat(1,Theo(:,1), Theo(:,2));
0367 Pipe = cat(1,Pipe(:,1), Pipe(:,2));
0368 pang = cat(1,pang(:,1), pang(:,2));
0369 AZ = cat(1,AZ(:,1), AZ(:,2));
0370 end
0371 
0372   end
0373