clc
close all
clear
tlength=300*100;
clength=300*50;
datalength=tlength-clength;
feedbacklevel=0.008;
alpha=6; T=100; P=5.2; cp=10; tao=10;
global k x c v b n;
x=alpha; c=T; v=P; b=cp; n=tao;
for i=1:length(feedbacklevel)
k=feedbacklevel(i);
[stable, instable]=stsol(k);
ini_2=stable{2};
delta_s=ini_2(1);
initial=PARAMETER(delta_s);
sol=dde23(@LK,tao,initial,[1 tlength]);
temp=deval(sol,1:tlength);
n_t=temp(3,:); % Carrier density
kphai_t=temp(2,:);% Electric field angular frequency
A_t=temp(1,:);% Amplitude of the electric field
% % Initial the variable for phase. Add zero condition for time before
% % starting.
delta_t=zeros(1,length(kphai_t));
kphai_t1=[zeros(1,tao) kphai_t];
%
% % Calculate the phase difference for each evaluated values for future
% % use
for j=tao+1:(length(kphai_t1))
delta_t(j-tao)=kphai_t1(j)-kphai_t1(j-tao)+b;
end
%
% %length of half the data collected, hope to kill all transient data
% ti=linspace(t0,tf,tlength);
% ti=ti(clength:tlength);
n_t=n_t(clength:tlength);
delta_t=delta_t(clength:tlength);
A_t=A_t(clength:tlength);
%
% %basic graph of changing Electric field amplitude with time
% subplot(221)
% plot(A_t)
% xlabel('time')
% ylabel('Electric field amplitude')
% title(['Time seriers for feedbacklevel value of ' num2str(k)])
%
% %phase space plot for the SL with feedback i.e. to look at its carrier
% %density .vs its phase delay
%
% subplot(222)
% plot(delta_t,n_t,'.','markersize',2); %plots linearly changing phase delay vs carrier density (with reference to time)
% maxh=max(n_t); minh=min(n_t);
% maxw=max(delta_t); minw=min(delta_t);
% xlabel('Delta');ylabel('Carrier Density');
% title(['Phase space for feedback value of ' num2str(k)])
%
% %Maps phase space plot onto the surface of a cylinder and dumps it onto a new figure
% % figure(2)
% % phi=delta_t+pi; %delta_t around base of cylinder 0 to 2pi n_t; % z
% % r=ones(1,length(n_t)); %r always equals one for cylinder
% % [Xvalue,Yvalue,Zvalue] = pol2cart(phi,r,n_t);
% % plot3(Xvalue,Yvalue,Zvalue,'.','markersize',2);
% % title(['Phase Space plot mapped onto cylinder for feedbacklevel value of ' num2str(k)])
% % zlabel('Carrier Density')
% % figure(1)
%
% %Poincare plot
% subplot(224)
% I=find(abs(A_t-mean(A_t))<=1e-5);
% for i=1:length(I)
%
% if A_t(I(i))-A_t(I(i)-10)<0
% I_new(:,i)=I(i);
% end
% end
% index=find(I_new~=0);
% I_new=I_new(index);
% n_tp=n_t(I_new);
% delta_tp=delta_t(I_new);
% plot(delta_tp,n_tp,'x')
% axis([minw maxw minh maxh])
% xlabel('Delta')
% ylabel('Carrier Density')
% title(['Poincare plot for feedbacklevel value of ' num2str(k)])
%
%
% %fft power spectra graph
% subplot(2,2,3)
% ti=linspace(t0,tf,tlength);
% FFTA_t=fft(A_t);
% PowerFFTA_t=FFTA_t.*conj(FFTA_t)./(tlength-clength);
% PowerFFTA_t(1,1)=0; %sets offset term to 0
% samplerate=datalength./(ti(datalength)-ti(1)); %assumes sample rate over time is linear
% freqaxis=samplerate*(0:(datalength/2-1))/datalength; %dividing n by 2 so that the reflection is not show in the power spec
% [m,n] = size(freqaxis);
% plot(freqaxis,PowerFFTA_t(1:n));
% xlim([0 0.3]) %sets axis limits to an appropriate range to view data
% xlabel('Frequency Hz')
% ylabel('Log(Power spectral density)')
% title(['Power spectrum for feedbacklevel value of ' num2str(k)])
%
% % This part will give the 3-Dimensional View of the chaotic behaviour for the current feedback level.
%
%figure
plot3(A_t,delta_t,n_t)
xlabel('Amplitude of Electric Field')
ylabel('Phase Difference Delta')
zlabel('Carrier Density')
A_t1=[zeros(1,tao) A_t];
A_t1=A_t1(1:length(A_t1)-tao);
end
