calcOffScanMel

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

  function [obs ide off] = calcOffScan(d, [plotparams, field, maindir])


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

0001 function [obs ide off out] = calcOffScan(d, plotparams, field, maindir)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function [obs ide off] = calcOffScan(d, [plotparams, field, maindir])
0006 %
0007 %
0008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0009 
0010 if nargin==1
0011   plotparams.plot = 1;
0012   plotparams.interactive = 1;
0013   plotparams.save = 0;
0014   field = [];
0015   maindir = [];
0016 end
0017 
0018 setPlotDisplay(plotparams.plot);
0019 
0020 % load up the current model
0021 om(:,1)=d.antenna0.tracker.flexure(:,1);
0022 om(:,2)=d.antenna0.tracker.flexure(:,2);
0023 om(:,3)=d.antenna0.tracker.tilts(:,1);
0024 om(:,4)=d.antenna0.tracker.tilts(:,2);
0025 om(:,5)=d.antenna0.tracker.tilts(:,3);
0026 om(:,6)=d.antenna0.tracker.fixedCollimation(:,1);
0027 om(:,7)=d.antenna0.tracker.fixedCollimation(:,2);
0028 om(:,8)=d.antenna0.tracker.encoder_off(:,1);
0029 om(:,9)=d.antenna0.tracker.encoder_off(:,2);
0030 om = mean(om);
0031 
0032 % I. split them up into respective scans
0033 % Ia.  find the start and stop index of each scan.
0034 [s e] = findStartStop(d.index.radio_point_scan.slow);
0035 betaAzAll = nan(length(s), 2, 3);
0036 betaElAll = nan(length(s), 2, 3);
0037 
0038 
0039 if(isempty(e) | isempty(s))
0040   obs.az = nan;
0041   obs.el = nan;
0042   ide.az = nan;
0043   ide.el = nan;
0044   off.az = nan;
0045   off.el = nan;
0046   goodPoint = [0 0];
0047   display('No crosses');
0048   s = [];
0049   return;
0050   
0051 end
0052 
0053 if(e(1)<s(1))
0054   obs.az = nan;
0055   obs.el = nan;
0056   ide.az = nan;
0057   ide.el = nan;
0058   off.az = nan;
0059   off.el = nan;
0060   goodPoint = [0 0];
0061   display('No crosses');
0062   s = [];
0063   return;
0064 end
0065 
0066 
0067 % determine the index of the intensity channel
0068 switch (size(d.antenna0.receiver.data,2))
0069   case 6
0070     intIndex = [1 6];
0071 
0072   case 8
0073     intIndex = [1 8];
0074     
0075   case 10
0076     intIndex = [1 9];
0077 end
0078 
0079 plotIndex = 0;
0080 
0081 % let's check if all the powers are negative.  if they are, we do an
0082 % absolute value.
0083 if(all((all(d.antenna0.receiver.data(:,intIndex)<0))))
0084   display('WARNING:  Your data are negative');
0085   display('Taking the absolute value only for fitting offsets');
0086   display('This will ONLY change the sign of the data in the pointing plots');
0087   display('All other analysis will NOT be affected');
0088   d.antenna0.receiver.data = abs(d.antenna0.receiver.data);
0089 end
0090 
0091 for m=1:length(s)
0092   doFit = 1;
0093   
0094   % cut the data for each scan
0095   ind = zeros(size(d.array.frame.features));
0096   ind(s(m):e(m)) = 1;
0097   ind = logical(ind);
0098   dc  = framecut(d, ind);
0099   aa  = unique(dc.antenna0.tracker.source);
0100   sourceName{m} = aa{1};
0101   timeVal(m)    = dc.array.frame.utc(1);
0102   
0103   % values for the pointing model
0104   idevals = dc.antenna0.tracker.horiz_topo(1,:);
0105   [obsvals(1) obsvals(2)] = pointing_model(om, idevals(1), idevals(2));
0106   ide.az(m) = idevals(1);
0107   ide.el(m) = idevals(2);
0108   obs.az(m) = obsvals(1);
0109   obs.el(m) = obsvals(2);
0110   
0111   % This one fits for 2D, so we keep the az/el scans and fit together.
0112 
0113   % get the offsets from the control system.
0114   azApp = interp1(dc.antenna0.tracker.utc, ...
0115       dc.antenna0.tracker.horiz_topo(:,1), dc.antenna0.receiver.utc);
0116   azOffSave = azApp - 180/pi*unwrap(dc.antenna0.servo.apparent(:,1)*pi/180);
0117   azOffSave = -azOffSave;
0118   
0119   elApp = interp1(dc.antenna0.tracker.utc, ...
0120       dc.antenna0.tracker.horiz_topo(:,2), dc.antenna0.receiver.utc);
0121   elOffSave = elApp - dc.antenna0.servo.apparent(:,2);
0122   elOffSave = -elOffSave;
0123   
0124   % put it in sky angles.
0125   azOffSave = azOffSave.*cosd(obs.el(m));
0126   
0127   % assuming the slowest we would ever scan is at 0.2 degrees per second, a
0128   % az scan happens when the az speed is mostly constant, increasing, and
0129   % greater than 0.2 degrees/s, which is
0130   azInd = (deriv(dc.antenna0.servo.az))>0.2/100;
0131   [startaz endaz] = findStartStop(azInd);
0132   azInd = zeros(size(azInd));
0133   if(~isempty(startaz))
0134     f = find( (endaz - startaz) == max(endaz - startaz));
0135     azInd(startaz(f):endaz(f)) = 1;
0136   else
0137     doFit = 0;
0138   end
0139   azInd = logical(azInd);
0140 
0141   % same goes for elevation
0142   elInd = deriv(dc.antenna0.servo.el)>0.2/100;
0143   [startel endel] = findStartStop(elInd);
0144   elInd = zeros(size(elInd));
0145   if(~isempty(startel))
0146     % no good el scan
0147     f = find( (endel-startel) == max(endel - startel));
0148     elInd(startel(f):endel(f)) = 1;
0149   else
0150     doFit = 0;
0151   end
0152   elInd = logical(elInd);
0153 
0154   % combine the two that we're interested in
0155   fitInd   = azInd | elInd;
0156   fitFlags = dc.flags.fast(fitInd, [1 3]);
0157   azFlags  = dc.flags.fast(azInd, [1 3]);
0158   elFlags  = dc.flags.fast(elInd, [1 3]);
0159   
0160   % setup mel's variables
0161   timeAz = dc.antenna0.receiver.utc(azInd);
0162   if(~isempty(timeAz))
0163     timeAz = (timeAz - timeAz(1)).*24*60;
0164   end
0165   skyAz  = azOffSave(azInd);
0166   newDataAz = dc.antenna0.receiver.data(azInd, intIndex);
0167   
0168   timeEl = dc.antenna0.receiver.utc(elInd);
0169   if(~isempty(timeEl))
0170     timeEl = (timeEl - timeEl(1)).*24*60;
0171   end
0172   skyEl  = elOffSave(elInd);
0173   newDataEl = dc.antenna0.receiver.data(elInd, intIndex);
0174   
0175   sigmaFit = 0.73/(2*sqrt(2*log(2)));
0176 
0177   % we should check we actually have a scan.
0178   if(length(find(elInd)) < 1500  | length(find(azInd)) < 1500)
0179     doFit = 0;
0180   end
0181 
0182   % check how much of the data are flagged
0183   elflag1  = length(find(elFlags(:,1)))./length(find(elInd));
0184   elflag2  = length(find(elFlags(:,2)))./length(find(elInd));
0185   azflag1  = length(find(azFlags(:,1)))./length(find(azInd));
0186   azflag2  = length(find(azFlags(:,2)))./length(find(azInd));
0187   
0188   if(elflag1 > 0.7 & elflag2 > 0.7)
0189     doFit = 0;
0190   end
0191   if(azflag1 > 0.7 & azflag2 > 0.7)
0192     doFit = 0;
0193   end
0194   
0195   
0196   if(doFit)
0197     % Calculate the baseline
0198     offEl = remove_background(timeEl, newDataEl, 'StepSize', 0.1, 'WindowSize', 1);
0199     f1 = find(offEl(:,1) == max(offEl(:,1)));
0200     f2 = find(offEl(:,2) == max(offEl(:,2)));
0201     elOff = mean(skyEl([f1 f2]));
0202     straightEl = newDataEl - offEl;
0203     % get the slope of the background signal and other parameters
0204     Pel(1,:) = polyfit([skyEl(1) last(skyEl)], [straightEl(1,1) ...
0205       straightEl(length(timeEl),1)], 1);
0206     Pel(2,:) = polyfit([skyEl(1) last(skyEl)], [straightEl(1,2) ...
0207       straightEl(length(timeEl),2)], 1);
0208     gradEl = Pel(:,1)';
0209     offsetEl = min(newDataEl);
0210     f = find(abs(skyEl)==min(abs(skyEl)));
0211     peakEl = newDataEl(f,:);
0212   
0213     
0214     
0215     offAz = remove_background(timeAz, newDataAz, 'StepSize', 0.1, 'WindowSize', 1);
0216     f1 = find(offAz(:,1) == max(offAz(:,1)));
0217     f2 = find(offAz(:,2) == max(offAz(:,2)));
0218     az_off = mean(skyAz([f1 f2]));
0219     straightAz = newDataAz - offAz;
0220     % get the slope of the background signal
0221     Paz(1,:) = polyfit([skyAz(1) last(skyAz)], [straightAz(1,1) ...
0222       straightAz(length(timeAz),1)], 1);
0223     Paz(2,:) = polyfit([skyAz(1) last(skyAz)], [straightAz(1,2) ...
0224       straightAz(length(timeAz),2)], 1);
0225     gradAz = Paz(:,1)';
0226     offsetAz = min(newDataAz);
0227     f = find(abs(skyAz)==min(abs(skyAz)));
0228     peakAz = newDataAz(f,:);
0229     azOff  = skyAz(f);
0230     
0231     % set up the parameters for the fit
0232     peakVal = mean([peakAz; peakEl]);
0233     
0234     % beta has two rows, one for each channel
0235     Beta = [peakVal', [azOff; azOff], [sigmaFit; sigmaFit], [elOff;elOff], [sigmaFit; sigmaFit], gradAz', gradEl', ...
0236       offsetAz', offsetEl'];
0237     Beta(isnan(Beta)) = 0;
0238     
0239     % first let's get our data ready
0240     skylocs = [skyAz; skyEl];
0241     lenAz   = length(skyAz);
0242     dataFit = [newDataAz; newDataEl];
0243     
0244     % we don't want sidelobes in the fit
0245     indLobes = abs(skylocs)>0.5 & abs(skylocs)<1.7;
0246     % if it's CygA, we need to check for CygX
0247     indCygx = zeros(size(indLobes));
0248     if(sourceCorrespondance(unique(dc.antenna0.tracker.source))==1)
0249       % check the mean on either side of the peak have similar values
0250       azDat = newDataAz;
0251       azDat(azFlags) = nan;
0252       f1       = find(azDat(:,1) == max(azDat(:,1)));
0253       if(min(skyAz)< (skyAz(f1)-1.3))
0254     ind_no_lobe = skyAz < (skyAz(f1) - 1.1);
0255       else
0256     ind_no_lobe = skyAz < skyAz(f1);
0257       end
0258       m1_first = sqrt(nanvar(azDat(ind_no_lobe,:)));
0259       if(max(skyAz)> (skyAz(f1)+ 1.3))
0260     ind_no_lobe = skyAz > (skyAz(f1) + 1.1);
0261       else
0262     ind_no_lobe = skyAz > skyAz(f1);
0263       end
0264       m1_sec   = sqrt(nanvar(azDat(ind_no_lobe,:)));
0265       m1_diff  = ((m1_first-m1_sec)./m1_first)*100;
0266       
0267       elDat = newDataEl;
0268       elDat(elFlags) = nan;
0269       f2       = find(elDat(:,1) == max(elDat(:,1)));
0270       if(min(skyEl)< (skyEl(f2)-1.3))
0271     ind_no_lobe = skyEl < (skyEl(f2) - 1.1);
0272       else
0273     ind_no_lobe = skyEl < skyEl(f2);
0274       end
0275       m2_first = sqrt(nanvar(elDat(ind_no_lobe,:)));
0276       if(max(skyEl) > (skyEl(f2) + 1.3))
0277     ind_no_lobe = skyEl > (skyEl(f2) + 1.1);
0278       else
0279     ind_no_lobe = skyEl > skyEl(f2);
0280       end
0281       m2_sec   = sqrt(nanvar(elDat(ind_no_lobe,:)));
0282       m2_diff  = ((m2_first-m2_sec)./m2_first)*100;      
0283       
0284       indCygaz = zeros(size(skyAz));
0285       indCygel = zeros(size(skyEl));
0286       if(any(abs(m1_diff) > 35))
0287     % cygnux X is present in your az scan
0288     if(mean(m1_diff) > 0)
0289       % it's present in first half
0290       % don't get rid of the lobe itself.
0291       ff = find(skyAz > (skyAz(f1) - 0.5), 1); 
0292       indCygaz(1:ff) = 1;
0293     else
0294       ff = find(skyAz > (skyAz(f1) + 0.5),1);
0295       indCygaz(ff:lenAz) = 1;
0296     end
0297       end
0298       if(any(abs(m2_diff) > 35))
0299     % Cygnus X is present in your el scan
0300     if(mean(m2_diff) > 0)
0301       % present in first half
0302       ff = find(skyEl > (skyEl(f2)-0.5),1);
0303       indCygel(1:ff) = 1;
0304     else
0305       ff = find(skyEl > (skyEl(f2)+0.5),1);
0306       indCygel(ff:length(skyEl)) = 1;
0307     end
0308       end
0309       indCygx = [indCygaz; indCygel];
0310     else
0311       indCygx = zeros(size(indLobes));
0312     end
0313   end
0314 
0315   % Mel's Method
0316   for mm=1:2
0317     if(doFit)
0318       % do not select flagged data
0319       indNoFit = indCygx | indLobes | fitFlags(:,mm);
0320       indFit   = abs(skylocs)<5 & ~indNoFit;
0321     else
0322       indFit = 0;
0323     end
0324     if(length(find(indFit)) > 1000)
0325       % we fit if less than half the data are flagged
0326       thisData = dataFit(:,mm);
0327 %            [betanew(mm,:), resid1, J1, cov1, mse] = nlinfit([lenAz; skylocs; indFit], thisData(indFit), ...
0328 %            @gauss2D_fit2, Beta(mm,:));
0329 %            errbeta(mm,:) = sqrt(abs(diag(cov1)));
0330 %      LB = [0 -3 0.1 -3 0.1 -inf -inf -inf -inf];
0331 %      UB = [inf 3 2 3 2 inf inf inf inf];
0332 %
0333 %      [betanew(mm,:), resnorm, resid1, exitflag, output, lambda, J1] = ...
0334 %      lsqcurvefit('gauss2D_fit2', Beta(mm,:), [lenAz; skylocs; indFit], ...
0335 %      thisData(indFit), LB, UB);
0336 %      errbeta = ones(size(betanew));
0337       display('calculating scan off');
0338       % let's try it with matmin because it gives us errors on the fits!!!!
0339       [betanew(mm,:), errbeta(mm,:), gof, stat] = matmin('chisq', Beta(mm,:), [], ...
0340       'gauss2D_fit2', thisData(indFit), ones(size(thisData(indFit))), ...
0341       [lenAz; skylocs; indFit]);
0342 
0343       fitModel = gauss2D_fit2(betanew(mm,:), [lenAz; skylocs; indFit]);
0344       resid1   = thisData(indFit) - fitModel;
0345       % calculate our models
0346       azScanModel = gauss2D_draw(betanew(mm,:), skyAz, zeros(size(skyAz)));
0347       elScanModel = gauss2D_draw(betanew(mm,:), zeros(size(skyEl)), skyEl);
0348       % flag if it's bad
0349       % we'll flag on three conditions:
0350       % 1.  the x/y offset is more than 1 degree
0351       % 2.  the residuals on the fit data are high
0352       % 3.  the width of the gaussian doesn't match our beam to within 20%
0353       thisPointBad = 0;
0354       badString = 'GOOD';
0355       if(any(abs(betanew(mm,[2 4])) >1))
0356     display('Flagging point');
0357     display('Offset too great');
0358     thisPointBad = 1;
0359       end
0360       residPercent = mean(abs(resid1))./min(thisData)*100;
0361       if(mean(abs(resid1)) > 0.3)
0362     display('Flagging point');
0363     display('Residuals too high');
0364     mean(abs(resid1))
0365     thisPointBad = 1;
0366       end
0367       beamRatio = mean(betanew(mm, [3 5]))/sigmaFit;
0368       if( beamRatio < 0.3 | beamRatio > 2)
0369     display('Flagging point');
0370     display('Width does not match expected beam');
0371     display(sprintf('beam ratio is: %3.2f', beamRatio));
0372 %    keyboard;
0373     thisPointBad = 1;
0374       end
0375       
0376       % select the information we want to keep
0377       indAzFit     = indFit(1:lenAz);
0378       indElFit     = indFit(lenAz+1:length(indFit));
0379       thisDataAz   = newDataAz(:,mm);
0380       thisDataEl   = newDataEl(:,mm);
0381       
0382       %  we need the az and el offset (az as in real azimuth, not sky angle)
0383       %  need the x and y widths
0384       %  need the peak value, baseline (min) value and rms of residuals
0385       xwidth(m,mm) = betanew(mm,3);
0386       ywidth(m,mm) = betanew(mm,5);
0387       xoff(m,mm)   = betanew(mm,2);
0388       yoff(m,mm)   = betanew(mm,4);
0389       peakheight(m,mm) = betanew(mm,1);
0390       errpeak(m,mm)= errbeta(mm,1);
0391       errxoff(m,mm)= errbeta(mm,2);
0392       erryoff(m,mm)= errbeta(mm,4);
0393       eloff(m,mm)  = betanew(mm,4);
0394       azoff(m,mm)  = betanew(mm,2)./cosd(obs.el(m));
0395       faz          = find( abs(skyAz - xoff(m,mm)) == min(abs(skyAz - xoff(m,mm))));
0396       azpeak(m,mm) = thisDataAz(faz);
0397       thisDataAz(~indAzFit) = nan;
0398       azmin(m,mm)  = min(thisDataAz);
0399       azrms(m,mm)  = rms(thisDataAz - azScanModel);
0400       fel          = find( abs(skyEl - yoff(m,mm)) == min(abs(skyEl - yoff(m,mm))));
0401       elpeak(m,mm) = thisDataEl(fel);
0402       thisDataEl(~indElFit) = nan;
0403       elmin(m,mm)  = min(thisDataEl);
0404       elrms(m,mm)  = rms(thisDataEl - elScanModel);
0405       
0406       % next we plot the point
0407       display(sprintf('Plotting Point %d of %d', m, length(s)));
0408       clf
0409       skywrap = unwrap(skyAz);
0410       subplot(2,1,1)
0411       plot(skywrap, newDataAz(:,mm));
0412       xlabel('Sky Az Offset (degrees)');
0413       ylabel('Intensity');
0414       if(1)
0415     hold on; plot(skywrap, azScanModel, 'r'); hold off
0416     axis([min([skywrap; skyEl]) max([skyAz; skyEl]) min(newDataAz(:,mm)) max(newDataAz(:,mm))]);
0417       end
0418       subplot(2,1,2)
0419       plot(skyEl, newDataEl(:,mm));
0420       xlabel('El Offset (degrees)');
0421       ylabel('Intensity');
0422       if(1)
0423     hold on; plot(skyEl, elScanModel, 'r'); hold off;
0424     axis([min([skyAz; skyEl]) max([skyAz; skyEl]) min(newDataEl(:,mm)) max(newDataEl(:,mm))]);
0425       end
0426       eval(sprintf('title(''Intensity Channel # %d'');', mm));
0427       eval(sprintf('gtitle(''scan on source %s:'', 0.96, 2);', ...
0428       dc.antenna0.tracker.source{1}));
0429       
0430       if(plotparams.interactive)
0431     display(sprintf('Source Name: %s', dc.antenna0.tracker.source{1}));
0432     r = input('Keep This point? [y/n]: ', 's');
0433     if(strcmp(r, 'y') || strcmp(r, 'Y'))
0434       badPoint(m,mm) = 0;
0435       badString = 'GOOD';
0436     elseif(strcmp(r, 'n') ||  strcmp(r, 'N'))
0437       badPoint(m,mm) = 1;
0438       badString = 'BAD';
0439     elseif(strcmp(r, 'k'))
0440       keyboard;
0441     end
0442       else
0443     badPoint(m,mm) = thisPointBad;
0444     if(badPoint(m,mm))
0445       badString = 'BAD ';
0446     else
0447       badString = 'GOOD';
0448     end
0449       end
0450       hold off;
0451       eval(sprintf('gtitle(''%s'', 0.96, 3);',...
0452       badString));
0453       
0454       if(plotparams.save)
0455     plotIndex = plotIndex+1;
0456     if(isempty(maindir))
0457       maindir = getMainDir(d, field);
0458     end
0459     dbclear if error
0460     set(gcf,'paperposition',[0 0 6.0 9.0])
0461     filename = sprintf('%s/pointing/fig%d.png', maindir, plotIndex);
0462     eval(sprintf('print -dpng -r70 %s', filename));
0463     dbstop if error
0464       end
0465       
0466     else  % if we're not fitting because we don't have enough data.
0467       display('The majority of the data in this cross was flagged');
0468       xwidth(m,mm) = nan;
0469       ywidth(m,mm) = nan;
0470       xoff(m,mm)   = nan;
0471       yoff(m,mm)   = nan;
0472       peakheight(m,mm) = nan;
0473       errpeak(m,mm)= nan;
0474       errxoff(m,mm)= nan;
0475       erryoff(m,mm)= nan;
0476       eloff(m,mm)  = nan;
0477       azoff(m,mm)  = nan;
0478       azmin(m,mm)  = nan;
0479       azpeak(m,mm) = nan;
0480       azrms(m,mm)  = nan;
0481       elmin(m,mm)  = nan;
0482       elpeak(m,mm) = nan;
0483       elrms(m,mm)  = nan;
0484       badPoint(m,mm) = 1;
0485     end  % with check of whether we should fit
0486   
0487   end   % of loop over intensity channels
0488   pause(0.1);
0489   
0490 end  % of loop over source scans
0491 
0492 % first we assign a weight with the source.
0493 azSig = (azpeak - azmin)./azrms;
0494 azErr = 1./azSig;
0495 
0496 elSig = (elpeak - elmin)./elrms;
0497 elErr = 1./elSig;
0498 
0499 % throw out values where the error is large or negative.
0500 elBad = elErr < 0 | elErr > 1;
0501 azBad = azErr < 0 | azErr > 1;
0502 
0503 indBad = logical(badPoint);
0504 indBad(elBad) = 1;
0505 indBad(azBad) = 1;
0506 
0507 eloff(indBad) = nan;
0508 azoff(indBad) = nan;
0509 
0510 % if hte values for both channels aren't within 10 arcmin, toss them.
0511 dEl = abs(diff(eloff, [], 2));
0512 dAz = abs(diff(azoff, [], 2).*cosd(obs.el)');
0513 f = find(dEl>15/60 | dAz>15/60);
0514 indBad(f,:) = 1;
0515 
0516 % save all the data
0517 out.time = timeVal;
0518 out.name = sourceName;
0519 out.az   = ide.az;
0520 out.el   = ide.el;
0521 
0522 badPoints = indBad;
0523 
0524 indBad = sum(indBad,2)==2;
0525 
0526 display(sprintf('Throwing out %d points', length(find(indBad))));
0527 azoff(badPoints) = nan;
0528 eloff(badPoints) = nan;
0529 
0530 [off.az  off.azErr] = vwstat(azoff, azErr, 2);
0531 [off.el  off.elErr] = vwstat(eloff, elErr, 2);
0532 
0533 display('Done with all crosses');
0534 
0535 % add final offests to the observed az/el
0536 obs.az = obs.az + off.az';
0537 obs.el = obs.el + off.el';
0538 
0539 % throw out the data that is bad
0540 ind = isnan(obs.az);
0541 obs.az(ind) = [];
0542 obs.el(ind) = [];
0543 obs.name    = sourceName(~ind);
0544 obs.timeVal = timeVal(~ind);
0545 obs.index   = find(~ind);
0546 ide.az(ind) = [];
0547 ide.el(ind) = [];
0548 off.az(ind) = [];
0549 off.el(ind) = [];
0550 off.azErr(ind) = [];
0551 off.elErr(ind) = [];
0552 
0553 
0554 % all the outputs we might want.
0555 out.peakheight = peakheight;
0556 out.xwidth = xwidth;  
0557 out.ywidth = ywidth; 
0558 out.xoff   = xoff;   
0559 out.yoff   = yoff;
0560 out.errpeak= errpeak;
0561 out.errxoff= errxoff;
0562 out.erryoff= erryoff;
0563 out.eloff  = eloff;  
0564 out.azoff  = azoff;  
0565 out.azmin  = azmin;  
0566 out.azpeak = azpeak; 
0567 out.azrms  = azrms;  
0568 out.elmin  = elmin;  
0569 out.elpeak = elpeak; 
0570 out.elrms  = elrms;  
0571 out.azsig  = azSig;
0572 out.elsig  = elSig;
0573 out.badpt  = badPoints;
0574 
0575 
0576 return;
0577