calcOffScan

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

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


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

0001 function [obs ide off out] = calcOffScan(d, plotparams, field, maindir)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function [obs ide off out] = 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   % cut the data for each scan
0093   ind = zeros(size(d.array.frame.features));
0094   ind(s(m):e(m)) = 1;
0095   ind = logical(ind);
0096   dc  = framecut(d, ind);
0097 
0098   
0099   
0100   % find the observation az/el
0101 %  indOnSource = dc.antenna0.tracker.offSource==0;
0102 %  indAzOff    = dc.antenna0.tracker.horiz_off(:,1) == 0 & ...
0103 %      dc.antenna0.tracker.scan_off(:,1) ==0;
0104 %  indElOff    = dc.antenna0.tracker.horiz_off(:,2) == 0 & ...
0105 %      dc.antenna0.tracker.scan_off(:,2) ==0;
0106 %  f = find(indOnSource & indAzOff & indElOff);
0107 %  if(length(f)==0)
0108 %    f = find(indOnSource);
0109 %  end
0110 %  if(length(f)>1)
0111 %    obsvals = (dc.antenna0.tracker.horiz_mount(f(1),:)) + dc.antenna0.tracker.encoder_off(f(1),:);
0112 %    idevals = (dc.antenna0.tracker.horiz_topo(f(1),:));
0113 %
0114 %  elseif(length(f) == 1)
0115 %    obsvals = dc.antenna0.tracker.horiz_mount(f,:) + dc.antenna0.tracker.encoder_off(f,:);
0116 %    idevals = (dc.antenna0.tracker.horiz_topo(f,:));
0117 %  end
0118 
0119 
0120   idevals = dc.antenna0.tracker.horiz_topo(1,:);
0121   [obsvals(1) obsvals(2)] = pointing_model(om, idevals(1), idevals(2));
0122   ide.az(m) = idevals(1);
0123   ide.el(m) = idevals(2);
0124   obs.az(m) = obsvals(1);
0125   obs.el(m) = obsvals(2);
0126   
0127   
0128   
0129   
0130   
0131   % next we split it into the azimuth/elevation scan
0132   % for the elevation scan, we fit a line and remove the slope (effectively
0133   % a sky dip), and fit a gaussian to the remainder
0134   % for the azimuth scan, we just fit a gaussian to it.
0135 
0136   % get the offsets from the control system.
0137   azApp = interp1(dc.antenna0.tracker.utc, ...
0138       dc.antenna0.tracker.horiz_topo(:,1), dc.antenna0.receiver.utc);
0139   azOffSave = azApp - (180/pi)*(unwrap(pi/180*dc.antenna0.servo.apparent(:,1)));
0140   
0141   %  [aa, ee ] =calcAzElCs(dc, 1);
0142 %  ao = aa - dc.antenna0.servo.az;
0143 %  azOffSave = ao;
0144 %  azOffSave(azOffSave<-180) = azOffSave(azOffSave<-180) + 360;
0145 %  azOffSave(azOffSave>180) = azOffSave(azOffSave>180) - 360;
0146    azOffSave = -azOffSave;
0147 
0148   elApp = interp1(dc.antenna0.tracker.utc, ...
0149       dc.antenna0.tracker.horiz_topo(:,2), dc.antenna0.receiver.utc);
0150   elOffSave = elApp - dc.antenna0.servo.apparent(:,2);
0151   %elOffSave = ee - dc.antenna0.servo.el;
0152   elOffSave = -elOffSave;
0153   
0154   % assuming the slowest we would ever scan is at 0.2 degrees per second, a
0155   % az scan happens when the az speed is mostly constant, increasing, and
0156   % greater than 0.2 degrees/s, which is
0157   azInd = (deriv(dc.antenna0.servo.az))>0.07/100;
0158   [startaz endaz] = findStartStop(azInd);
0159   f = find( (endaz - startaz) == max(endaz - startaz));
0160   azInd = zeros(size(azInd));
0161   azInd(startaz(f):endaz(f)) = 1;
0162   azInd = logical(azInd);
0163   azIs  = dc.antenna0.receiver.data(azInd,intIndex);
0164   azIs(:,1) = smooth(azIs(:,1));
0165   azIs(:,2) = smooth(azIs(:,2));
0166 
0167   azFlags       = dc.flags.fast(azInd,[1 3]);
0168   azIs(azFlags) = nan;
0169   
0170   
0171   % if amplitude is negative, no peak in data
0172   
0173   % same goes for elevation
0174   elInd = deriv(dc.antenna0.servo.el)>0.07/100;
0175   [startel endel] = findStartStop(elInd);
0176   f = find( (endel-startel) == max(endel - startel));
0177   elInd = zeros(size(elInd));
0178   elInd(startel(f):endel(f)) = 1;
0179   elInd = logical(elInd);
0180   elIs  = dc.antenna0.receiver.data(elInd,intIndex);
0181   elIs(:,1) = smooth(elIs(:,1));
0182   elIs(:,2) = smooth(elIs(:,2));
0183   
0184   elFlags       = dc.flags.fast(elInd, [1 3]);
0185   elIs(elFlags) = nan;
0186   
0187   % if most of the scan is flagged, don't bother fitting.
0188   azFlagPer = length(find(azFlags))./length(azIs(:));
0189   elFlagPer = length(find(elFlags))./length(elIs(:));
0190   if(azFlagPer < 0.15 & elFlagPer < 0.15)
0191     doFit = 1;
0192   else
0193     doFit = 0;
0194   end
0195 
0196   aa = unique(dc.antenna0.tracker.source);
0197   sourceName{m} = aa{1};
0198   timeVal(m)    = dc.array.frame.utc(1);
0199   
0200   if(doFit)
0201       % Set the minimum, step and maximum number of data points to fit
0202     testVals = [560:30:1000];
0203     for mm=1:2
0204       
0205       if(~isempty(find(azInd)))
0206     % fit a gaussian about the max in both cases
0207     azOff = azOffSave;
0208     azI = azIs(:,mm);
0209     f = find(azI == max(azI));
0210     azOff = azOff(azInd);
0211     
0212     stopAz = 0;
0213 
0214     azIndex = 1;
0215     while(stopAz==0)
0216       azFit = find(azInd);
0217       if(f<=testVals(azIndex) | (f+testVals(azIndex))>=length(azI))
0218         badAzPoint = 1;
0219         azFit = [];
0220         xoff = nan;
0221         azIndex = azIndex+1;
0222         betaAz(1:3) = nan;
0223         ampVal = nan;
0224       else
0225         azFit = azOff((f-testVals(azIndex):f+testVals(azIndex)));
0226         azIFit = azI(f-testVals(azIndex):f+testVals(azIndex));
0227         indgood= zeros(size(azI));
0228         indgood(f-testVals(azIndex):f+testVals(azIndex)) = 1;
0229         
0230         
0231 %        display(sprintf('testing value: %d', testVals(azIndex)));
0232         
0233         OPTIONS = optimset('Display', 'off', 'MaxIter', 200000, ...
0234         'MaxFunEvals', 20000, 'TolX', 0.00005, 'TolFun', 0.00005);
0235     
0236         % MP 9-Jan-2015 Changed to use nlinfit and gaussfit_sig, plus
0237         % updated boundaries
0238         %LB = [0 azFit(testVals(azIndex))-0.2  -10 -inf -inf];
0239         %UB = [max(azIFit) azFit(testVals(azIndex))+0.2 10 inf inf];
0240         %LB = [-inf -inf  -inf -inf -inf];
0241         %UB = [inf inf inf inf inf];
0242         %x0 = [1 azFit(testVals(azIndex)-1) 1 1 1];
0243         ampguess = max(azIFit)-min(azIFit);
0244         %[betaAz resnorm residual exitflag] = lsqnonlin(@(x) mchisq2(azIFit, gaussfit(x, azFit)), x0, LB, UB, OPTIONS);
0245         [betaAz resnorm residual exitflag] = nlinfit(azFit, azIFit, @gaussfit_sig, [ampguess azFit(testVals(azIndex)) 0.3 min(azIFit) 0.0 ]);
0246         
0247         amplitude = betaAz(1);
0248         xoff      = betaAz(2);
0249         sigma     = betaAz(3); 
0250         xoffIndex = find(abs(azFit - xoff) == min(abs(azFit - xoff)));
0251         ampVal    = azIFit(xoffIndex);
0252         
0253         % other terms remove a baseline
0254         if (amplitude<0 | abs(xoff)>8 | resnorm>500)
0255           badAzPoint = 1;
0256         else
0257           badAzPoint = 0;
0258         end 
0259         azFit = gaussfit_sig(betaAz, azOff);
0260 
0261         if(abs(xoff)<4)
0262           stopAz = 1;
0263         else
0264           azIndex = azIndex+1;
0265         end
0266       end
0267       
0268       if(azIndex>length(testVals))
0269         stopAz = 1;
0270       end
0271     end
0272     % save for weighting
0273     if(azIndex >= length(testVals))
0274       azIndex = azIndex-1;
0275     end
0276     azWidth(m,mm) = testVals(azIndex);
0277     azMax(m,mm)   = ampVal;
0278     if(~badAzPoint)
0279       azBase(m,mm)  = nanmedian(azI(~indgood));
0280       azRms(m,mm)   = nanstd(azI(~indgood));
0281     else
0282       azBase(m,mm)  = nan;
0283       azRms(m,mm)   = nan;
0284     end
0285     betaAzAll(m,mm,:) = betaAz(1:3);
0286       
0287     % check one last time if it's a good point.
0288     sigDet  = (azMax(m,mm) - azBase(m,mm) )./azRms(m,mm);
0289     if(sigDet < 3 | isnan(sigDet))
0290       badAzPoint = 1;
0291       azBase(m,mm) = nan;
0292       azRms(m,mm) = nan;
0293     end
0294       end
0295       
0296       
0297       if(~isempty(elInd))
0298     elOff = elOffSave;
0299     elI = elIs(:,mm);
0300     % first we remove the offset from slewing the atmosphere
0301     % from the sky dip code
0302     x = 1./(sind(dc.antenna0.servo.el(elInd)));
0303     y = abs(elI);
0304     x(isnan(y)) = [];
0305     y(isnan(y)) = [];
0306     [tau Tground] = linfit(x,y);
0307     % remove the effect
0308     elI = elI - Tground(1) - tau(1)./sind(dc.antenna0.servo.el(elInd));
0309     
0310     % fit a gaussian about the max
0311     frange = elOff(elInd);
0312     indrange = abs(frange)<4;
0313     f = find(elI(indrange) == max(elI(indrange))) + find(indrange,1);
0314     elOff2 = elOff(elInd);
0315     
0316     stopEl = 0;
0317     elIndex = 1;
0318     while(stopEl==0)
0319       elFit = find(elInd);
0320       if(f<=testVals(elIndex) | (f+testVals(elIndex))>=length(elI))
0321         badElPoint = 1;
0322         elFit = [];
0323         yoff = nan;
0324         elIndex = elIndex+1;
0325         betaEl = [nan nan nan];
0326         ampVal = nan;
0327       else
0328         elFit = elOff2((f-testVals(elIndex):f+testVals(elIndex)));
0329         elIFit = elI(f-testVals(elIndex):f+testVals(elIndex));
0330         indgood= zeros(size(elI));
0331         indgood(f-testVals(elIndex):f+testVals(elIndex)) = 1;
0332         elguess = max(elIFit)-min(elIFit);
0333         
0334 %        display(sprintf('el testing value: %d', testVals(elIndex)));
0335         [betaEl] = nlinfit(elFit, elIFit, @gaussfit_sig,  [elguess elFit(testVals(elIndex)) 0.3 min(elIFit) 0.0 ]);
0336         amplitude = betaEl(1);
0337         yoff      = betaEl(2);
0338         sigma     = betaEl(3);
0339         yoffIndex = find(abs(elFit - yoff) == min(abs(elFit - yoff)));
0340         ampVal    = elIFit(yoffIndex);        
0341         
0342         % other terms remove a baseline
0343         if (amplitude<0 | abs(yoff)>8)
0344           badElPoint = 1;
0345         else
0346           badElPoint = 0;
0347         end 
0348         elFit = gaussfit_sig(betaEl, elOff2);
0349         
0350         if(abs(yoff)<1)
0351           stopEl = 1;
0352         else
0353           elIndex = elIndex+1;
0354         end
0355       end
0356       if(elIndex>length(testVals))
0357         stopEl = 1;
0358       end
0359       
0360     end
0361     % save for weighting
0362     if(elIndex >= length(testVals))
0363       elIndex = elIndex - 1;
0364     end
0365     elWidth(m,mm) = testVals(elIndex);
0366     elMax(m,mm)   = ampVal;
0367     if(~badElPoint)
0368       elBase(m,mm)   = nanmedian(elI(~indgood));
0369       elRms(m,mm)    = nanstd(elI(~indgood));
0370     else
0371       elBase(m,mm)   = nan;
0372       elRms(m,mm)    = nan;
0373     end
0374     betaElAll(m,mm,:) = betaEl(1:3);      
0375       end
0376       
0377 %      display('here');
0378       
0379       display(sprintf('Point %d of %d', m, length(s)));
0380       
0381       
0382       clf
0383       %if(~any(isnan(betaAzAll(m,mm,:))))
0384       subplot(2,1,1)
0385       plot( (azOff)*cos(dc.antenna0.servo.el(1)*pi/180), azI');
0386       xlabel('Az Offset (degrees)');
0387       ylabel('Intensity');
0388       if(~isempty(azFit))
0389     hold on; plot(azOff*cos(dc.antenna0.servo.el(1)*pi/180), azFit, 'r'); hold off
0390     axis([min(azOff) max(azOff) min(azI) max(azI)]);
0391       end
0392       
0393       %end
0394       
0395       %if(~any(isnan(betaElAll(m,mm,:))))
0396       subplot(2,1,2)
0397       plot(elOff(elInd), elI);
0398       xlabel('El Offset (degrees)');
0399       ylabel('Intensity');
0400       if(~isempty(elFit))
0401     hold on; plot(elOff(elInd), elFit, 'r'); hold off
0402     axis([min(elOff) max(elOff) min(elI) max(elI)]);
0403       end
0404       %end
0405 
0406       eval(sprintf('gtitle(''scan on source %s:'', 0.96, 2);', ...
0407       dc.antenna0.tracker.source{1}));
0408       
0409       if(plotparams.interactive)
0410     display(sprintf('Source Name: %s', dc.antenna0.tracker.source{1}));
0411     r = input('Keep This point? [y/n]: ', 's');
0412     
0413     if(strcmp(r, 'y') || strcmp(r, 'Y'))
0414       badPoint(m,mm) = 0;
0415       badString = 'GOOD';
0416     elseif(strcmp(r, 'n') ||  strcmp(r, 'N'))
0417       badPoint(m,mm) = 1;
0418       badString = 'BAD';
0419     elseif(strcmp(r, 'k'))
0420       keyboard;
0421     end
0422       else
0423     badPoint(m,mm) = badAzPoint | badElPoint;
0424     if(badPoint(m,mm))
0425       badString = 'BAD ';
0426     else
0427       badString = 'GOOD';
0428     end
0429       end
0430       hold off;
0431       eval(sprintf('gtitle(''%s'', 0.96, 3);',...
0432       badString));
0433       
0434       if(plotparams.save)
0435     plotIndex = plotIndex+1;
0436     if(isempty(maindir))
0437       maindir = getMainDir(d, field);
0438     end
0439     dbclear if error
0440     set(gcf,'paperposition',[0 0 6.0 9.0])
0441     filename = sprintf('%s/pointing/fig%d.png', maindir, plotIndex);
0442     eval(sprintf('print -dpng -r70 %s', filename));
0443     dbstop if error
0444       end
0445     end    
0446     clear azIs;
0447     clear elIs;
0448     
0449     pause(0.1);
0450   else
0451     display('All data in this cross have been flagged');
0452     azWidth(m,1:2) = nan;
0453     azMax(m,1:2)   = nan;
0454     azBase(m,1:2)  = nan;
0455     azRms(m,1:2)   = nan;
0456     betaAzAll(m,1:2,1:3) = nan;
0457     elWidth(m,1:2) = nan;
0458     elMax(m,1:2)   = nan;
0459     elBase(m,1:2)  = nan;
0460     elRms(m,1:2)   = nan;
0461     betaElAll(m,1:2,1:3) = nan;
0462     badPoint(m,1:2) = 1;
0463   end
0464 
0465 end
0466 
0467 
0468 badPoint = logical(badPoint);
0469 elOff = betaElAll(:,:,2);
0470 azOff = betaAzAll(:,:,2);
0471 elOff(badPoint) = nan;
0472 azOff(badPoint) = nan;
0473 indBad = badPoint;
0474 
0475 % first we assign a weight with the source.
0476 azSig = (azMax - azBase)./azRms;
0477 azWidth = abs(betaAzAll(:,:,3)/2);
0478 azErr = azWidth./azSig;
0479 
0480 
0481 elSig = (elMax - elBase)./elRms;
0482 elWidth = abs(betaElAll(:,:,3)/2);
0483 elErr   = elWidth./elSig;
0484 
0485 
0486 % throw out values where the error is large or negative.
0487 elBad = elErr < 0 | elErr > 1;
0488 azBad = azErr < 0 | azErr > 1;
0489 
0490 elOff(elBad) = nan;
0491 azOff(azBad) = nan;
0492 elErr(elBad) = nan;
0493 azErr(elBad) = nan;
0494 indBad(elBad) = 1;
0495 indBad(azBad) = 1;
0496 
0497 
0498 % if hte values for both channels aren't within 10 arcmin, toss them.
0499 dEl = abs(diff(elOff, [], 2));
0500 dAz = abs(diff(azOff, [], 2).*cos(obs.el*pi/180)');
0501 
0502 f = find(dEl>20/60 | dAz>20/60);
0503 elOff(f,:) = nan;
0504 azOff(f,:) = nan;
0505 elErr(f,:) = nan;
0506 azErr(f,:) = nan;
0507 indBad(f,:) = 1;
0508 
0509 badPoint = indBad;
0510 
0511 indBad = sum(indBad,2)==2;
0512 
0513 display(sprintf('Throwing out %d points', length(find(indBad))));
0514 
0515 [off.az  off.azErr] = vwstat(azOff, azErr, 2);
0516 [off.el  off.elErr] = vwstat(elOff, elErr, 2);
0517 
0518 display('Done with all crosses');
0519 
0520 % add final offests to the observed az/el
0521 obs.az = obs.az + off.az';
0522 obs.el = obs.el + off.el';
0523 
0524 obsel = obs.el;
0525 
0526 out.time   = timeVal;
0527 out.name   = sourceName;
0528 out.az     = obs.az;
0529 out.el     = obs.el;
0530 
0531 % throw out the data that is bad
0532 ind = isnan(obs.az);
0533 obs.az(ind) = [];
0534 obs.el(ind) = [];
0535 obs.name    = sourceName(~ind);
0536 obs.timeVal = timeVal(~ind);
0537 obs.index   = find(~ind);
0538 ide.az(ind) = [];
0539 ide.el(ind) = [];
0540 off.az(ind) = [];
0541 off.el(ind) = [];
0542 off.azErr(ind) = [];
0543 off.elErr(ind) = [];
0544 
0545 
0546 %azOff(ind,:) = [];
0547 %azErr(ind,:) = [];
0548 
0549 
0550 out.xwidth = nan(size(azOff));
0551 out.ywidth = nan(size(azOff));
0552 out.xoff   = azOff./sind([obsel' obsel']);
0553 out.yoff   = elOff;
0554 out.errpeak= nan(size(azErr));
0555 out.errxoff= azErr./sind([obsel' obsel']);
0556 out.erryoff= elErr;
0557 out.eloff  = elOff;
0558 out.azoff  = azOff;
0559 out.azmin  = azBase;
0560 out.azpeak = azMax;
0561 out.azrms  = azRms;
0562 out.elmin  = elBase;
0563 out.elpeak = elMax;
0564 out.elrms  = elRms;
0565 out.azsig  = azSig;
0566 out.elsig  = elSig;
0567 out.badpt  = badPoint | isnan(azOff);
0568 
0569 
0570 return;
0571 
0572 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0573 %function [s e] = findStartStop(features)
0574 %
0575 %% start indices are when the difference in the feature value increases, stop
0576 %% is when it decreases
0577 %s = find(diff(double(features))>0);
0578 %
0579 %% end is when it decreases
0580 %e = find(diff(double(features))<0);
0581 %
0582 %% last one
0583 %lastOneBad = 0;
0584 %if(length(e)<length(s))
0585 %  % the dimensions don't match
0586 %  if(s(1)<e(1))
0587 %    % if start is less than end, that means we're usually missing hte last
0588 %    % endpoint
0589 %    e(length(s)) = length(features);
0590 %    lastOneBad = 1;
0591 %  else
0592 %    s(1) = [];
0593 %  end
0594 %elseif(length(s)<length(e))
0595 %  s(1) = [];
0596 %end
0597 %s = s+1;
0598 %
0599 %goodPoint = zeros(size(s));
0600 %for m=1:length(s)
0601 %  thisFeat = unique(features(s(m):e(m)));
0602 %  if(thisFeat == 129 | thisFeat ==128)
0603 %    goodPoint(m) = 1;
0604 %  end
0605 %end
0606 %
0607 %if(lastOneBad==1)
0608 %  goodPoint(length(s)) = 0;
0609 %end
0610 %
0611 %
0612 %s = s(logical(goodPoint));
0613 %e = e(logical(goodPoint));
0614 %
0615 %return;
0616 
0617 
0618 function x = mchisq2(dataVals, modelVals)
0619 
0620 x = dataVals - modelVals;
0621 
0622 f = find(isnan(x));
0623 x(f) = [];
0624 f = find(isinf(x));
0625 x(f) = [];
0626 
0627 if(isempty(x))
0628   display('damn it all to hell')
0629   keyboard;
0630 end
0631 
0632 return;