%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [obs, off, ide] = radioPointPlots(d, plotparams, type, field, maindir) pipeline radio pointing analysis. d is the data structure type is either 'raster' -- schedule did a full raster on the source (az/el) 'scan' -- schedule did an az/el scan 'cross' -- schedule did discrete steps and integrated gen is whether to generate the plots or not. off are the offsets. sjcm modified: 3/23/2011 - sjcm: made it compatible with data from before October 2010. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 function [obs, off, ide, out] = radioPointPlots(d, plotparams, type, field, ... 0002 maindir, useMel) 0003 0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 % 0006 % function [obs, off, ide] = radioPointPlots(d, plotparams, type, field, maindir) 0007 % 0008 % pipeline radio pointing analysis. 0009 % 0010 % d is the data structure 0011 % type is either 0012 % 'raster' -- schedule did a full raster on the source (az/el) 0013 % 'scan' -- schedule did an az/el scan 0014 % 'cross' -- schedule did discrete steps and integrated 0015 % gen is whether to generate the plots or not. 0016 % 0017 % off are the offsets. 0018 % sjcm 0019 % 0020 % modified: 3/23/2011 - sjcm: made it compatible with data from before 0021 % October 2010. 0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0023 0024 if(nargin<5) 0025 maindir = []; 0026 useMel = 0; 0027 end 0028 if(nargin<6) 0029 useMel = 0; 0030 end 0031 0032 0033 0034 % we must determine the az/el offset between the source and its "thought-of" 0035 % location 0036 % find data where the tracker utc is not 0037 if(length(d.antenna0.tracker.utc)~=length(unique(d.antenna0.tracker.utc))) 0038 d.antenna0.tracker.utc = d.array.frame.utc + 1/60/60/24; 0039 end 0040 0041 numCols = size(d.antenna0.receiver.data,2); 0042 0043 if( (mean(d.antenna0.thermal.coldLoad) > 20) & (useMel == 0) ) 0044 % isn't needed with Mel's stuff. 0045 % I channels are negative 0046 I1min = min(d.antenna0.receiver.data(:,1)); 0047 I2min = min(d.antenna0.receiver.data(:,numCols)); 0048 d.antenna0.receiver.data(:,1) = d.antenna0.receiver.data(:,1) - I1min; 0049 d.antenna0.receiver.data(:,numCols) = d.antenna0.receiver.data(:,numCols) - I2min; 0050 end 0051 0052 0053 0054 switch type 0055 case 'raster' 0056 % [obs ide] = calcOffRaster(d, gen, field); 0057 0058 case 'scan' 0059 if(useMel) 0060 [obs ide off out] = calcOffScanMel(d, plotparams, field, maindir); 0061 else 0062 [obs ide off out] = calcOffScan(d, plotparams, field, maindir); 0063 end 0064 0065 0066 0067 case 'cross' 0068 [obs ide off] = calcOffCross(d, plotparams, field); 0069 end 0070 0071 % Put az values into 0-360 range 0072 ind = obs.az<0; obs.az(ind) = obs.az(ind)+360; 0073 ind = obs.az>360; obs.az(ind) = obs.az(ind)-360; 0074 0075 % now we just return the obs, ide, off 0076 0077 0078 0079 return; 0080 0081 0082 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0083 %function [obs ide off] = calcOffScan(d, plotparams, field, maindir) 0084 % 0085 %setPlotDisplay(plotparams.plot); 0086 % 0087 %% load up the current model 0088 %om(:,1)=d.antenna0.tracker.flexure(:,1); 0089 %om(:,2)=d.antenna0.tracker.flexure(:,2); 0090 %om(:,3)=d.antenna0.tracker.tilts(:,1); 0091 %om(:,4)=d.antenna0.tracker.tilts(:,2); 0092 %om(:,5)=d.antenna0.tracker.tilts(:,3); 0093 %om(:,6)=d.antenna0.tracker.fixedCollimation(:,1); 0094 %om(:,7)=d.antenna0.tracker.fixedCollimation(:,2); 0095 %om(:,8)=d.antenna0.tracker.encoder_off(:,1); 0096 %om(:,9)=d.antenna0.tracker.encoder_off(:,2); 0097 %om = mean(om); 0098 % 0099 %% I. split them up into respective scans 0100 %% Ia. find the start and stop index of each scan. 0101 %[s e] = findStartStop(d.index.radio_point_scan.slow); 0102 % 0103 %if(isempty(e) | isempty(s)) 0104 % obs.az = nan; 0105 % obs.el = nan; 0106 % ide.az = nan; 0107 % ide.el = nan; 0108 % off.az = nan; 0109 % off.el = nan; 0110 % goodPoint = [0 0]; 0111 % display('No crosses'); 0112 % s = []; 0113 % return; 0114 % 0115 %end 0116 % 0117 %if(e(1)<s(1)) 0118 % obs.az = nan; 0119 % obs.el = nan; 0120 % ide.az = nan; 0121 % ide.el = nan; 0122 % off.az = nan; 0123 % off.el = nan; 0124 % goodPoint = [0 0]; 0125 % display('No crosses'); 0126 % s = []; 0127 % return; 0128 %end 0129 % 0130 %betaAzAll = nan(length(s), 2, 3); 0131 %betaElAll = nan(length(s), 2, 3); 0132 % 0133 %plotIndex = 0; 0134 %for m=1:length(s) 0135 % % cut the data for each scan 0136 % ind = zeros(size(d.array.frame.features)); 0137 % ind(s(m):e(m)) = 1; 0138 % ind = logical(ind); 0139 % dc = framecut(d, ind); 0140 % 0141 % idevals = dc.antenna0.tracker.horiz_topo(1,:); 0142 % [obsvals(1) obsvals(2)] = pointing_model(om, idevals(1), idevals(2)); 0143 % ide.az(m) = idevals(1); 0144 % ide.el(m) = idevals(2); 0145 % obs.az(m) = obsvals(1); 0146 % obs.el(m) = obsvals(2); 0147 % 0148 % % next we split it into the azimuth/elevation scan 0149 % % for the elevation scan, we fit a line and remove the slope (effectively 0150 % % a sky dip), and fit a gaussian to the remainder 0151 % % for the azimuth scan, we just fit a gaussian to it. 0152 % 0153 % % get the offsets from the control system. 0154 % azApp = interp1(dc.antenna0.tracker.utc, ... 0155 % dc.antenna0.tracker.horiz_topo(:,1), dc.antenna0.receiver.utc); 0156 % azOffSave = azApp - dc.antenna0.servo.apparent(:,1); 0157 % azOffSave = -azOffSave; 0158 % 0159 % elApp = interp1(dc.antenna0.tracker.utc, ... 0160 % dc.antenna0.tracker.horiz_topo(:,2), dc.antenna0.receiver.utc); 0161 % elOffSave = elApp - dc.antenna0.servo.apparent(:,2); 0162 % elOffSave = -elOffSave; 0163 % 0164 % % assuming the slowest we would ever scan is at 0.2 degrees per second, a 0165 % % az scan happens when the az speed is mostly constant, increasing, and 0166 % % greater than 0.2 degrees/s, which is 0167 % if(dc.array.frame.utc(1)<date2mjd(2010,10,1,0,0,0)) 0168 % azInd = (deriv(dc.antenna0.servo.az))<-0.2/100; 0169 % else 0170 % azInd = (deriv(dc.antenna0.servo.az))>0.1/100; 0171 % end 0172 % azIs = dc.antenna0.receiver.data(azInd,[1 6]); 0173 % azIs(:,1) = smooth(azIs(:,1)); 0174 % azIs(:,2) = smooth(azIs(:,2)); 0175 % 0176 % % if amplitude is negative, no peak in data 0177 % 0178 % % same goes for elevation 0179 % if(dc.array.frame.utc(1)<date2mjd(2010,10,1,0,0,0)) 0180 % elInd = deriv(dc.antenna0.servo.el)<-0.2/100; 0181 % else 0182 % elInd = deriv(dc.antenna0.servo.el)>0.1/100; 0183 % end 0184 % elIs = dc.antenna0.receiver.data(elInd,[1 6]); 0185 % elIs(:,1) = smooth(elIs(:,1)); 0186 % elIs(:,2) = smooth(elIs(:,2)); 0187 % 0188 % 0189 % testVals = [120:30:800]; 0190 % for mm=1:2 0191 % 0192 % if(~isempty(find(azInd))) 0193 % % fit a gaussian about the max in both cases 0194 % azOff = azOffSave; 0195 % azI = azIs(:,mm); 0196 % f = find(azI == max(azI)); 0197 % azOff = azOff(azInd); 0198 % 0199 % stopAz = 0; 0200 % 0201 % azIndex = 1; 0202 % while(stopAz==0) 0203 % azFit = find(azInd); 0204 % if(f<=testVals(azIndex) | (f+testVals(azIndex))>=length(azI)) 0205 % badAzPoint = 1; 0206 % azFit = []; 0207 % xoff = nan; 0208 % azIndex = azIndex+1; 0209 % betaAz(1:3) = nan; 0210 % else 0211 % azFit = azOff((f-testVals(azIndex):f+testVals(azIndex))); 0212 % azIFit = azI(f-testVals(azIndex):f+testVals(azIndex)); 0213 % 0214 %% display(sprintf('testing value: %d', testVals(azIndex))); 0215 % 0216 % OPTIONS = optimset('Display', 'off', 'MaxIter', 200000, ... 0217 % 'MaxFunEvals', 20000, 'TolX', 0.00005, 'TolFun', 0.00005); 0218 % LB = [0 azFit(testVals(azIndex))-0.2 -10 -inf -inf]; 0219 % UB = [max(azIFit) azFit(testVals(azIndex))+0.2 10 inf inf]; 0220 % x0 = [1 azFit(testVals(azIndex)-1) 1 1 1]; 0221 % [betaAz resnorm residual exitflag] = lsqnonlin(@(x) mchisq2(azIFit, gaussfit(x, azFit)), x0, LB, UB, OPTIONS); 0222 % 0223 % amplitude = betaAz(1); 0224 % xoff = betaAz(2); 0225 % sigma = betaAz(3); 0226 % 0227 % % other terms remove a baseline 0228 % if (amplitude<0 | abs(xoff)>1 | resnorm>500) 0229 % badAzPoint = 1; 0230 % else 0231 % badAzPoint = 0; 0232 % end 0233 % azFit = gaussfit(betaAz, azOff); 0234 % 0235 % if(abs(xoff)<2) 0236 % stopAz = 1; 0237 % else 0238 % azIndex = azIndex+1; 0239 % end 0240 % end 0241 % 0242 % if(azIndex>length(testVals)) 0243 % stopAz = 1; 0244 % end 0245 % end 0246 % 0247 % betaAzAll(m,mm,:) = betaAz(1:3); 0248 % end 0249 % 0250 % 0251 % if(~isempty(elInd)) 0252 % elOff = elOffSave; 0253 % elI = elIs(:,mm); 0254 % % first we remove the offset from slewing the atmosphere 0255 % % from the sky dip code 0256 % x = 1./(sind(dc.antenna0.servo.el(elInd))); 0257 % y = abs(elI); 0258 % [tau Tground] = linfit(x,y); 0259 % % remove the effect 0260 % elI = elI - Tground(1) - tau(1)./sind(dc.antenna0.servo.el(elInd)); 0261 % 0262 % % fit a gaussian about the max 0263 % f = find(elI == max(elI)); 0264 % elOff2 = elOff(elInd); 0265 % 0266 % stopEl = 0; 0267 % elIndex = 1; 0268 % while(stopEl==0) 0269 % elFit = find(elInd); 0270 % if(f<=testVals(elIndex) | (f+testVals(elIndex))>=length(elI)) 0271 % badElPoint = 1; 0272 % elFit = []; 0273 % yoff = nan; 0274 % elIndex = elIndex+1; 0275 % betaEl = [nan nan nan]; 0276 % else 0277 % elFit = elOff2((f-testVals(elIndex):f+testVals(elIndex))); 0278 % elIFit = elI(f-testVals(elIndex):f+testVals(elIndex)); 0279 % 0280 %% display(sprintf('el testing value: %d', testVals(elIndex))); 0281 % [betaEl] = nlinfit(elFit, elIFit, @gaussfit, [1 elFit(testVals(elIndex)) 1 1 1]); 0282 % amplitude = betaEl(1); 0283 % yoff = betaEl(2); 0284 % sigma = betaEl(3); 0285 % % other terms remove a baseline 0286 % if (amplitude<0 | abs(yoff)>1) 0287 % badElPoint = 1; 0288 % else 0289 % badElPoint = 0; 0290 % end 0291 % elFit = gaussfit(betaEl, elOff2); 0292 % 0293 % if(abs(yoff)<1) 0294 % stopEl = 1; 0295 % else 0296 % elIndex = elIndex+1; 0297 % end 0298 % end 0299 % if(elIndex>length(testVals)) 0300 % stopEl = 1; 0301 % end 0302 % 0303 % end 0304 % betaElAll(m,mm,:) = betaEl(1:3); 0305 % end 0306 % 0307 %% display('here'); 0308 % badPoint(m,mm) = badAzPoint | badElPoint; 0309 % 0310 % if(~badPoint(m,mm)) 0311 % % then we plot it 0312 % 0313 % if(~any(isnan(betaAzAll(m,mm,:)))) 0314 % 0315 % subplot(2,1,1) 0316 % plot((azOff)*cos(dc.antenna0.servo.el(1)*pi/180), azI'); 0317 % xlabel('Az Offset (degrees)'); 0318 % ylabel('Intensity'); 0319 % if(~isempty(azFit)) 0320 % hold on; plot(azOff*cos(dc.antenna0.servo.el(1)*pi/180), azFit, 'r'); hold off 0321 % axis([min(azOff) max(azOff) min(azI) max(azI)]); 0322 % end 0323 % 0324 % end 0325 % 0326 % if(~any(isnan(betaElAll(m,mm,:)))) 0327 % subplot(2,1,2) 0328 % plot(elOff(elInd), elI); 0329 % xlabel('El Offset (degrees)'); 0330 % ylabel('Intensity'); 0331 % if(~isempty(elFit)) 0332 % hold on; plot(elOff(elInd), elFit, 'r'); hold off 0333 % axis([min(elOff) max(elOff) min(elI) max(elI)]); 0334 % end 0335 % end 0336 % eval(sprintf('gtitle(''Pointing scan on %s - Chan %d'');', ... 0337 % dc.antenna0.tracker.source{1}, mm)); 0338 % 0339 % display(sprintf('Source Name: %s', dc.antenna0.tracker.source{1})); 0340 % 0341 % if(plotparams.save==1) 0342 % plotIndex = plotIndex+1; 0343 % if(isempty(maindir)) 0344 % maindir = getMainDir(d, field); 0345 % end 0346 % dbclear if error 0347 % set(gcf,'paperposition',[0 0 6.0 9.0]) 0348 % filename = sprintf('%s/pointing/fig%d.png', maindir, plotIndex); 0349 % eval(sprintf('print -dpng -r70 %s', filename)); 0350 % dbstop if error 0351 % end 0352 % end 0353 % end 0354 % 0355 % clear azIs; 0356 % clear elIs; 0357 % pause(0.1); 0358 % sources{m} = unique(dc.antenna0.tracker.source); 0359 % timeVal(m) = dc.antenna0.receiver.utc(1); 0360 %end 0361 % 0362 %badPoint = logical(badPoint); 0363 %elOff = betaElAll(:,:,2); 0364 %azOff = betaAzAll(:,:,2); 0365 %elOff(badPoint) = nan; 0366 %azOff(badPoint) = nan; 0367 % 0368 % 0369 %% if hte values for both channels don't agree, toss them. 0370 %dEl = abs(diff(elOff, [], 2)); 0371 %dAz = abs(diff(azOff, [], 2).*cos(obs.el*pi/180)'); 0372 % 0373 %% if they differ by more than 20 arcminutes, throw them out 0374 %f = find(dEl>20/60 | dAz>20/60); 0375 %elOff(f,:) = nan; 0376 %azOff(f,:) = nan; 0377 % 0378 %off.az = nanmean(azOff, 2); 0379 %off.el = nanmean(elOff, 2); 0380 %off.azErr = nanstd(azOff, [], 2); 0381 %off.elErr = nanstd(elOff, [], 2); 0382 %% if there's only one sample, we need to set a weight equal to the other 0383 %% ones. 0384 %indSingle = off.azErr == 0; 0385 %off.azErr(indSingle) = median(off.azErr(~indSingle)); 0386 %off.elErr(indSingle) = median(off.elErr(~indSingle)); 0387 % 0388 %display('Done with all crosses'); 0389 % 0390 %% add final offests to the observed az/el 0391 %obs.az = obs.az + off.az'; 0392 %obs.el = obs.el + off.el'; 0393 % 0394 % 0395 %% throw out the data that is bad 0396 %ind = isnan(obs.az); 0397 %obs.az(ind) = []; 0398 %obs.el(ind) = []; 0399 %ide.az(ind) = []; 0400 %ide.el(ind) = []; 0401 %off.az(ind) = []; 0402 %off.el(ind) = []; 0403 %off.azErr(ind) = []; 0404 %off.elErr(ind) = []; 0405 %off.sources = sources(~ind); 0406 %off.timeVal = timeVal(~ind); 0407 % 0408 % 0409 %return; 0410 % 0411 0412 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0413 function [obs ide] = calcOffRaster(d, gen) 0414 0415 % I. split them up into respective scans 0416 % Ia. find the start and stop index of each scan. 0417 [s e] = findStartStop(d.array.frame.features); 0418 0419 0420 for m=1:length(s) 0421 % cut the data for each scan 0422 ind = zeros(size(d.array.frame.features)); 0423 ind(s(m):e(m)) = 1; 0424 ind = logical(ind); 0425 dc = framecut(d, ind); 0426 0427 % calculate the offset from nominal for this scan 0428 [obs.az(m) obs.el(m)] = sourceScanMap(dc, 0.5, 0); 0429 0430 % calculate the az/el for the source 0431 ide.az(m) = mean(dc.antenna0.tracker.horiz_mount(:,1)); 0432 ide.el(m) = mean(dc.antenna0.tracker.horiz_mount(:,2)); 0433 end 0434 0435 0436 return; 0437 0438 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0439 %function [s e] = findStartStop(features) 0440 % 0441 %% start indices are when the difference in the feature value increases, stop 0442 %% is when it decreases 0443 %s = find(diff(double(features))>0); 0444 % 0445 %% end is when it decreases 0446 %e = find(diff(double(features))<0); 0447 % 0448 %if(length(e)<length(s)) 0449 % % the dimensions don't match 0450 % if(s(1)<e(1)) 0451 % % if start is less than end, that means we're usually missing hte last 0452 % % endpoint 0453 % e(length(s)) = length(features); 0454 % else 0455 % s(1) = []; 0456 % end 0457 %elseif(length(s)<length(e)) 0458 % s(1) = []; 0459 %end 0460 %s = s+1; 0461 % 0462 %goodPoint = ones(size(s)); 0463 %for m=1:length(s) 0464 % thisFeat = unique(features(s(m):e(m))); 0465 % if(thisFeat ~= 1) 0466 % goodPoint(m) = 0; 0467 % end 0468 %end 0469 % 0470 %s = s(logical(goodPoint)); 0471 %e = e(logical(goodPoint)); 0472 % 0473 %return; 0474 0475 0476 function x = mchisq2(dataVals, modelVals) 0477 0478 x = dataVals - modelVals; 0479 0480 f = find(isnan(x)); 0481 x(f) = []; 0482 f = find(isinf(x)); 0483 x(f) = []; 0484 0485 if(isempty(x)) 0486 display('damn it all to hell') 0487 keyboard; 0488 end 0489 0490 return; 0491 0492 0493 0494 0495 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0496 function [obs ide off] = calcOffCross(d, plotparams, field); 0497 0498 setPlotDisplay(plotparams.plot); 0499 0500 % load up the current model 0501 om(:,1)=d.antenna0.tracker.flexure(:,1); 0502 om(:,2)=d.antenna0.tracker.flexure(:,2); 0503 om(:,3)=d.antenna0.tracker.tilts(:,1); 0504 om(:,4)=d.antenna0.tracker.tilts(:,2); 0505 om(:,5)=d.antenna0.tracker.tilts(:,3); 0506 om(:,6)=d.antenna0.tracker.fixedCollimation(:,1); 0507 om(:,7)=d.antenna0.tracker.fixedCollimation(:,2); 0508 om(:,8)=d.antenna0.tracker.encoder_off(:,1); 0509 om(:,9)=d.antenna0.tracker.encoder_off(:,2); 0510 om = mean(om); 0511 0512 % calculates the offsets when doing a pointing cross: fits for a gaussian 0513 % in x and y offsets, converts to offsets in az/el, and calculates the 0514 % observed location of the source 0515 % should be very similar to the sza 0516 0517 [s e] = findCrossStartStop(d); 0518 0519 % now for each of these fields, we split them up into the x and y offsets, 0520 % fit a gaussian, and find the offsets. convert offsets to az/el offsets, 0521 % and be done. 0522 for m=1:length(s) 0523 % cut the data for each scan 0524 ind = zeros(size(d.array.frame.features)); 0525 ind(s(m):e(m)) = 1; 0526 ind = logical(ind); 0527 dc = framecut(d, ind); 0528 0529 [intStart intStop] = findCrossStepStartStop(dc); 0530 if(length(intStart)~=26) 0531 display('You may have some mixed up data. not the right number of points'); 0532 obs = []; 0533 ide = []; 0534 off = []; 0535 return; 0536 end 0537 0538 for n=1:length(intStart) 0539 ind = zeros(size(dc.array.frame.features)); 0540 ind(intStart(n):intStop(n)) = 1; 0541 ind = logical(ind); 0542 dcc = framecut(dc, ind); 0543 0544 offsets(n,:) = median(dcc.antenna0.tracker.sky_xy_off); 0545 meanIs(n,:) = mean(dcc.antenna0.receiver.data(:,[1 6])); 0546 stdIs(n,:) = std(dcc.antenna0.receiver.data(:,[1 6])); 0547 0548 ideal(n,:) = dcc.antenna0.tracker.horiz_topo(2,:); 0549 observed(n,:) = dcc.antenna0.tracker.horiz_mount(2,:) + ... 0550 dcc.antenna0.tracker.encoder_off(2,:); 0551 end 0552 fy = find(offsets(:,2)~=0); 0553 yIndex = fy(1):last(fy); 0554 fx = find(offsets(:,1)~=0); 0555 xIndex = fx(1):last(fx); 0556 0557 % fit a gaussian to both data sets. 0558 xOffVals = offsets(xIndex,1); 0559 yOffVals = offsets(yIndex,2); 0560 xIntVals = abs(meanIs(xIndex,:)); 0561 yIntVals = abs(meanIs(yIndex,:)); 0562 xErrVals = stdIs(xIndex,:); 0563 yErrVals = stdIs(yIndex,:); 0564 0565 % find the x offset 0566 x0 = [max(xIntVals(:,1)) 0 0.2]; 0567 ind = abs(xOffVals)<=0.5; 0568 [px(1,:) pex(1,:)] = matmin('chisq', x0, [], 'gauss', ... 0569 xIntVals(ind,1), xErrVals(ind,1), xOffVals(ind)); 0570 [px(2,:) pex(2,:)] = matmin('chisq', x0, [], 'gauss', ... 0571 xIntVals(ind,2), xErrVals(ind,2), xOffVals(ind)); 0572 fullX = min(xOffVals):0.01:max(xOffVals); 0573 aa(:,1) = gauss(px(1,:), fullX); 0574 aa(:,2) = gauss(px(2,:), fullX); 0575 offset.x(m) = vwstat(px(:,2), pex(:,2).^2); 0576 0577 elFit = xOffVals(ind); 0578 elIFit = xIntVals(ind,2); 0579 0580 [betaEl] = nlinfit(elFit, elIFit, @gaussfit, [1 0 1 1 1]); 0581 amplitude = betaEl(1); 0582 yoff = betaEl(2); 0583 sigma = betaEl(3); 0584 xoffvals = min(xOffVals):0.01:max(xOffVals); 0585 elFitFin = gaussfit(betaEl, xoffvals); 0586 0587 y0 = [max(yIntVals(:,1)) 0 0.2]; 0588 ind = abs(yOffVals)<=0.5; 0589 [py(1,:) pey(1,:)] = matmin('chisq', y0, [], 'gauss', ... 0590 yIntVals(ind,1), yErrVals(ind,1), yOffVals(ind)); 0591 [py(2,:) pey(2,:)] = matmin('chisq', y0, [], 'gauss', ... 0592 yIntVals(ind,2), yErrVals(ind,2), yOffVals(ind)); 0593 fullY = min(yOffVals):0.01:max(yOffVals); 0594 bb(:,1) = gauss(py(1,:), fullY); 0595 bb(:,2) = gauss(py(2,:), fullY); 0596 offset.y(m) = vwstat(py(:,2), pey(:,2).^2); 0597 0598 if(any(abs(px(:,2)) > max(xOffVals))) 0599 display('Suspect cross'); 0600 badXPoint = 1; 0601 else 0602 badXPoint = 0; 0603 end 0604 if(abs(diff(px(:,2)))>0.05) 0605 display('Best Fit X offset does not agree between channels'); 0606 badXPoint = 1; 0607 end 0608 0609 0610 if(any(abs(py(:,2)) > max(xOffVals))) 0611 display('Suspect cross'); 0612 badYPoint = 1; 0613 else 0614 badYPoint = 0; 0615 end 0616 if(abs(diff(py(:,2)))>0.05) 0617 display('Best Fit Y offset does not agree between channels'); 0618 badYPoint = 1; 0619 end 0620 0621 0622 0623 % plot the results to see the fit. 0624 meanX = mean(xIntVals); 0625 xIntVals(:,1) = xIntVals(:,1) - meanX(1); 0626 xIntVals(:,2) = xIntVals(:,2) - meanX(2); 0627 aa(:,1) = aa(:,1) - meanX(1); 0628 aa(:,2) = aa(:,2) - meanX(2); 0629 0630 meanY = mean(yIntVals); 0631 yIntVals(:,1) = yIntVals(:,1) - meanY(1); 0632 yIntVals(:,2) = yIntVals(:,2) - meanY(2); 0633 bb(:,1) = bb(:,1) - meanY(1); 0634 bb(:,2) = bb(:,2) - meanY(2); 0635 0636 subplot(2,1,1) 0637 plot(xOffVals, xIntVals, '*-'); 0638 hold on 0639 plot(fullX, aa, 'r') 0640 axis([min(xOffVals) max(xOffVals) min(xIntVals(:)) max(xIntVals(:))]); 0641 legend('Channel 1', 'Channel 2', 'Fit Gaussian'); 0642 xlabel('X offset'); 0643 ylabel('Counts'); 0644 hold off; 0645 0646 subplot(2,1,2) 0647 plot(yOffVals, yIntVals, '*-'); 0648 hold on 0649 plot(fullY, bb, 'r') 0650 axis([min(yOffVals) max(yOffVals) min(yIntVals(:)) max(yIntVals(:))]); 0651 legend('Channel 1', 'Channel 2', 'Fit Gaussian'); 0652 xlabel('Y offset'); 0653 ylabel('Counts'); 0654 hold off; 0655 0656 if(plotflag.save==1) 0657 plotIndex = plotIndex+1; 0658 maindir = getMainDir(d, field); 0659 dbclear if error 0660 set(gcf,'paperposition',[0 0 6.0 9.0]) 0661 filename = sprintf('%s/pointing/fig%d.png', maindir, plotIndex); 0662 eval(sprintf('print -dpng -r70 %s', filename)); 0663 dbstop if error 0664 end 0665 0666 display(sprintf('Source Name: %s', dc.antenna0.tracker.source{1})); 0667 eval(sprintf('display(''X offset: %3.4f, Y offset: %3.4f'');', ... 0668 offset.x(m), offset.y(m))); 0669 if(badXPoint | badYPoint) 0670 badPoint = 1; 0671 display('Cross Suspect'); 0672 end 0673 0674 0675 % check if we want to keep this source 0676 % r = input('Keep This point? [y/n]: ', 's'); 0677 % 0678 % if(strcmp(r, 'y') || strcmp(r, 'Y')) 0679 % badPoint = 0; 0680 % elseif(strcmp(r, 'n') || strcmp(r, 'N')) 0681 % badPoint = 1; 0682 % elseif(strcmp(r, 'k')) 0683 % keyboard; 0684 % end 0685 0686 if(~badPoint) 0687 % now for the ideal and observed 0688 f = find(offsets(:,1)==0 & offsets(:,2)==0); 0689 0690 ide.az(m) = ideal(f(1),1); 0691 ide.el(m) = ideal(f(1),2); 0692 obs.az(m) = observed(f(1),1); 0693 obs.el(m) = observed(f(1),2); 0694 0695 % convert great circle offsets into az/el offsets 0696 0697 a.az = obs.az(m); 0698 a.el = obs.el(m); 0699 b = xyoff_to_azel(a, offset.x(m), offset.y(m)); 0700 off.x(m) = offset.x(m); 0701 off.y(m) = offset.y(m); 0702 obs.az(m) = b.az; 0703 obs.el(m) = b.el; 0704 0705 else 0706 ide.az(m) = nan; 0707 ide.el(m) = nan; 0708 obs.az(m) = nan; 0709 obs.el(m) = nan; 0710 off.x(m) = offset.x(m); 0711 off.y(m) = offset.y(m); 0712 end 0713 end 0714 0715 0716 return; 0717 0718 function [s e] = findCrossStartStop(d) 0719 % check for when the radio pointing flag switches. 0720 0721 index = d.index.radio_point_cross.slow; 0722 0723 s = find(diff(index)>0)+1; 0724 e = find(diff(index)<0); 0725 0726 % that's it. 0727 if(length(e)<length(s)) 0728 % the dimensions don't match 0729 if(s(1)<e(1)) 0730 % if start is less than end, that means we're usually missing hte last 0731 % endpoint 0732 % e(length(s)) = length(d.array.frame.features); 0733 % cut that last start out. 0734 s(length(s)) = []; 0735 else 0736 s(1) = []; 0737 end 0738 elseif(length(s)<length(e)) 0739 s(1) = []; 0740 end 0741 0742 return; 0743 0744 0745 function [s e] = findCrossStepStartStop(d) 0746 % check for when the radio pointing flag switches. 0747 0748 index = d.index.source.slow; 0749 0750 s = find(diff(index)>0)+1; 0751 e = find(diff(index)<0); 0752 0753 % that's it. 0754 if(length(e)<length(s)) 0755 % the dimensions don't match 0756 if(s(1)<e(1)) 0757 % if start is less than end, that means we're usually missing hte last 0758 % endpoint 0759 e(length(s)) = length(d.array.frame.features); 0760 else 0761 s(1) = []; 0762 end 0763 elseif(length(s)<length(e)) 0764 s(1) = []; 0765 end 0766 0767 return; 0768 0769 0770 %%%%%%%%%%%%%%% 0771 function azel=xyoff_to_azel(azel,x,y) 0772 % Apply xy sky offset to az,el coord 0773 % following Martin Shepherd in tracker.c 0774 % Convert all to radians 0775 d2r=pi/180; 0776 azel.az=azel.az*d2r; azel.el=azel.el*d2r; 0777 x=x*d2r; y=y*d2r; 0778 0779 % Convert x,y to r,theta on surface of sphere 0780 t=atan2(sin(x),tan(y)); 0781 r=acos(cos(x).*cos(y)); 0782 0783 % Apply offset 0784 e=azel.el; 0785 top=sin(t).*sin(r); 0786 bot=cos(e).*cos(r)-sin(e).*sin(r).*cos(t); 0787 azel.az=azel.az+atan2(top,bot); 0788 azel.el=asin(sin(e).*cos(r)+cos(e).*sin(r).*cos(t)); 0789 0790 % Back to deg 0791 azel.az=azel.az/d2r; azel.el=azel.el/d2r; 0792 0793 return