apertures.m | searchcode

/ATF2/FlightSim/watchdogApps/apertures/apertures.m

http://atf2flightsim.googlecode.com/
MATLAB | 499 lines | 449 code | 7 blank | 43 comment | 77 complexity | c8613204055078bc0902fb67896a9b81 MD5 | raw file
Possible License(s): BSD-2-Clause, LGPL-2.0, IPL-1.0, BSD-3-Clause

function [stat output] = apertures(cmd,varargin)
% [stat output] = apertures(cmd)
% Watchdog function for watching orbit in apertures
%
% 100510 JLN change data.alarm=false to data.alarm=0, so this watchdog
%    is consistent with the others
persistent pars data lastupdate
stat{1}=1; output=[];

% Make sure initialised
if isempty(pars) && ~isequal(cmd,'init')
  [stat pars]=initPars; if stat{1}~=1; return; end;
  data=initData;
end

% Methods
switch lower(cmd)
  case 'init'
    [stat pars]=initPars; if stat{1}~=1; return; end;
    data=initData;
  case 'reset'
    [stat pars]=initPars('reset'); if stat{1}~=1; return; end;
    data=initData('reset');
  case 'isalarm'
    if isempty(lastupdate) || etime(clock,lastupdate)>pars.updateRate
      [stat data]=updateStatus(pars,data);
      if stat{1}~=1; return; end;
      lastupdate=clock;
    end
    output=data.alarm;
  case 'getpars'
    output=pars;
  case 'setpars'
    if nargin>0
      stat=checkPars(varargin{1});
      if stat{1}==1
        pars=setPars(pars,varargin{1});
      end
    end
  case 'update'
    pars.forceUpdate=true;
    [stat data]=updateStatus(pars,data);
    pars.forceUpdate=false;
    if stat{1}~=1; return; end;
    output=data;
%     data.alarm
  case 'getdata'
    output=data;
  case 'save'
    if isdir(fullfile('userData','guiData'))
      save(fullfile('userData','guiData','apertures.mat'),'pars','data');
    end
  case 'new'
    if nargin==3
      [stat pars data]=genNewAp(pars,data,varargin{1},varargin{2});
    end
  case 'bpmcheck'
    if nargin~=2
      error('need 2 input arguments')
    end
    [bpmind instrind]=bpmCheck(varargin{1});
    output=[bpmind; instrind];
  otherwise
    stat{1}=-1; stat{2}='Unknown command';
    return
end

%% Internal functions
function stat=checkPars(pars) %#ok<INUSD>
stat{1}=1;
function [stat pars data]=genNewAp(pars,data,bl_ind,name)
global BEAMLINE
stat{1}=1;
try
  pars.active=true;
  pars.(name).active=false;
  bpmind=findcells(BEAMLINE,'Class','MONI');
  if bl_ind<=0 || bl_ind>length(BEAMLINE)
    stat{1}=-1; stat{2}='Invalid beamline index';
    return
  end
  pars.(name).targetInd=bl_ind;
  pars.(name).bpmind=find(bpmind>pars.(name).targetInd,pars.nbpm);
  [pars.(name).bpmind pars.(name).bpminstr] =bpmCheck(pars.(name).bpmind);
  pars.(name).reforbit='default';
  pars.(name).shape='ellipse';
  pars.(name).xoff=0; pars.(name).yoff=0;
  pars.(name).xpoff=0; pars.(name).ypoff=0;
  pars.(name).xwid=1e-3; pars.(name).xpwid=1e-3;
  pars.(name).ywid=1e-3; pars.(name).ypwid=1e-3;
  pars.(name).xrot=0; pars.(name).yrot=0;
  pars=apCalc(pars,name);
  data.(name).alarm=0;
catch
  stat{1}=-1; stat{2}=sprintf('Error creating new aperture: %s',lasterr);
end
function pars=setPars(pars,parsin)
global INSTR
bpmind_change={}; bpminstr_change={}; apchange={};
apchangelist={'xoff' 'xpoff' 'yoff' 'ypoff' 'xwid' 'ywid' 'xpwid' 'ypwid' 'xrot' 'yrot' 'shape'};
% Update new pars fields
if isstruct(parsin)
  pf=fieldnames(parsin);
  for ipf=1:length(pf)
    if isstruct(parsin.(pf{ipf}))
      pf2=fieldnames(parsin.(pf{ipf}));
      for ipf2=1:length(pf2)
        if strcmp(pf2{ipf2},'bpmind')
          bpmind_change{end+1}=pf{ipf};
        elseif strcmp(pf2{ipf2},'bpminstr')
          bpminstr_change{end+1}=pf{ipf};
        elseif ismember(pf2{ipf2},apchangelist)
          apchange{end+1}=pf{ipf};
        end
        pars.(pf{ipf}).(pf2{ipf2})=parsin.(pf{ipf}).(pf2{ipf2});
      end
    else
      pars.(pf{ipf})=parsin.(pf{ipf});
    end
  end
end
% Check for changes in bpm indicies
if ~isempty(bpmind_change)
  for iap=1:length(bpmind_change)
    [pars.(bpmind_change{iap}).bpmind pars.(bpmind_change{iap}).bpminstr]=bpmCheck(pars.(bpmind_change{iap}).bpmind);
  end
end
if ~isempty(bpminstr_change)
  for iap=1:length(bpminstr_change)
    pars.(bpminstr_change{iap}).bpmind=arrayfun(@(x) INSTR{x}.Index,pars.(bpminstr_change{iap}).bpminstr);
    [pars.(bpminstr_change{iap}).bpmind pars.(bpmind_change{iap}).bpminstr]=bpmCheck(pars.(bpminstr_change{iap}).bpmind);
  end
end
% recalc apertures if needed
if ~isempty(apchange)
  for iap=1:length(apchange)
    pars=apCalc(pars,apchange{iap});
  end
end
function [stat pars]=initPars(varargin)
global BEAMLINE INSTR
stat{1}=1;
% Test for reset command
if nargin && isequal(varargin{1},'reset')
  reset=true;
else
  reset=false;
end
% Load pars from last session if existing
if isdir(fullfile('userData','watchdogData')) && exist(fullfile('userData','watchdogData','apertures.mat'),'file') && ~reset
  load(fullfile('userData','watchdogData','apertures.mat'),'pars');
else % initialise here
  pars.active=true;
  pars.guiActive=false;
  pars.bpmave=10;
  pars.nbpm=5;
  pars.maxrms=3;
  pars.minq=0.5;
  pars.updateRate=60; % s
  pars.whichICT='ICTDUMP';
  pars.forceUpdate=false;
  bh=findcells(BEAMLINE,'Name','BH*XB');
  bpmind=findcells(BEAMLINE,'Class','MONI');
  nohw=arrayfun(@(x) INSTR{x}.Index,findcells(INSTR,'Class','NOHW'));
  bpmind=bpmind(~ismember(bpmind,nohw));
  pars.bpmind=bpmind;
  % ============== Septum
  pars.Septum.active=false;
  pars.Septum.targetInd=findcells(BEAMLINE,'Name','BS1XB');
  pars.Septum.bpmind=bpmind(find(bpmind>pars.Septum.targetInd,pars.nbpm));
  [pars.Septum.bpmind pars.Septum.bpminstr]=bpmCheck(pars.Septum.bpmind);
  pars.Septum.reforbit='default';
  pars.Septum.shape='ellipse';
  pars.Septum.xoff=0; pars.Septum.yoff=0;
  pars.Septum.xpoff=0; pars.Septum.ypoff=0;
  pars.Septum.xwid=1e-3; pars.Septum.xpwid=1e-3;
  pars.Septum.ywid=1e-3; pars.Septum.ypwid=1e-3;
  pars.Septum.xrot=0; pars.Septum.yrot=0;
  pars=apCalc(pars,'Septum');
  % ============== BH1X
  pars.BH1X.active=false;
  pars.BH1X.targetInd=bh(1);
  pars.BH1X.bpmind=bpmind(find(bpmind>pars.BH1X.targetInd,pars.nbpm));
  [pars.BH1X.bpmind pars.BH1X.bpminstr]=bpmCheck(pars.BH1X.bpmind);
  pars.BH1X.reforbit='default';
  pars.BH1X.shape='ellipse';
  pars.BH1X.xoff=0; pars.BH1X.yoff=0;
  pars.BH1X.xpoff=0; pars.BH1X.ypoff=0;
  pars.BH1X.xwid=1e-3; pars.BH1X.xpwid=1e-3;
  pars.BH1X.ywid=1e-3; pars.BH1X.ypwid=1e-3;
  pars.BH1X.xrot=0; pars.BH1X.yrot=0;  
  pars=apCalc(pars,'BH1X');
  % ============== BH2X
  pars.BH2X.active=false;
  pars.BH2X.targetInd=bh(2);
  pars.BH2X.bpmind=bpmind(find(bpmind>pars.BH2X.targetInd,pars.nbpm));
  [pars.BH2X.bpmind pars.BH2X.bpminstr]=bpmCheck(pars.BH2X.bpmind);
  pars.BH2X.reforbit='default';
  pars.BH2X.shape='ellipse';
  pars.BH2X.xoff=0; pars.BH2X.yoff=0;
  pars.BH2X.xpoff=0; pars.BH2X.ypoff=0;
  pars.BH2X.xwid=1e-3; pars.BH2X.xpwid=1e-3;
  pars.BH2X.ywid=1e-3; pars.BH2X.ypwid=1e-3;
  pars.BH2X.xrot=0; pars.BH2X.yrot=0;
  pars=apCalc(pars,'BH2X');
  % ============== BH3X
  pars.BH3X.active=false;
  pars.BH3X.targetInd=bh(3);
  pars.BH3X.bpmind=bpmind(find(bpmind>pars.BH3X.targetInd,pars.nbpm));
  [pars.BH3X.bpmind pars.BH3X.bpminstr]=bpmCheck(pars.BH3X.bpmind);
  pars.BH3X.reforbit='default';
  pars.BH3X.shape='ellipse';
  pars.BH3X.xoff=0; pars.BH3X.yoff=0;
  pars.BH3X.xpoff=0; pars.BH3X.ypoff=0;
  pars.BH3X.xwid=1e-3; pars.BH3X.xpwid=1e-3;
  pars.BH3X.ywid=1e-3; pars.BH3X.ypwid=1e-3;
  pars.BH3X.xrot=0; pars.BH3X.yrot=0;
  pars=apCalc(pars,'BH3X');
  % ============= KEX2
  pars.KEX2.active=false;
  pars.KEX2.targetInd=findcells(BEAMLINE,'Name','IP02');
  pars.KEX2.bpmind=bpmind(find(bpmind>pars.KEX2.targetInd,pars.nbpm));
  [pars.KEX2.bpmind pars.KEX2.bpminstr]=bpmCheck(pars.KEX2.bpmind);
  pars.KEX2.reforbit='default';
  pars.KEX2.shape='ellipse';
  pars.KEX2.xoff=0; pars.KEX2.yoff=0;
  pars.KEX2.xpoff=0; pars.KEX2.ypoff=0;
  pars.KEX2.xwid=1e-3; pars.KEX2.xpwid=1e-3;
  pars.KEX2.ywid=1e-3; pars.KEX2.ypwid=1e-3;
  pars.KEX2.xrot=0; pars.KEX2.yrot=0;
  pars=apCalc(pars,'KEX2');
  % ============= Form response matricies
  [stat pars]=getRmat(pars); if stat{1}~=1; return; end;
end
function data=initData(varargin)
% Test for reset command
if nargin && isequal(varargin{1},'reset')
  reset=true;
else
  reset=false;
end
% Load data from last session if existing
if isdir(fullfile('userData','watchdogData')) && exist(fullfile('userData','watchdogData','apertures.mat'),'file') && ~reset
  load(fullfile('userData','watchdogData','apertures.mat'),'data');
else % initialise here
end
%
% 100510 JLN change all of these from false (logical) to 0 (double)
%
data.alarm=0;%false; % start with alarm state being false
data.Septum.alarm=0;%false;
data.BH1X.alarm=0;%false;
data.BH2X.alarm=0;%false;
data.BH3X.alarm=0;%false;
data.KEX2.alarm=0;%false;
    
function [stat data]=updateStatus(pars,data)
global INSTR
stat{1}=1; %#ok<NASGU>
% Update response matrix
[stat pars]=getRmat(pars); if stat{1}~=1; return; end;
% Get BPM readings
if strcmp(pars.whichICT,'ICTDUMP')
  w_ict=0;
elseif strcmp(pars.whichICT,'ICT1X')
  w_ict=1;
else
  stat{1}=-1; stat{2}='Invalid ICT choice'; return;
end
% Check there is enough bpm data
[stat bsize]=FlHwUpdate('buffersize');
if stat{1}~=1 || bsize<pars.bpmave
  stat{2}='Not enough buffered bpm data or not available';
  return
end
% Get averaged bpm readings
[stat bpmdata]=FlHwUpdate('bpmave',pars.bpmave,0,[pars.minq,pars.maxrms,w_ict,0]); if stat{1}~=1; return; end;
% Loop over all apertures
pf=fieldnames(pars); data.alarm=0;
for ipf=1:length(pf)
  if isfield(pars.(pf{ipf}),'targetInd') && (pars.(pf{ipf}).active || pars.forceUpdate)
    % Treat upstream bpms seperately
    us=pars.(pf{ipf}).bpmind<pars.(pf{ipf}).targetInd;
    ds=~us;
    % Pull out BPM data of interest and get difference from gold orbit
    try
      if isequal(pars.(pf{ipf}).reforbit,'default')
        bpmx_ref=arrayfun(@(x) INSTR{x}.ref(1),pars.(pf{ipf}).bpminstr(ds))';
        bpmsx_ref=arrayfun(@(x) INSTR{x}.referr(1),pars.(pf{ipf}).bpminstr(ds))';
        bpmy_ref=arrayfun(@(x) INSTR{x}.ref(2),pars.(pf{ipf}).bpminstr(ds))';
        bpmsy_ref=arrayfun(@(x) INSTR{x}.referr(2),pars.(pf{ipf}).bpminstr(ds))';
      else
        bpmx_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(ds),1);
        bpmy_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(ds),2);
        bpmsx_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(ds),4);
        bpmsy_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(ds),5);
      end
    catch
      stat{1}=-1; stat{2}='No reference orbit defined'; return;
    end
    if any(isnan(bpmx_ref)) || any(isnan(bpmy_ref))
      stat{1}=-1; stat{2}='bad reference orbit';
      return
    end
    bpmx=bpmdata{1}(1,pars.(pf{ipf}).bpminstr(ds))'-bpmx_ref; bpmsx=sqrt(bpmdata{1}(4,pars.(pf{ipf}).bpminstr(ds))'.^2+bpmsx_ref.^2);
    bpmy=bpmdata{1}(2,pars.(pf{ipf}).bpminstr(ds))'-bpmy_ref; bpmsy=sqrt(bpmdata{1}(5,pars.(pf{ipf}).bpminstr(ds))'.^2+bpmsy_ref.^2);
    if any(isnan(bpmx)) || any(isnan(bpmy))
      stat{1}=-1; stat{2}='bad bpm readings';
      return
    end
    % Store bpm data
    data.(pf{ipf}).bpm=[bpmx bpmy]; data.(pf{ipf}).bpm_error=[bpmsx bpmsy];
    % Get estimated orbit at aperture(s) including error estimate based on bpm
    % readings RMS
    [data.(pf{ipf}).pos,data.(pf{ipf}).error,mse]=lscov([pars.(pf{ipf}).Rx(ds,:);pars.(pf{ipf}).Ry(ds,:)],[bpmx;bpmy],1./[bpmsx;bpmsy].^2);
    data.(pf{ipf}).error=data.(pf{ipf}).error*sqrt(1/mse);
    % Update alarm (true if outside aperture) % major alarm if pos outside,
    % minor if error bar outside
    ax=data.(pf{ipf}).pos(1); apx=data.(pf{ipf}).pos(2);
    asx=data.(pf{ipf}).error(1); apsx=data.(pf{ipf}).error(2);
    ay=data.(pf{ipf}).pos(3); apy=data.(pf{ipf}).pos(4);
    asy=data.(pf{ipf}).error(3); apsy=data.(pf{ipf}).error(4);
    data.(pf{ipf}).alarm=0;
    if pars.(pf{ipf}).active 
      if ~inpolygon(ax,apx,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ay,apy,pars.(pf{ipf}).aper_y,pars.(pf{ipf}).aper_yp)
        data.(pf{ipf}).alarm=2;
        data.alarm=2;
      elseif ~inpolygon(ax+asx,apx+apsx,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ax-asx,apx-apsx,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ax-asx,apx+apsx,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ax+asx,apx-apsx,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ay+asy,apy+apsy,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ay-asy,apy-apsy,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ay+asy,apy-apsy,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) || ...
          ~inpolygon(ay-asy,apy+apsy,pars.(pf{ipf}).aper_x,pars.(pf{ipf}).aper_xp) && data.(pf{ipf}).alarm~=2
        data.(pf{ipf}).alarm=1;
        if data.alarm~=2; data.alarm=1; end;
      end
    else
      data.alarm=0;
    end
    % Deal with upstream fitting
    if any(us)
      dsbpm=find(pars.(pf{ipf}).bpmind<pars.(pf{ipf}).targetInd, 1, 'last' );
      if isempty(dsbpm)
        stat{1}=-1; stat{2}='Must be one downstream bpm provided!';
        return
      end      
      try
        if isequal(pars.(pf{ipf}).reforbit,'default')
          bpmx_ref=arrayfun(@(x) INSTR{x}.ref(1),pars.(pf{ipf}).bpminstr(us))';
          bpmsx_ref=arrayfun(@(x) INSTR{x}.referr(1),pars.(pf{ipf}).bpminstr(us))';
          bpmy_ref=arrayfun(@(x) INSTR{x}.ref(2),pars.(pf{ipf}).bpminstr(us))';
          bpmsy_ref=arrayfun(@(x) INSTR{x}.referr(2),pars.(pf{ipf}).bpminstr(us))';
          dsbpmx_ref=arrayfun(@(x) INSTR{x}.ref(1),pars.(pf{ipf}).bpminstr(dsbpm))';
          dsbpmsx_ref=arrayfun(@(x) INSTR{x}.referr(1),pars.(pf{ipf}).bpminstr(dsbpm))';
          dsbpmy_ref=arrayfun(@(x) INSTR{x}.ref(2),pars.(pf{ipf}).bpminstr(dsbpm))';
          dsbpmsy_ref=arrayfun(@(x) INSTR{x}.referr(2),pars.(pf{ipf}).bpminstr(dsbpm))';
        else
          bpmx_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(us),1);
          bpmsx_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(us),4);
          bpmy_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(us),2);
          bpmsy_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(us),5);
          dsbpmx_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(dsbpm),1);
          dsbpmsx_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(dsbpm),4);
          dsbpmy_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(dsbpm),2);
          dsbpmsy_ref=pars.(pf{ipf}).reforbit(pars.(pf{ipf}).bpminstr(dsbpm),5);
        end
      catch
        stat{1}=-1; stat{2}='No reference orbit defined'; return;
       end
      bpmx=bpmdata{1}(1,pars.(pf{ipf}).bpminstr(us))'-bpmx_ref; bpmsx=sqrt(bpmdata{1}(4,pars.(pf{ipf}).bpminstr(us))'.^2+bpmsx_ref.^2);
      bpmy=bpmdata{1}(2,pars.(pf{ipf}).bpminstr(us))'-bpmy_ref; bpmsy=sqrt(bpmdata{1}(5,pars.(pf{ipf}).bpminstr(us))'.^2+bpmsy_ref.^2);
      dsbpmx=bpmdata{1}(1,pars.(pf{ipf}).bpminstr(dsbpm))'-dsbpmx_ref; dsbpmsx=sqrt(bpmdata{1}(4,pars.(pf{ipf}).bpminstr(dsbpm))'.^2+dsbpmsx_ref.^2);
      dsbpmy=bpmdata{1}(2,pars.(pf{ipf}).bpminstr(dsbpm))'-dsbpmy_ref; dsbpmsy=sqrt(bpmdata{1}(5,pars.(pf{ipf}).bpminstr(dsbpm))'.^2+dsbpmsy_ref.^2);
      % Store bpm data
      data.(pf{ipf}).bpm_us=[bpmx bpmy]; data.(pf{ipf}).bpm_error_us=[bpmsx bpmsy];
      % Get estimated orbit at aperture(s) including error estimate
      [data.(pf{ipf}).pos_us,data.(pf{ipf}).error_us]=fwdprop(bpmx,bpmsx,bpmy,bpmsy,dsbpmx,dsbpmsx,dsbpmy,dsbpmsy,pars.(pf{ipf}).Rij,pars.(pf{ipf}).Rik);
    end
  elseif isfield(pars.(pf{ipf}),'targetInd')
    data.(pf{ipf}).alarm=0;
  end
end
function [Xj,sigmaXj]=fwdprop(xi,sxi,yi,syi,xk,sxk,yk,syk,Rij,Rik)
N = length(xi); Xj=zeros(4,1); sigmaXj=zeros(4,1);
for i=1:length(xi)
  a=Rik{i}(1,2)-(Rik{i}(1,4)*Rik{i}(3,2))/Rik{i}(3,4);
  b=Rik{i}(3,4)-(Rik{i}(3,2)*Rik{i}(1,4))/Rik{i}(1,2);
  c=(Rik{i}(3,1)*Rik{i}(1,4))/Rik{i}(3,4)-Rik{i}(1,1);
  d=(Rik{i}(1,1)*Rik{i}(3,2))/Rik{i}(1,2)-Rik{i}(3,1);
  e=(Rik{i}(1,4)*Rik{i}(3,3))/Rik{i}(3,4)-Rik{i}(1,3);
  f=(Rik{i}(1,3)*Rik{i}(3,2))/Rik{i}(1,2)-Rik{i}(3,3);
  g=Rik{i}(3,2)/(b*Rik{i}(1,2));
  h=Rik{i}(1,4)/(a*Rik{i}(3,4));
  for i_ind=1:4
    Mi(i_ind,1)=Rij{i}(i_ind,1)+(Rij{i}(i_ind,2)/a)*c+(Rij{i}(i_ind,4)/b)*d;
    Mi(i_ind,2)=Rij{i}(i_ind,3)+(Rij{i}(i_ind,2)/a)*e+(Rij{i}(i_ind,4)/b)*f;
    Mi(i_ind,3)=Rij{i}(i_ind,2)/a - Rij{i}(i_ind,4)*g;
    Mi(i_ind,4)=Rij{i}(i_ind,4)/b - Rij{i}(i_ind,2)*h;
  end
  Xj = Xj + Mi * [xi(i);yi(i);xk;yk] ;
  sigmaXj = sigmaXj + Mi.^2 * [sxi(i);syi(i);sxk;syk].^2 ;
end
Xj=Xj/N;
sigmaXj = sqrt(sigmaXj/N);

function [stat pars]=getRmat(pars)

stat{1}=1;
pf=fieldnames(pars);
for ipf=1:length(pf)
  if isstruct(pars.(pf{ipf})) && isfield(pars.(pf{ipf}),'targetInd') && (pars.(pf{ipf}).active || pars.forceUpdate)
    % Using specific list?
    if ~isfield(pars.(pf{ipf}),'bpmind')
      error('bpmind not defined for aperture: %s',pf{ipf})
    end
    pars.(pf{ipf}).Rx=zeros(length(pars.(pf{ipf}).bpmind),4);
    pars.(pf{ipf}).Ry=zeros(length(pars.(pf{ipf}).bpmind),4);
    for ibpm=1:length(pars.(pf{ipf}).bpmind)
      if pars.(pf{ipf}).targetInd<pars.(pf{ipf}).bpmind(ibpm)
        [stat,R]=RmatAtoB(pars.(pf{ipf}).targetInd,pars.(pf{ipf}).bpmind(ibpm));
        if stat{1}~=1
          return
        end
        pars.(pf{ipf}).Rx(ibpm,:)=[R(1,1) R(1,2) R(1,3) R(1,4)];
        pars.(pf{ipf}).Ry(ibpm,:)=[R(3,1) R(3,2) R(3,3) R(3,4)];
      else
        dsbpm=pars.(pf{ipf}).bpmind(find(pars.(pf{ipf}).bpmind>pars.(pf{ipf}).targetInd,1));
        if isempty(dsbpm)
          stat{1}=-1; stat{2}='Must be one downstream bpm provided!';
          return
        end
        [stat,R]=RmatAtoB(pars.(pf{ipf}).bpmind(ibpm),pars.(pf{ipf}).targetInd);
        if stat{1}~=1
          return
        end
        [stat,R2]=RmatAtoB(pars.(pf{ipf}).bpmind(ibpm),dsbpm);
        if stat{1}~=1
          return
        end     
        pars.(pf{ipf}).Rij{ibpm}=R;
        pars.(pf{ipf}).Rik{ibpm}=R2;
      end
    end
  end
end
function pars=apCalc(pars,apName)
pars.(apName).aper_x=[];
pars.(apName).aper_xp=[];
pars.(apName).aper_y=[];
pars.(apName).aper_yp=[];
if strcmp(pars.(apName).shape,'ellipse')
  a1=1/pars.(apName).xwid^2; a2=1/pars.(apName).ywid^2;
  b1=1/pars.(apName).xpwid^2; b2=1/pars.(apName).ypwid^2;
  c1=-(pars.(apName).xrot/45)*sqrt(a1*b1)/2;
  c2=-(pars.(apName).yrot/45)*sqrt(a2*b2)/2;
  [pars.(apName).aper_x,pars.(apName).aper_xp]=ellipse([a1 c1;c1 b1],0,0);
  [pars.(apName).aper_y,pars.(apName).aper_yp]=ellipse([a2 c2;c2 b2],0,0);
elseif strcmp(pars.(apName).shape,'rectangle')
  pars.(apName).aper_x=[linspace(-pars.(apName).xwid,pars.(apName).xwid,1000) ones(1,1000).*pars.(apName).xwid ...
    linspace(pars.(apName).xwid,-pars.(apName).xwid,1000) -ones(1,1000).*pars.(apName).xwid];
  pars.(apName).aper_xp=[ones(1,1000).*pars.(apName).xpwid linspace(pars.(apName).xpwid,-pars.(apName).xpwid,1000) ...
    ones(1,1000).*-pars.(apName).xpwid linspace(-pars.(apName).xpwid,pars.(apName).xpwid,1000)];
  pars.(apName).aper_y=[linspace(-pars.(apName).ywid,pars.(apName).ywid,1000) ones(1,1000).*pars.(apName).ywid ...
    linspace(pars.(apName).ywid,-pars.(apName).ywid,1000) -ones(1,1000).*pars.(apName).ywid];
  pars.(apName).aper_yp=[ones(1,1000).*pars.(apName).ypwid linspace(pars.(apName).ypwid,-pars.(apName).ypwid,1000) ...
    ones(1,1000).*-pars.(apName).ypwid linspace(-pars.(apName).ypwid,pars.(apName).ypwid,1000)];  
  xang=(pars.(apName).xrot/180)*pi; yang=(pars.(apName).yrot/180)*pi;
  ROTx=[cos(xang) -sin(xang); sin(xang) cos(xang)];
  ROTy=[cos(yang) -sin(yang); sin(yang) cos(yang)];
  newAper_x=ROTx*[pars.(apName).aper_x; pars.(apName).aper_xp];
  pars.(apName).aper_x=newAper_x(1,:); pars.(apName).aper_xp=newAper_x(2,:);
  newAper_y=ROTy*[pars.(apName).aper_y; pars.(apName).aper_yp];
  pars.(apName).aper_y=newAper_y(1,:); pars.(apName).aper_yp=newAper_y(2,:);
end
pars.(apName).aper_x=pars.(apName).aper_x+pars.(apName).xoff;
pars.(apName).aper_xp=pars.(apName).aper_xp+pars.(apName).xpoff;
pars.(apName).aper_y=pars.(apName).aper_y+pars.(apName).yoff;
pars.(apName).aper_yp=pars.(apName).aper_yp+pars.(apName).ypoff;

function [bpmind instrind]=bpmCheck(bpmind)
global INSTR
% check HW exists and want to use bpm list
% get bpm info from BPM Tool
instrind=[];
[stat data] = FlBpmToolFn('GetData');
if stat{1}~=1
  return
end
goodbpm=data.global.hw & data.global.use;
goodinst=ismember(cellfun(@(x) x.Index,INSTR),bpmind);
instrind=find(goodbpm & goodinst);
if isempty(instrind)
  bpmind=[];
else
  bpmind=arrayfun(@(x) INSTR{x}.Index,instrind);
end