function cp = coupling(xm, ym, relationship)
    mapsize = size(xm);
    switch relationship
        case 'h2v'
            midcol = ceil(mapsize(2)/2);
            titi = [ym(:,1:midcol-1)-ym(:,2:midcol), zeros(mapsize(1),1), ym(:,midcol+1:end)-ym(:,midcol:end-1)];
            toto = [xm(:,1:midcol-1)-xm(:,2:midcol), zeros(mapsize(1),1), xm(:,midcol+1:end)-xm(:,midcol:end-1)];
        case 'v2h'
            midrow = ceil(mapsize(1)/2);
            titi = [xm(1:midrow-1,:)-xm(2:midrow,:); zeros(1,mapsize(2)); xm(midrow+1:end,:)-xm(midrow:end-1,:)]; 
            toto = [ym(1:midrow-1,:)-ym(2:midrow,:); zeros(1,mapsize(2)); ym(midrow+1:end,:)-ym(midrow:end-1,:)]; 
    end
    cp = abs(100 * titi ./ toto );
    cp(isnan(cp)) = 0;
end


% ---- Andriy's initial code, April 2018 ----
% - function call -
%dosCoupH2V = coupling(yr, ky*y_dos_grid, bpm.map_dos_x(2));
%dosCoupV2H = coupling(xr, kx*x_dos_grid, bpm.map_dos_y(2));

% - function -
%function cp = coupling(xr,xm,stepsize)
%cp = abs(100 * (xr-xm)./stepsize); 




%couplingH2V = couplingKBS(kx*x_dos_grid, ky*y_dos_grid, 'h2v');
%couplingV2H = couplingKBS(kx*x_dos_grid, ky*y_dos_grid, 'v2h');

% ==== Kees on 18.06.2018 ==== 
%titi = (simset.ky_manual*simset.y_dos_grid(:,2:61)-simset.ky_manual*simset.y_dos_grid(:,1:60));
%toto = (simset.kx_manual*simset.x_dos_grid(:,2:61)-simset.kx_manual*simset.x_dos_grid(:,1:60));
%CH2V=abs(100*titi./toto);
%image([-15:0.5:14.5],[-8:0.5:7.5],CH2V);
%colormap(jet(100));
%colorbar
%axis([-10 10 -10 10]);

% function cp = couplingKBS(xm, ym, relationship)
%     switch relationship
%         case 'h2v'
%             titi = ym(:,2:end) - ym(:,1:end-1); 
%             toto = xm(:,2:end) - xm(:,1:end-1); 
%             cp = abs(100 * titi ./ toto );
%         case 'v2h'
%             titi = ym(2:end,:) - ym(1:end-1,:); 
%             toto = xm(2:end,:) - xm(1:end-1,:); 
%             cp = abs(100 * toto ./ titi );
%     end
%     
% end



% ==== Kees on 16.04.2018 ====
%at any given point in the BPM geometry, move the beam with a step (e.g. 100um) in one plane only,
%and asses the difference between the 2 corresponding values in the other plane 
%(as provided by whatever calculation method, DoS or Polyn.).
%Coupling [%] is 100x that difference divided by step (e.g. 100um).
%couplingAN(yr, ky*y_dos_grid, bpm.map_dos_x(2));
%couplingAN(xr, kx*x_dos_grid, bpm.map_dos_y(2));

% function cp = couplingAN(xr,xm,stepsize)
% 
%     cp = abs(100 * (xr-xm)./stepsize); 
% end
