% Probability of Intercept Simulation
% ------------------------------------

clear;clf;clc;

% Input Window Parameters

s=[];
T1=6.13;    	% Emitter 1 Period in s	
T2=1.05;    	% Emitter 2 Period in s
t1=.0333;   	% Intercept 1 Pulse in s
t2=.00278;  	% Intercept 2 Pulse in s
df1=100*t1/T1;	% Duty Factor 1 - %
df2=100*t2/T2;	% Duty Factor 2 - %
n=1000;        % Number of Trials

% Simulate Probability of Intercept

for i=1:n;
p=[];
c=0;
r=rand;
while length(p)<1,
   t=linspace(c*T2,c*T2+t2,100);
   x=square(2*pi*(t/T1+r),df1);
   z=(x+1)/2;
   p=find(z);
   tot=c*T2;
   if tot>10000,break,end;
   c=c+1;
end;

% Find Probability of Intercept

s(i)=c*T2;
phase(i)=r;
end;
[hits,time]=hist(s,600);
cumulative=cumsum(hits);
y=cumulative/n;

% Plot Probability of Intercept


x1=time;
stairs(x1,y);
axis([0,500,0,1]);
grid;
title('Probability of Main Beam Intercept');
xlabel('Intercept Duration-seconds');
ylabel('Cumulative Probability of Intercept');
text(350,.5,['T1=',num2str(T1),'sec']);
text(350,.45,['t1=',num2str(t1),' sec']);
text(350,.4,['T2=',num2str(T2),' sec']);
text(350,.35,['t2=',num2str(t2),' sec']);
hold on;

% Calculate Probability of Intercept

POI95=.95;
To=T1*T2/(t1+t2);	      % Mean Intercept Time
Tj=-To*log(1-POI95);    % Time POI=.95
tx=0:1:Tj;
POI=1-exp(-tx/To);
plot(tx,POI);
hold off;