Scan2x2DinCOD

[bestinputcoord]=ScanPosAngle(...
     ropen,... lattice 
     inCOD,... initial trajectory [x,xp,y,yp]
     ngridstep,... grid steps in 1D (51 by default)
     gridsize,... grid size (5mm by default)
     figlabel)     figure label name. if not given, no figure

 scans initial coordinates for maximum number of turns

see also: ringpass

0001 function [bestinputcoord]=Scan2x2DinCOD(ropen,inCOD,ngridstep,gridsize,figlabel)
0002 %[bestinputcoord]=ScanPosAngle(...
0003 %     ropen,... lattice
0004 %     inCOD,... initial trajectory [x,xp,y,yp]
0005 %     ngridstep,... grid steps in 1D (51 by default)
0006 %     gridsize,... grid size (5mm by default)
0007 %     figlabel)     figure label name. if not given, no figure
0008 %
0009 % scans initial coordinates for maximum number of turns
0010 %
0011 %see also: ringpass
0012 
0013 if nargin<5
0014 figlabel=[];
0015 end
0016 if nargin<4
0017     gridsize=5e-3;
0018 end
0019 if nargin<3
0020     ngridstep=51;
0021 end
0022 
0023 % maxexcurs=20e-3;
0024 % Xapert=maxexcurs*ones(size(ropen));
0025 % Yapert=maxexcurs*ones(size(ropen));
0026 % ropen=SetPhysicalAperture(ropen,Xapert/2,Yapert/2);
0027 
0028 bestinputcoord=inCOD*0;
0029 TotTurn=10;
0030 
0031 DX=linspace(-1,1,ngridstep)*gridsize;
0032 DXP=DX/5;
0033 DY=DX;
0034 DYP=DX;
0035 
0036 % initial coordinate matrix
0037 
0038 [x,xp]=meshgrid(DX,DXP);
0039 [y,yp]=meshgrid(DY,DYP);
0040 
0041 incoory=repmat(inCOD,1,length(x(:)))+...
0042     [repmat(bestinputcoord([1]),length(x(:)),1),...
0043     repmat(bestinputcoord([2]),length(x(:)),1),...
0044     y(:),yp(:),zeros(size(xp(:),1),2)]';
0045 
0046 [~,lost,NT,lossinfo]=ringpass(ropen,incoory,TotTurn);
0047 
0048 lossinfo.element(isnan(lossinfo.element))=0;% particle not lost!
0049 NT(lost==0)=TotTurn;
0050 maxturnsy=(max(NT)+1);
0051 elpasy=lossinfo.element+length(ropen).*NT;
0052 
0053 [~,indt]=max(elpasy);
0054 % if more then one point, take closest to center.
0055 guessbest=sort([indt find(lost==0)]);
0056 
0057 dist=sqrt(incoory([3],guessbest).^2+incoory([4],guessbest).^2);
0058 [~,indb]=min(dist);
0059 
0060 bestinputcoord([3,4])=incoory([3,4],guessbest(indb));
0061 
0062 
0063 incoorx=repmat(inCOD,1,length(x(:)))+...
0064     [x(:),xp(:),...
0065     repmat(bestinputcoord([3]),length(x(:)),1),...
0066     repmat(bestinputcoord([4]),length(x(:)),1),...
0067     zeros(size(xp(:),1),2)]';
0068 
0069 [~,lost,NT,lossinfo]=ringpass(ropen,incoorx,TotTurn);
0070 lossinfo.element(isnan(lossinfo.element))=0;% particle not lost!
0071 NT(lost==0)=TotTurn;
0072 maxturnsx=(max(NT)+1);
0073 
0074 elpasx=lossinfo.element+length(ropen).*NT;
0075 
0076 [~,indt]=max(elpasx);
0077 guessbest=sort([indt find(lost==0)]);
0078 
0079 dist=sqrt(incoorx([1],guessbest).^2+incoorx([2],guessbest).^2);
0080 [~,indb]=min(dist);
0081 
0082 bestinputcoord([1,2])=incoorx([1,2],guessbest(indb));
0083 
0084 
0085 %%
0086 if ~isempty(figlabel)
0087     maxturns=max([maxturnsx,maxturnsy]);
0088     
0089     figure('Units','normalized','Position',[0.1 0.2 0.6 0.5]);
0090     subplot(1,2,1);
0091     surf(x,xp,reshape(elpasx,length(DX),[]));
0092       hcb=colorbar;
0093     caxis([0,length(ropen)*maxturns]);
0094     set(hcb,'YTick',[1:maxturns]*length(ropen));
0095     set(hcb,'YTickLabel',arrayfun(@(a)num2str(a,'turn %d'),[1:maxturns],'un',0));
0096     xlabel('x');ylabel('xp');view(2);shading flat;
0097     %title({['Number of turns vs injection position'],['Initial: ' num2str(inCOD)],['Best: ' num2str(bestinputcoord) ]});
0098     title({[],[],[],[],[]})
0099     subplot(1,2,2);
0100     surf(y,yp,reshape(elpasy,length(DY),[]));
0101     xlabel('y');ylabel('yp');view(2);shading flat;
0102      hcb=colorbar;
0103     caxis([0,length(ropen)*maxturns]);
0104     set(hcb,'YTick',[1:maxturns]*length(ropen));
0105     set(hcb,'YTickLabel',arrayfun(@(a)num2str(a,'turn %d'),[1:maxturns],'un',0));
0106     title({[],[],[],[],[]})
0107     
0108     set(gcf,'NextPlot','add');
0109     axes;
0110     h = title({['Number of turns vs injection position'],['Initial: ' num2str(inCOD)],['Best: ' num2str(bestinputcoord) ]});
0111     set(gca,'Visible','off');
0112     set(h,'Visible','on');
0113 
0114     
0115     saveas(gca,[figlabel '.fig']);
0116     try
0117         export_fig([figlabel '.jpg']);
0118     catch
0119     disp('missing export_fig, no jpg')
0120     end
0121     
0122 end
0123 
0124 return