clear all;close all;
h=figure('Position',[289 205 1205 695]);hold on;
colors = [      0    0.4470    0.7410
    0.8500    0.3250    0.0980
    0.4660    0.6740    0.1880
    0.4940    0.1840    0.5560
    0.9290    0.6940    0.1250
    0.3010    0.7450    0.9330
    0.6350    0.0780    0.1840];

% TP Definition
B0=0.972;
E0=18e5;
Lb = 50e-3;
Oe = Lb+24e-3;
Delta = 7e-3;
Ddet = 493e-3;
Te=Ddet-Oe;
Db=Ddet-Lb;
ymax = 70e-3;
%Species
qe = 1.602e-19;
mp = 1.6726e-27;
Z = [1 6 5 4];
A = [1 12 12 12];
Eminvis = [5.0560 1.2640*12 0.8778*12 0.5618*12]*1e6*qe;
Efilter = [0.85 14.16 14.16 14.16]*1e6*qe;
EpN = 0.001:0.0005:100;

%Final plate design
le = 30; Le = 300; h0 = 25;  hp = 50.3; h1=46;%h1 = 46;
y = 0:h0; x = Le*ones(size(y)); plot(x,y,'--r','LineWidth',1.5);
alpha = abs((h1-h0)/(Le-le)); beta = (h1*(Le+Oe)-h0*(le+Oe))/(Le-le);
x = 0:Le; y = beta-alpha*x; plot(x,y,'--r','LineWidth',1.5);
hF=fill([x fliplr(x)],[y zeros(size(y))],'r','EdgeColor','none'); set(hF,'facealpha',.15);
hF=fill([x fliplr(x)],[hp*ones(size(x)) fliplr(y)],'k','EdgeColor','none'); set(hF,'facealpha',.08);
% % Show original plates
x=0:300; y=hp*ones(size(x));
plot(x,y,'--k','LineWidth',1.5);
y = 0:hp; x = 300*ones(size(y));plot(x,y,'--k','LineWidth',1.5);

%Traces
clc
for mm = 1:length(Z)
%     disp(sprintf('%d/%d', Z(mm),A(mm)));
    E = EpN*1e6*qe*A(mm);
    le = sqrt(4*E*Delta/(Z(mm)*qe*E0));
    y = Z(mm)*qe*Lb*(Oe+le-Lb/2)./sqrt(2*A(mm)*mp*E);
    mClip = find(y<50.3e-3); mClip = [min(mClip) max(mClip)];
%     disp('Rectangular Clip (MeV)')
%     disp(EpN(mClip(1))*A(mm))
    mClip2 = find(y<beta*1e-3-alpha*le); mClip22=max(mClip2);mClip2 = min(mClip2); 
%     disp('Trapezoidal Clip (MeV)')
%     disp([EpN(mClip2)*A(mm) EpN(mClip22)*A(mm)]);
    mClip3 = find(le<=300e-3); mClip3 = max(mClip3);
%     disp('Plate exit (MeV)')
%     disp(EpN(mClip3)*A(mm))
    y = y(le<=300);
    le = le(le<=300);
    handleLines(mm)=plot(le*1e3,y*1e3,'Color',colors(mm,:),'LineWidth',1.5);
    plot(le(mClip)*1e3,y(mClip)*1e3,'o','MarkerEdgeColor',colors(mm,:),'MarkerSize',16,'LineWidth',0.5,'MarkerFaceColor', colors(mm,:));
    plot(le(mClip2)*1e3,y(mClip2)*1e3,'o','MarkerEdgeColor',colors(mm,:),'MarkerSize',16,'LineWidth',0.5,'MarkerFaceColor', colors(mm,:));
    plot(le(mClip22)*1e3,y(mClip22)*1e3,'s','MarkerEdgeColor',colors(mm,:),'MarkerSize',18,'LineWidth',0.5,'MarkerFaceColor', colors(mm,:));
    plot(le(mClip3)*1e3,y(mClip3)*1e3,'s','MarkerEdgeColor',colors(mm,:),'MarkerSize',18,'LineWidth',0.5,'MarkerFaceColor', colors(mm,:));
   
    leOld = le; yOld = y;
    
    E = (Z(mm)*qe*B0*Lb*(Ddet-Lb/2)/(sqrt(2*A(mm)*mp)*ymax))^2;
    le = sqrt(4*E*Delta/(Z(mm)*qe*E0));
    y = Z(mm)*qe*Lb*(Oe+le-Lb/2)./sqrt(2*A(mm)*mp*E);
    plot(le*1e3,y*1e3,'d','MarkerEdgeColor',colors(mm,:),'MarkerFaceColor', colors(mm,:),'MarkerSize',18,'LineWidth',0.5);
    plot(le*1e3,y*1e3,'d','MarkerEdgeColor',colors(mm,:)*0.5,'MarkerSize',18,'LineWidth',1.5);
    
    le = sqrt(4*Efilter(mm)*Delta/(Z(mm)*qe*E0));
    y = Z(mm)*qe*Lb*(Oe+le-Lb/2)./sqrt(2*A(mm)*mp*Efilter(mm));
    plot(le*1e3,y*1e3,'p','MarkerEdgeColor',colors(mm,:),'MarkerFaceColor',colors(mm,:),'MarkerSize',20,'LineWidth',0.5);
    plot(le*1e3,y*1e3,'p','MarkerEdgeColor',colors(mm,:)*0.5,'MarkerSize',20,'LineWidth',1.5);
    
    le = leOld; y=yOld;
    plot(le(mClip)*1e3,y(mClip)*1e3,'o','MarkerEdgeColor',colors(mm,:)*.5,'MarkerSize',16,'LineWidth',1.5);
    plot(le(mClip2)*1e3,y(mClip2)*1e3,'o','MarkerEdgeColor',colors(mm,:)*.5,'MarkerSize',16,'LineWidth',1.5);
    plot(le(mClip22)*1e3,y(mClip22)*1e3,'s','MarkerEdgeColor',colors(mm,:)*.4,'MarkerSize',18,'LineWidth',1.5);
    plot(le(mClip3)*1e3,y(mClip3)*1e3,'s','MarkerEdgeColor',colors(mm,:)*.5,'MarkerSize',18,'LineWidth',1.5);
end

daspect([3.5,1,1]);
ylim([20 70]);  xlim([0 300]);
set(gca,'FontSize',22);
legend(handleLines,{'p^+','C^{6+}','C^{5+}','C^{4+}'},'Location', 'SouthWest');
hsp1 = get(gca, 'OuterPosition');

%%
h2=figure();hold on;
daspect([1,1,1]);
Emaxy = [5.0560 1.2640 0.8778 0.5618];
Efilter = [0.85 14.16/12 14.16/12 14.16/12];

Delta = 7e-3;
le = 30e-3; Le = 300e-3; h0 = 25e-3; h1 = 46e-3;
alpha = abs((h1-h0)/(Le-le)); beta = (h1*(Le+Oe)-h0*(le+Oe))/(Le-le);
grey=0.7*ones(1,3);

for mm = 1:length(Z) %unclipped plates traces
    EpN=max([Emaxy(mm) Efilter(mm)]):0.005:25;
    E = EpN*1e6*qe*A(mm);
    Eclip = Z(mm)*qe*E0*Le^2./(4*Delta);
    E=E(E>Eclip);
    x = Z(mm)*qe*E0*Le*(Te-Le/2)./(2*E);
    y = Z(mm)*qe*B0*Lb*(Db+Lb/2)./sqrt(2*A(mm)*mp*E);
    handleLines(mm)=plot(x*1e3,y*1e3,'LineWidth',1.5,'Color',colors(mm,:)); 
    E = EpN*1e6*qe*A(mm);
    E=E(E<Eclip);
    x = Z(mm)*qe*E0*Le*(Te-Le/2)./(2*E);
    y = Z(mm)*qe*B0*Lb*(Db+Lb/2)./sqrt(2*A(mm)*mp*E);
    plot(x*1e3,y*1e3,'--','Color',colors(mm,:),'LineWidth',1.5); 
end
plot(0,0,'ok','MarkerSize',8);
plot(0,0,'.k','MarkerSize',8);
xlim([-1 30]);  ylim([-1 70]);
xlabel('x (mm)','FontWeight','bold');
ylabel('y (mm)','FontWeight','bold');
legend(handleLines,{'p^+','C^{6+}','C^{5+}','C^{4+}'},'Location','southeast');
set(gca,'FontSize',18);
%%
h3=figure();hold on;
EpN=0.3:0.01:50;
zInPlates = 0:0.001:Le; yInPlates = h1*ones(size(zInPlates)); yInPlates (zInPlates>le)= -alpha*(zInPlates(zInPlates>le)+Oe)+beta;
for mm = 1:length(Z) %clipped plates traces
    EpN=max([Emaxy(mm) Efilter(mm)]):0.005:30;
    leClipped=[];
    for nn=1:length(EpN)
        E = EpN(nn)*1e6*qe*A(mm);
        y = Z(mm)*qe*B0*Lb*(zInPlates+Oe-Lb/2)./sqrt(2*A(mm)*mp*E);
        temp=zInPlates(y>yInPlates);
        if ~isempty(temp)
            leClipped(nn)=temp(1);
        else
            leClipped(nn)=Le;
        end
    end
    E = EpN*1e6*qe*A(mm);
    Eclip = Z(mm)*qe*E0*leClipped.^2./(4*Delta);
    min(E(E>Eclip))/(qe*1e6*A(mm));
    x = Z(mm)*qe*E0.*leClipped(E>Eclip).*(Te-leClipped(E>Eclip)/2)./(2*E(E>Eclip));
    y = Z(mm)*qe*B0*Lb*(Db+Lb/2)./sqrt(2*A(mm)*mp*E(E>Eclip));
    handleLines(mm)=plot(x*1e3,y*1e3,'LineWidth',1.5,'Color',colors(mm,:)); hold on;
    
    E = EpN*1e6*qe*A(mm);
    x = Z(mm)*qe*E0.*leClipped(E<Eclip).*(Te-leClipped(E<Eclip)/2)./(2*E(E<Eclip));
    y = Z(mm)*qe*B0*Lb*(Db+Lb/2)./sqrt(2*A(mm)*mp*E(E<Eclip));
    plot(x*1e3,y*1e3,'Color',colors(mm,:),'LineWidth',1.5); 
end

plot(0,0,'ok','MarkerSize',8);
plot(0,0,'.k','MarkerSize',8);
xlim([-1 30]);  ylim([-1 70]);
daspect([1,1,1]);

xlabel('x (mm)','FontWeight','bold');
ylabel('y (mm)','FontWeight','bold');

set(gca,'FontSize',18);