radioPointing

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

  function par = radioPointing(d, type, filename)

  radio pointing analysis schedule.

   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
   filename -- string that will be used to name your plots.

  sjcm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function par = radioPointing(d, type, filename)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function par = radioPointing(d, type, filename)
0006 %
0007 %  radio pointing analysis schedule.
0008 %
0009 %   d is the data structure
0010 %   type is either
0011 %   'raster' -- schedule did a full raster on the source (az/el)
0012 %   'scan'  -- schedule did an az/el scan
0013 %   'cross'   -- schedule did discrete steps and integrated
0014 %   filename -- string that will be used to name your plots.
0015 %
0016 %  sjcm
0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0018 
0019 % we must determine the az/el offset between the source and its "thought-of"
0020 % location
0021 
0022 % find data where the tracker utc is not
0023 if(length(d.antenna0.tracker.utc)~=length(unique(d.antenna0.tracker.utc)))
0024   d.antenna0.tracker.utc = d.array.frame.utc + 1/60/60/24;
0025 end
0026 
0027 
0028 if(d.array.frame.utc(1) < date2mjd(2010, 6, 15))
0029   d.antenna0.receiver.data = abs(d.antenna0.receiver.data);
0030 end
0031 
0032 
0033 switch type
0034   case 'raster'
0035     [obs ide] = calcOffRaster(d);
0036  
0037   case 'scan'
0038       
0039     [obs ide off] = calcOffScan(d);
0040 %   [obs1 ide1 off1] = calcOffScanMel(d);
0041     
0042   case 'cross'
0043     [obs ide off] = calcOffCross(d);    
0044 end
0045 
0046 % Put az values into 0-360 range
0047 ind = obs.az<0;   obs.az(ind) = obs.az(ind)+360;
0048 ind = obs.az>360; obs.az(ind) = obs.az(ind)-360;
0049 
0050 sources = unique(d.antenna0.tracker.source);
0051 disp(sprintf('Read in %d data points',length(obs.az)));
0052 disp(sprintf('  ... %d unique sources', length(sources)));
0053     
0054 % Fit the model
0055 % we only want to fit the collimation and flexure
0056 om(:,1)=d.antenna0.tracker.flexure(:,1);
0057 om(:,2)=d.antenna0.tracker.flexure(:,2);
0058 om(:,3)=d.antenna0.tracker.tilts(:,1);
0059 om(:,4)=d.antenna0.tracker.tilts(:,2);
0060 om(:,5)=d.antenna0.tracker.tilts(:,3);
0061 om(:,6)=d.antenna0.tracker.fixedCollimation(:,1);
0062 om(:,7)=d.antenna0.tracker.fixedCollimation(:,2);
0063 om(:,8)=d.antenna0.tracker.encoder_off(:,1);
0064 om(:,9)=d.antenna0.tracker.encoder_off(:,2);
0065 om = mean(om);
0066 
0067 off.azErr(isnan(off.azErr)) = max(off.azErr);
0068 off.elErr(isnan(off.elErr)) = max(off.elErr);
0069 
0070 om = om(1,:);
0071 
0072 figure(1)
0073 [rmssao ind sa0] = plot_res(om, ide, obs, 'No Fit', 0);
0074 sa0 = mean(sa0);
0075 
0076 
0077 if(rmssao*60 < 6)
0078   display('**********ATTENTION*************');
0079   display('The RMS of the current pointing is very low');
0080   display('The model might be good with no update');
0081 else
0082   display('**********ATTENTION*************');
0083   display('The RMS of the current pointing is not very low');
0084   display('Will try to fit for the new solution');
0085 end
0086 
0087 old_m = om;
0088 %opt_m = [ -0.0500 -0.0438 -0.0243 -0.0291 -0.1331 0.1142 0 -0.6572 -0.1666];
0089 %opt_m = [ -0.0479 -0.0377 -0.0303 -0.0302 -0.0716 0.1943 0 -0.7072  -0.1594];
0090 %opt_m = [-0.0905   -0.0627   -0.0238   -0.0292   -0.0736    0.1910 0   -0.7956   -0.2366];
0091 
0092 %om = opt_m;
0093 
0094 
0095 figure(2)
0096 % only fit for collimation
0097 fitPars = ones(size(om));
0098 fitPars([1:5 8:9]) = 0;
0099 [par sa, gof] = do_fit(om, fitPars, obs, ide, off.azErr', off.elErr');
0100 plot_res(par, ide, obs, 'Collimation Fit Only', 0)
0101 
0102 figure(3)
0103 % only fit for flexure
0104 fitPars = ones(size(om));
0105 fitPars(3:9) = 0;
0106 [par sa, gof] = do_fit(om, fitPars, obs, ide, off.azErr', off.elErr');
0107 plot_res(par, ide, obs, 'Flexure Fit Only', 0)
0108 
0109 figure(4)
0110 % fit both collimation and flexure
0111 fitPars = ones(size(om));
0112 fitPars([3:5 8:9]) = 0;
0113 [par sa, gof] = do_fit(om, fitPars, obs, ide, off.azErr', off.elErr');
0114 [par sa, gof] = do_fit(om, fitPars, obs, ide, ones(size(off.azErr')), ones(size(off.elErr')));
0115 [rmssaf ind saf] = plot_res(par, ide, obs, 'Collimation and Flexure', 0);
0116 %saf = mean(saf);
0117 
0118 par
0119 figure(5)
0120 [rmssaf ind saf] = plot_res(par, ide, obs, 'Feb Full Solution', 0);
0121 
0122 fm = par;
0123 
0124 if(rmssaf*60 > 10)
0125   display('**********WARNING*************');
0126   display('Fit to data not so great');
0127   display('Suggest scheduling an optical pointing run');
0128 end
0129 
0130 if(rmssaf > rmssao)
0131   display('**********WARNING*************');
0132   display('Old model better than new fit');
0133 end
0134 
0135 
0136 if(any(abs(par)>7))
0137   display('**********WARNING*************');
0138   display('Some pointing model terms fit are quite large');
0139   display('Be suspect of the resulting model');
0140   display('i.e., CHECK the solution before committing to the control system');
0141 end
0142 
0143 
0144 if(any(abs(om-fm) > 1.5))
0145   display('**********WARNING*************');
0146   display('Some pointing model terms changed by over 1.5 degree');
0147   display('Be suspect of the resulting model');
0148   display('i.e., CHECK the solution before committing to the control system');  
0149 end
0150 
0151 
0152 
0153 display('Pointing solution derived');
0154 display('Add these lines to bottom of pointing.init');
0155 display(sprintf('encoder_zeros %3.4f, %3.4f, 0', fm(8:9)));
0156 display(sprintf('tilts %3.4f, %3.4f, %3.4f', fm(3:5)));
0157 display(sprintf('collimate_fixed radio, %3.4f, %3.4f', ...
0158     fm(6:7)));
0159 display(sprintf('flexure radio, %3.4f, %3.4f', fm(1:2)));
0160 
0161 
0162 
0163   
0164 % that does it.
0165 return;
0166 
0167 
0168 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0169 function [obs ide off] = calcOffScan2(d)
0170 
0171 % load up the current model
0172 om(:,1)=d.antenna0.tracker.flexure(:,1);
0173 om(:,2)=d.antenna0.tracker.flexure(:,2);
0174 om(:,3)=d.antenna0.tracker.tilts(:,1);
0175 om(:,4)=d.antenna0.tracker.tilts(:,2);
0176 om(:,5)=d.antenna0.tracker.tilts(:,3);
0177 om(:,6)=d.antenna0.tracker.fixedCollimation(:,1);
0178 om(:,7)=d.antenna0.tracker.fixedCollimation(:,2);
0179 om(:,8)=d.antenna0.tracker.encoder_off(:,1);
0180 om(:,9)=d.antenna0.tracker.encoder_off(:,2);
0181 om = mean(om);
0182 
0183 % I. split them up into respective scans
0184 % Ia.  find the start and stop index of each scan.
0185 [s e] = findStartStop(d.array.frame.features);
0186 betaAzAll = nan(length(s), 2, 3);
0187 betaElAll = nan(length(s), 2, 3);
0188   
0189 for m=1:length(s)
0190   % cut the data for each scan
0191   ind = zeros(size(d.array.frame.features));
0192   ind(s(m):e(m)) = 1;
0193   ind = logical(ind);
0194   dc  = framecut(d, ind);
0195 
0196   % find the observation az/el
0197 %  indOnSource = dc.antenna0.tracker.offSource==0;
0198 %  indAzOff    = dc.antenna0.tracker.horiz_off(:,1) == 0 & ...
0199 %      dc.antenna0.tracker.scan_off(:,1) ==0;
0200 %  indElOff    = dc.antenna0.tracker.horiz_off(:,2) == 0 & ...
0201 %      dc.antenna0.tracker.scan_off(:,2) ==0;
0202 %  f = find(indOnSource & indAzOff & indElOff);
0203 %  if(length(f)==0)
0204 %    f = find(indOnSource);
0205 %  end
0206 %  if(length(f)>1)
0207 %    obsvals = (dc.antenna0.tracker.horiz_mount(f(1),:)) + dc.antenna0.tracker.encoder_off(f(1),:);
0208 %    idevals = (dc.antenna0.tracker.horiz_topo(f(1),:));
0209 %
0210 %  elseif(length(f) == 1)
0211 %    obsvals = dc.antenna0.tracker.horiz_mount(f,:) + dc.antenna0.tracker.encoder_off(f,:);
0212 %    idevals = (dc.antenna0.tracker.horiz_topo(f,:));
0213 %  end
0214 
0215 
0216   idevals = dc.antenna0.tracker.horiz_topo(1,:);
0217   [obsvals(1) obsvals(2)] = pointing_model(om, idevals(1), idevals(2));
0218   ide.az(m) = idevals(1);
0219   ide.el(m) = idevals(2);
0220   obs.az(m) = obsvals(1);
0221   obs.el(m) = obsvals(2);
0222   
0223   
0224   
0225   
0226   
0227   % next we split it into the azimuth/elevation scan
0228   % for the elevation scan, we fit a line and remove the slope (effectively
0229   % a sky dip), and fit a gaussian to the remainder
0230   % for the azimuth scan, we just fit a gaussian to it.
0231 
0232   % get the offsets from the control system.
0233   azApp = interp1(dc.antenna0.tracker.utc, ...
0234       dc.antenna0.tracker.horiz_topo(:,1), dc.antenna0.receiver.utc);
0235   azOffSave = azApp - dc.antenna0.servo.apparent(:,1);
0236 %  [aa, ee ] =calcAzElCs(dc, 1);
0237 %  ao = aa - dc.antenna0.servo.az;
0238 %  azOffSave = ao;
0239 %  azOffSave(azOffSave<-180) = azOffSave(azOffSave<-180) + 360;
0240 %  azOffSave(azOffSave>180) = azOffSave(azOffSave>180) - 360;
0241    azOffSave = -azOffSave;
0242 
0243   elApp = interp1(dc.antenna0.tracker.utc, ...
0244       dc.antenna0.tracker.horiz_topo(:,2), dc.antenna0.receiver.utc);
0245   elOffSave = elApp - dc.antenna0.servo.apparent(:,2);
0246   %elOffSave = ee - dc.antenna0.servo.el;
0247   elOffSave = -elOffSave;
0248   
0249   % assuming the slowest we would ever scan is at 0.2 degrees per second, a
0250   % az scan happens when the az speed is mostly constant, increasing, and
0251   % greater than 0.2 degrees/s, which is
0252   azInd = (deriv(dc.antenna0.servo.az))>0.1/100;
0253   if(isempty(find(azInd)))
0254     % fudge it so it doesn't crash
0255     azInd(1:2) = 1;
0256     azInd = logical(azInd);
0257     azIs  = dc.antenna0.receiver.data(azInd,[1 6]);
0258     azIs(:) = nan;
0259     azFlags = ones(size(azIs));
0260   else
0261     [startaz endaz] = findStartStop(azInd);
0262     f = find( (endaz - startaz) == max(endaz - startaz));
0263     azInd = zeros(size(azInd));
0264     azInd(startaz(f):endaz(f)) = 1;
0265     azInd = logical(azInd);
0266     azIs  = dc.antenna0.receiver.data(azInd,[1 6]);
0267     azIs(:,1) = smooth(azIs(:,1));
0268     azIs(:,2) = smooth(azIs(:,2));
0269     
0270     azFlags       = dc.flags.fast(azInd,[1 3]);
0271     azIs(azFlags) = nan;
0272   end
0273   
0274   % if amplitude is negative, no peak in data
0275   
0276   % same goes for elevation
0277   elInd = deriv(dc.antenna0.servo.el)>0.1/100;
0278   if(isempty(find(elInd)))
0279     elInd(1:2) = 1;
0280     elInd = logical(elInd);
0281     elIs  = dc.antenna0.receiver.data(elInd,[1 6]);
0282     elIs(:) = nan;
0283     elFlags = ones(size(elIs));    
0284   else
0285     [startel endel] = findStartStop(elInd);
0286     f = find( (endel-startel) == max(endel - startel));
0287     elInd = zeros(size(elInd));
0288     elInd(startel(f):endel(f)) = 1;
0289     elInd = logical(elInd);
0290     elIs  = dc.antenna0.receiver.data(elInd,[1 6]);
0291     elIs(:,1) = smooth(elIs(:,1));
0292     elIs(:,2) = smooth(elIs(:,2));
0293     
0294     elFlags       = dc.flags.fast(elInd, [1 3]);
0295     elIs(elFlags) = nan;
0296   end
0297     
0298   % if most of the scan is flagged, don't bother fitting.
0299   azFlagPer = length(find(azFlags))./length(azIs(:));
0300   elFlagPer = length(find(elFlags))./length(elIs(:));
0301   if(azFlagPer < 0.15 & elFlagPer < 0.15)
0302     doFit = 1;
0303   else
0304     doFit = 0;
0305   end
0306 
0307 
0308   if(doFit)
0309     testVals = [120:30:800];
0310     for mm=1:2
0311       
0312       if(~isempty(find(azInd)))
0313     % fit a gaussian about the max in both cases
0314     azOff = azOffSave;
0315     azI = azIs(:,mm);
0316     f = find(azI == max(azI));
0317     azOff = azOff(azInd);
0318     ampVal = nan;
0319     stopAz = 0;
0320     
0321     azIndex = 1;
0322     while(stopAz==0)
0323       azFit = find(azInd);
0324       if(f<=testVals(azIndex) | (f+testVals(azIndex))>=length(azI))
0325         badAzPoint = 1;
0326         azFit = [];
0327         xoff = nan;
0328         azIndex = azIndex+1;
0329         betaAz(1:3) = nan;
0330       else
0331         azFit = azOff((f-testVals(azIndex):f+testVals(azIndex)));
0332         azIFit = azI(f-testVals(azIndex):f+testVals(azIndex));
0333         indgood= zeros(size(azI));
0334         indgood(f-testVals(azIndex):f+testVals(azIndex)) = 1;
0335         
0336         
0337         display(sprintf('testing value: %d', testVals(azIndex)));
0338         
0339         OPTIONS = optimset('Display', 'off', 'MaxIter', 200000, ...
0340         'MaxFunEvals', 20000, 'TolX', 0.00005, 'TolFun', 0.00005);
0341         LB = [0 azFit(testVals(azIndex))-0.2  -10 -inf -inf];
0342         UB = [max(azIFit) azFit(testVals(azIndex))+0.2 10 inf inf];
0343         x0 = [1 azFit(testVals(azIndex)-1) 1 1 1];
0344         [betaAz resnorm residual exitflag] = lsqnonlin(@(x) mchisq2(azIFit, gaussfit(x, azFit)), x0, LB, UB, OPTIONS);
0345         
0346         amplitude = betaAz(1);
0347         xoff      = betaAz(2);
0348         sigma     = betaAz(3);
0349         xoffIndex = find(abs(azFit - xoff) == min(abs(azFit - xoff)));
0350         ampVal    = azIFit(xoffIndex);
0351         
0352         
0353         % other terms remove a baseline
0354         if (amplitude<0 | abs(xoff)>8 | resnorm>500)
0355           badAzPoint = 1;
0356         else
0357           badAzPoint = 0;
0358         end 
0359         azFit = gaussfit(betaAz, azOff);
0360         
0361         if(abs(xoff)<4)
0362           stopAz = 1;
0363         else
0364           azIndex = azIndex+1;
0365         end
0366       end
0367       
0368       if(azIndex>length(testVals))
0369         stopAz = 1;
0370       end
0371     end
0372     % save for weighting
0373     if(azIndex >= length(testVals))
0374       azIndex = azIndex-1;
0375     end
0376     azWidth(m,mm) = testVals(azIndex);
0377     azMax(m,mm)   = ampVal;
0378     if(~badAzPoint)
0379       azBase(m,mm)  = median(azI(~indgood));
0380       azRms(m,mm)   = std(azI(~indgood));
0381     else
0382       azBase(m,mm)  = nan;
0383       azRms(m,mm)   = nan;
0384     end
0385     betaAzAll(m,mm,:) = betaAz(1:3);
0386       end
0387       
0388       
0389       if(~isempty(elInd))
0390     elOff = elOffSave;
0391     elI = elIs(:,mm);
0392     % first we remove the offset from slewing the atmosphere
0393     % from the sky dip code
0394     x = 1./(sind(dc.antenna0.servo.el(elInd)));
0395     y = abs(elI);
0396     x(isnan(y)) = [];
0397     y(isnan(y)) = [];
0398     [tau Tground] = linfit(x,y);
0399     % remove the effect
0400     elI = elI - Tground(1) - tau(1)./sind(dc.antenna0.servo.el(elInd));
0401     
0402     % fit a gaussian about the max
0403     frange = elOff(elInd);
0404     indrange = abs(frange)<4;
0405     f = find(elI(indrange) == max(elI(indrange))) + find(indrange,1);
0406     elOff2 = elOff(elInd);
0407     
0408     stopEl = 0;
0409     elIndex = 1;
0410     while(stopEl==0)
0411       elFit = find(elInd);
0412       if(f<=testVals(elIndex) | (f+testVals(elIndex))>=length(elI))
0413         badElPoint = 1;
0414         elFit = [];
0415         yoff = nan;
0416         elIndex = elIndex+1;
0417         betaEl = [nan nan nan];
0418       else
0419         elFit = elOff2((f-testVals(elIndex):f+testVals(elIndex)));
0420         elIFit = elI(f-testVals(elIndex):f+testVals(elIndex));
0421         indgood= zeros(size(elI));
0422         indgood(f-testVals(elIndex):f+testVals(elIndex)) = 1;
0423         
0424         display(sprintf('el testing value: %d', testVals(elIndex)));
0425         [betaEl] = nlinfit(elFit, elIFit, @gaussfit, [1 elFit(testVals(elIndex)) 1 1 1]);
0426         amplitude = betaEl(1);
0427         yoff      = betaEl(2);
0428         sigma     = betaEl(3);
0429         yoffIndex = find(abs(elFit - yoff) == min(abs(elFit - yoff)));
0430         ampVal    = elIFit(yoffIndex);        
0431         
0432         % other terms remove a baseline
0433         if (amplitude<0 | abs(yoff)>8)
0434           badElPoint = 1;
0435         else
0436           badElPoint = 0;
0437         end 
0438         elFit = gaussfit(betaEl, elOff2);
0439         
0440         if(abs(yoff)<1)
0441           stopEl = 1;
0442         else
0443           elIndex = elIndex+1;
0444         end
0445       end
0446       if(elIndex>length(testVals))
0447         stopEl = 1;
0448       end
0449       
0450     end
0451     % save for weighting
0452     if(elIndex >= length(testVals))
0453       elIndex = elIndex - 1;
0454     end
0455     elWidth(m,mm) = testVals(elIndex);
0456     elMax(m,mm)   = ampVal;
0457     if(~badElPoint)
0458       elBase(m,mm)   = median(elI(~indgood));
0459       elRms(m,mm)    = std(elI(~indgood));
0460     else
0461       elBase(m,mm)   = nan;
0462       elRms(m,mm)    = nan;
0463     end
0464     betaElAll(m,mm,:) = betaEl(1:3);      
0465       end
0466       
0467       display('here');
0468       
0469       display(sprintf('Point %d of %d', m, length(s)));
0470       
0471       
0472       clf
0473       %if(~any(isnan(betaAzAll(m,mm,:))))
0474       
0475       subplot(2,1,1)
0476       plot((azOff)*cos(dc.antenna0.servo.el(1)*pi/180), azI');
0477       xlabel('Az Offset (degrees)');
0478       ylabel('Intensity');
0479       if(~isempty(azFit))
0480     hold on; plot(azOff*cos(dc.antenna0.servo.el(1)*pi/180), azFit, 'r'); hold off
0481     axis([min(azOff) max(azOff) min(azI) max(azI)]);
0482       end
0483       
0484       %end
0485       
0486       %if(~any(isnan(betaElAll(m,mm,:))))
0487       subplot(2,1,2)
0488       plot(elOff(elInd), elI);
0489       xlabel('El Offset (degrees)');
0490       ylabel('Intensity');
0491       if(~isempty(elFit))
0492     hold on; plot(elOff(elInd), elFit, 'r'); hold off
0493     axis([min(elOff) max(elOff) min(elI) max(elI)]);
0494       end
0495       %end
0496       eval(sprintf('gtitle(''Pointing scan on source %s'');', dc.antenna0.tracker.source{1}));
0497       
0498       
0499       
0500       display(sprintf('Source Name: %s', dc.antenna0.tracker.source{1}));
0501       r = input('Keep This point? [y/n]: ', 's');
0502       
0503       if(strcmp(r, 'y') || strcmp(r, 'Y'))
0504     badPoint(m,mm) = 0;
0505       elseif(strcmp(r, 'n') ||  strcmp(r, 'N'))
0506     badPoint(m,mm) = 1;
0507       elseif(strcmp(r, 'k'))
0508     keyboard;
0509       end
0510     end
0511     display('here 3');
0512     
0513     clear azIs;
0514     clear elIs;
0515     
0516     pause(0.1);
0517   else
0518     display('All data in this cross have been flagged');
0519     azWidth(m,1:2) = nan;
0520     azMax(m,1:2)   = nan;
0521     azBase(m,1:2)  = nan;
0522     azRms(m,1:2)   = nan;
0523     betaAzAll(m,1:2,1:3) = nan;
0524     elWidth(m,1:2) = nan;
0525     elMax(m,1:2)   = nan;
0526     elBase(m,1:2)  = nan;
0527     elRms(m,1:2)   = nan;
0528     betaElAll(m,1:2,1:3) = nan;
0529   end
0530 end
0531 
0532 badPoint = logical(badPoint);
0533 elOff = betaElAll(:,:,2);
0534 azOff = betaAzAll(:,:,2);
0535 elOff(badPoint) = nan;
0536 azOff(badPoint) = nan;
0537 
0538 % first we assign a weight with the source.
0539 azSig = (azMax - azBase)./azRms;
0540 azWidth = abs(betaAzAll(:,:,3)/2);
0541 azErr = azWidth./azSig;
0542 
0543 
0544 elSig = (elMax - elBase)./elRms;
0545 elWidth = abs(betaElAll(:,:,3)/2);
0546 elErr   = elWidth./elSig;
0547 
0548 % throw out values where the error is large or negative.
0549 elBad = elErr < 0 | elErr > 1;
0550 azBad = azErr < 0 | azErr > 1;
0551 
0552 elOff(elBad) = nan;
0553 azOff(azBad) = nan;
0554 elErr(elBad) = nan;
0555 azErr(elBad) = nan;
0556 
0557 % if hte values for both channels aren't within 10 arcmin, toss them.
0558 dEl = abs(diff(elOff, [], 2));
0559 dAz = abs(diff(azOff, [], 2).*cos(obs.el*pi/180)');
0560 
0561 f = find(dEl>10/60 | dAz>10/60);
0562 elOff(f,:) = nan;
0563 azOff(f,:) = nan;
0564 elErr(f,:) = nan;
0565 azErr(f,:) = nan;
0566 
0567 
0568 display(sprintf('Throwing out %d points', length(f)));
0569 
0570 [off.az  off.azErr] = vwstat(azOff, azErr, 2);
0571 [off.el  off.elErr] = vwstat(elOff, elErr, 2);
0572 
0573 display('Done with all crosses');
0574 
0575 % add final offests to the observed az/el
0576 obs.az = obs.az + off.az';
0577 obs.el = obs.el + off.el';
0578 
0579 % throw out the data that is bad
0580 ind = isnan(obs.az);
0581 obs.az(ind) = [];
0582 obs.el(ind) = [];
0583 ide.az(ind) = [];
0584 ide.el(ind) = [];
0585 off.az(ind) = [];
0586 off.el(ind) = [];
0587 off.azErr(ind) = [];
0588 off.elErr(ind) = [];
0589 
0590 
0591 return;
0592 
0593 
0594 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0595 function [obs ide] = calcOffRaster(d)
0596 
0597 % I. split them up into respective scans
0598 % Ia.  find the start and stop index of each scan.
0599 [s e] = findStartStop(d.array.frame.features);
0600 
0601 
0602 for m=1:length(s)
0603   % cut the data for each scan
0604   ind = zeros(size(d.array.frame.features));
0605   ind(s(m):e(m)) = 1;
0606   ind = logical(ind);
0607   dc  = framecut(d, ind);
0608 
0609   % calculate the offset from nominal for this scan
0610   [obs.az(m) obs.el(m)] = sourceScanMap(dc, 0.5, 0);
0611   
0612   % calculate the az/el for the source
0613   ide.az(m) = mean(dc.antenna0.tracker.horiz_mount(:,1));
0614   ide.el(m) = mean(dc.antenna0.tracker.horiz_mount(:,2));
0615 end
0616 
0617 
0618 return;
0619 
0620 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0621 function [s e] = findStartStop(features)
0622 
0623 % start indices are when the difference in the feature value increases, stop
0624 % is when it decreases
0625 s = find(diff(double(features))>0);
0626 
0627 % end is when it decreases
0628 e = find(diff(double(features))<0);
0629 
0630 % last one
0631 lastOneBad = 0;
0632 if(length(e)<length(s))
0633   % the dimensions don't match
0634   if(s(1)<e(1))
0635     % if start is less than end, that means we're usually missing hte last
0636     % endpoint
0637     e(length(s)) = length(features);
0638     lastOneBad = 1;
0639   else
0640     s(1) = [];
0641   end
0642 elseif(length(s)<length(e))
0643   s(1) = [];
0644 end
0645 s = s+1;
0646 
0647 goodPoint = zeros(size(s));
0648 for m=1:length(s)
0649   thisFeat = unique(features(s(m):e(m)));
0650   if(thisFeat == 129 | thisFeat ==128)
0651     goodPoint(m) = 1;
0652   end
0653 end  
0654 
0655 if(lastOneBad==1)
0656   goodPoint(length(s)) = 0;
0657 end
0658 
0659 
0660 s = s(logical(goodPoint));
0661 e = e(logical(goodPoint));
0662 
0663 return;
0664 
0665 
0666 function x = mchisq2(dataVals, modelVals)
0667 
0668 x = dataVals - modelVals;
0669 
0670 f = find(isnan(x));
0671 x(f) = [];
0672 f = find(isinf(x));
0673 x(f) = [];
0674 
0675 if(isempty(x))
0676   display('damn it all to hell')
0677   keyboard;
0678 end
0679 
0680 return;
0681 
0682         
0683 
0684 
0685 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0686 function [obs ide off] = calcOffCross(d)
0687 
0688 % load up the current model
0689 om(:,1)=d.antenna0.tracker.flexure(:,1);
0690 om(:,2)=d.antenna0.tracker.flexure(:,2);
0691 om(:,3)=d.antenna0.tracker.tilts(:,1);
0692 om(:,4)=d.antenna0.tracker.tilts(:,2);
0693 om(:,5)=d.antenna0.tracker.tilts(:,3);
0694 om(:,6)=d.antenna0.tracker.fixedCollimation(:,1);
0695 om(:,7)=d.antenna0.tracker.fixedCollimation(:,2);
0696 om(:,8)=d.antenna0.tracker.encoder_off(:,1);
0697 om(:,9)=d.antenna0.tracker.encoder_off(:,2);
0698 om = mean(om);
0699 
0700 % calculates the offsets when doing a pointing cross:  fits for a gaussian
0701 % in x and y offsets, converts to offsets in az/el, and calculates the
0702 % observed location of the source
0703 % should be very similar to the sza
0704 
0705 [s e] = findCrossStartStop(d);
0706 
0707 % now for each of these fields, we split them up into the x and y offsets,
0708 % fit a gaussian, and find the offsets.  convert offsets to az/el offsets,
0709 % and be done.
0710 for m=1:length(s)
0711   % cut the data for each scan
0712   ind = zeros(size(d.array.frame.features));
0713   ind(s(m):e(m)) = 1;
0714   ind = logical(ind);
0715   dc  = framecut(d, ind);
0716   
0717   [intStart intStop] = findCrossStepStartStop(dc);
0718   if(length(intStart)~=26)
0719     display('You may have some mixed up data.  not the right number of points');
0720     keyboard;
0721   end
0722   
0723   for n=1:length(intStart)
0724     ind = zeros(size(dc.array.frame.features));
0725     ind(intStart(n):intStop(n)) = 1;
0726     ind = logical(ind);
0727     dcc = framecut(dc, ind);
0728     
0729     offsets(n,:) = median(dcc.antenna0.tracker.sky_xy_off);
0730     meanIs(n,:) = mean(dcc.antenna0.receiver.data(:,[1 6]));
0731     stdIs(n,:)  = std(dcc.antenna0.receiver.data(:,[1 6]));
0732     
0733     ideal(n,:) = dcc.antenna0.tracker.horiz_topo(2,:);
0734     observed(n,:) = dcc.antenna0.tracker.horiz_mount(2,:) + ...
0735     dcc.antenna0.tracker.encoder_off(2,:);
0736   end
0737   fy = find(offsets(:,2)~=0);
0738   yIndex = fy(1):last(fy);
0739   fx = find(offsets(:,1)~=0);
0740   xIndex = fx(1):last(fx);
0741   
0742   % fit a gaussian to both data sets.
0743   xOffVals = offsets(xIndex,1);
0744   yOffVals = offsets(yIndex,2);
0745   xIntVals = abs(meanIs(xIndex,:));
0746   yIntVals = abs(meanIs(yIndex,:));
0747   xErrVals = stdIs(xIndex,:);
0748   yErrVals = stdIs(yIndex,:);
0749 
0750   % find the x offset
0751   x0 = [max(xIntVals(:,1)) 0 0.2];
0752   ind = abs(xOffVals)<=0.5;
0753   [px(1,:) pex(1,:)] = matmin('chisq', x0, [], 'gauss', ...
0754       xIntVals(ind,1), xErrVals(ind,1), xOffVals(ind));
0755   [px(2,:) pex(2,:)] = matmin('chisq', x0, [], 'gauss', ...
0756       xIntVals(ind,2), xErrVals(ind,2), xOffVals(ind));
0757   fullX = min(xOffVals):0.01:max(xOffVals);
0758   aa(:,1) = gauss(px(1,:), fullX);
0759   aa(:,2) = gauss(px(2,:), fullX);
0760   offset.x(m) = vwstat(px(:,2), pex(:,2).^2);
0761 
0762   elFit = xOffVals(ind);
0763   elIFit = xIntVals(ind,2);
0764 
0765   [betaEl] = nlinfit(elFit, elIFit, @gaussfit, [1 0 1 1 1]);
0766   amplitude = betaEl(1);
0767   yoff      = betaEl(2);
0768   sigma     = betaEl(3);
0769   xoffvals  = min(xOffVals):0.01:max(xOffVals);
0770   elFitFin = gaussfit(betaEl, xoffvals);
0771   
0772   y0 = [max(yIntVals(:,1)) 0 0.2];
0773   ind = abs(yOffVals)<=0.5;
0774   [py(1,:) pey(1,:)] = matmin('chisq', y0, [], 'gauss', ...
0775       yIntVals(ind,1), yErrVals(ind,1), yOffVals(ind));
0776   [py(2,:) pey(2,:)] = matmin('chisq', y0, [], 'gauss', ...
0777       yIntVals(ind,2), yErrVals(ind,2), yOffVals(ind));  
0778   fullY = min(yOffVals):0.01:max(yOffVals);
0779   bb(:,1) = gauss(py(1,:), fullY);
0780   bb(:,2) = gauss(py(2,:), fullY);  
0781   offset.y(m) = vwstat(py(:,2), pey(:,2).^2);
0782 
0783   if(any(abs(px(:,2)) > max(xOffVals)))
0784     display('Suspect cross');
0785     badXPoint = 1;
0786   else
0787     badXPoint = 0;
0788   end
0789   if(abs(diff(px(:,2)))>0.05)
0790     display('Best Fit X offset does not agree between channels');
0791     badXPoint = 1;
0792   end
0793   
0794 
0795   if(any(abs(py(:,2)) > max(xOffVals)))
0796     display('Suspect cross');
0797     badYPoint = 1;
0798   else
0799     badYPoint = 0;
0800   end
0801   if(abs(diff(py(:,2)))>0.05)
0802     display('Best Fit Y offset does not agree between channels');
0803     badYPoint = 1;
0804   end
0805   
0806   
0807   
0808   % plot the results to see the fit.
0809   meanX = mean(xIntVals);
0810   xIntVals(:,1) = xIntVals(:,1) - meanX(1);
0811   xIntVals(:,2) = xIntVals(:,2) - meanX(2);
0812   aa(:,1) = aa(:,1) - meanX(1);
0813   aa(:,2) = aa(:,2) - meanX(2);
0814   
0815   meanY = mean(yIntVals);
0816   yIntVals(:,1) = yIntVals(:,1) - meanY(1);
0817   yIntVals(:,2) = yIntVals(:,2) - meanY(2);
0818   bb(:,1) = bb(:,1) - meanY(1);
0819   bb(:,2) = bb(:,2) - meanY(2);
0820 
0821   figure(1)
0822   subplot(2,1,1)
0823   plot(xOffVals, xIntVals, '*-');
0824   hold on
0825   plot(fullX, aa, 'r')
0826   axis([min(xOffVals) max(xOffVals) min(xIntVals(:)) max(xIntVals(:))]);
0827   legend('Channel 1', 'Channel 2', 'Fit Gaussian');
0828   xlabel('X offset');
0829   ylabel('Counts');
0830   hold off;
0831   
0832   subplot(2,1,2)
0833   plot(yOffVals, yIntVals, '*-');
0834   hold on
0835   plot(fullY, bb, 'r')
0836   axis([min(yOffVals) max(yOffVals) min(yIntVals(:)) max(yIntVals(:))]);
0837   legend('Channel 1', 'Channel 2', 'Fit Gaussian');
0838   xlabel('Y offset');
0839   ylabel('Counts');
0840   hold off;
0841 
0842   display(sprintf('Source Name: %s', dc.antenna0.tracker.source{1}));
0843   eval(sprintf('display(''X offset: %3.4f, Y offset: %3.4f'');', ...
0844       offset.x(m), offset.y(m)));
0845   if(badXPoint | badYPoint)
0846     badPoint = 1;
0847     display('Cross Suspect');
0848   end
0849   
0850   
0851   % check if we want to keep this source
0852 %  r = input('Keep This point? [y/n]: ', 's');
0853 %
0854 %  if(strcmp(r, 'y') || strcmp(r, 'Y'))
0855 %    badPoint = 0;
0856 %  elseif(strcmp(r, 'n') ||  strcmp(r, 'N'))
0857 %    badPoint = 1;
0858 %  elseif(strcmp(r, 'k'))
0859 %    keyboard;
0860 %  end
0861   
0862   if(~badPoint)
0863     % now for the ideal and observed
0864     f = find(offsets(:,1)==0 & offsets(:,2)==0);
0865     
0866     ide.az(m) = ideal(f(1),1);
0867     ide.el(m) = ideal(f(1),2);
0868     obs.az(m) = observed(f(1),1);
0869     obs.el(m) = observed(f(1),2);
0870     
0871     % convert great circle offsets into az/el offsets
0872     
0873     a.az = obs.az(m);
0874     a.el = obs.el(m);
0875     b = xyoff_to_azel(a, offset.x(m), offset.y(m));
0876     off.x(m) = offset.x(m);
0877     off.y(m) = offset.y(m);
0878     obs.az(m) = b.az;
0879     obs.el(m) = b.el;
0880   
0881   else
0882     ide.az(m) = nan;
0883     ide.el(m) = nan;
0884     obs.az(m) = nan;
0885     obs.el(m) = nan;
0886     off.x(m) = offset.x(m);
0887     off.y(m) = offset.y(m);
0888   end
0889 end  
0890 
0891 
0892 return;
0893 
0894 function [s e] = findCrossStartStop(d)
0895 % check for when the radio pointing flag switches.
0896 
0897 index = d.index.radio_point_cross.slow;
0898 
0899 s = find(diff(index)>0)+1;
0900 e  = find(diff(index)<0);
0901 
0902 % that's it.
0903 if(length(e)<length(s))
0904   % the dimensions don't match
0905   if(s(1)<e(1))
0906     % if start is less than end, that means we're usually missing hte last
0907     % endpoint
0908     %    e(length(s)) = length(d.array.frame.features);
0909     % cut that last start out.
0910     s(length(s)) = [];
0911   else
0912     s(1) = [];
0913   end
0914 elseif(length(s)<length(e))
0915   s(1) = [];
0916 end
0917 
0918 return;
0919 
0920 
0921 function [s e] = findCrossStepStartStop(d)
0922 % check for when the radio pointing flag switches.
0923 
0924 index = d.index.source.slow;
0925 
0926 s = find(diff(index)>0)+1;
0927 e  = find(diff(index)<0);
0928 
0929 % that's it.
0930 if(length(e)<length(s))
0931   % the dimensions don't match
0932   if(s(1)<e(1))
0933     % if start is less than end, that means we're usually missing hte last
0934     % endpoint
0935     e(length(s)) = length(d.array.frame.features);
0936   else
0937     s(1) = [];
0938   end
0939 elseif(length(s)<length(e))
0940   s(1) = [];
0941 end
0942 
0943 return;
0944 
0945 
0946 %%%%%%%%%%%%%%%
0947 function azel=xyoff_to_azel(azel,x,y)
0948 % Apply xy sky offset to az,el coord
0949 % following Martin Shepherd in tracker.c
0950 % Convert all to radians
0951 d2r=pi/180;
0952 azel.az=azel.az*d2r; azel.el=azel.el*d2r;
0953 x=x*d2r; y=y*d2r;
0954 
0955 % Convert x,y to r,theta on surface of sphere
0956 t=atan2(sin(x),tan(y));
0957 r=acos(cos(x).*cos(y));
0958 
0959 % Apply offset
0960 e=azel.el;
0961 top=sin(t).*sin(r);
0962 bot=cos(e).*cos(r)-sin(e).*sin(r).*cos(t);
0963 azel.az=azel.az+atan2(top,bot);
0964 azel.el=asin(sin(e).*cos(r)+cos(e).*sin(r).*cos(t));
0965 
0966 % Back to deg
0967 azel.az=azel.az/d2r; azel.el=azel.el/d2r;
0968 
0969 return