% Fig_4_3.m 
% classical oscillator using the RK4 algorithm to see time evolution.
% The different physical condtions can be achived by adjusting the
% parameters. Up to 3rd-order non-linearity is included.
% Linear term is linear_coeff*(omegax)^2 and 3rd-order term is
% non_linear_3*(omegax)^2.

%plotting parameters + fontsizes
clear all;
clf;

FS      = 24;	%label fontsize 18
FSN     = 24;	%number fontsize 16
LW      = 2;	%linewidth 2

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', FSN)

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times New Roman')
set(0,'DefaultTextFontSize', FSN)


%%%%%%%%%%%%%%%% input parameters data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

omegax=pi*1e+06;    % unit:s-1, frequency of the oscillator. 3.2698e+06 Hz
domegax=omegax*0.10; % frequency shift (0, 0.01, 0.1, 1.0)
omegax1=(pi*1e+06)+domegax; % unit:s-1, frequency of the periodic force 
                            % applied to the oscillator. 3.2691e+06

gamma = omegax / 10;% unit:s-1, oscillator damping term (200, 500)

linear_coeff =+1; % linear simple harmonic term, +ve (1) for SHO 
                  % and -ve (-1, -1.5, -2.5) for Duffing oscillator

non_linear_3 =0.0; % coefficient of 3rd order nonlinear term. (0, 1, 10)


%Example of SHO use /200, 1.0, 0.0, 50e12
%Example frequency doubling use /200, -2.5, 1.0, 50e12
            
total_time=0.01; % unit:s, total time for which system is observed. 0.1

fraction1= 0.01; % time increment of integration = fraction1 / omegax
%it is a small fraction of the oscillatory period 
%set to smaller value for smaller time steps.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% enter initial values %%%%%%%%%%%%%%%
y1= 0; % initial VELOCITY of oscillator at t = 0.0
y2= 0; % initial POSITION of oscillator at t = 0.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

deltat1=(fraction1*2*pi/omegax);% unit:s, time increment of integration, 
                                % set to 1% of oscillation period
deltat2=(fraction1*2*pi/omegax1);% unit:s, time increment of integration, 
                                % set to 1% of oscillation period

%smaller of the time increments chosen for integretion.
arr = sort([deltat1 deltat2]);
deltat = arr(1);

extent=round(total_time/deltat) % total number of time steps determined 
                                % from time increment.


ss = 1;
if gamma>0
%ss = round((50/gamma)/deltat) % remove initial transient decay 
                              % from position array
end
F0(1:extent) = 50e12; % amplitude of periodic force applied to oscillator 
                      % in absence of noise. 1e11 - 1e14
% F0(ss+round(0.005*extent)-1:extent)=-6335e12;%-6335e12
% F0(ss+round(0.005*extent))=+6335e12;%-6335e12
% F0(ss+round(0.005*extent):extent)=-6340e12;%-6335e12
% F0(ss+round(0.005*extent)+50:extent)=0.0;
%%%%%%%%%%%%% initialising storage matrices%%%%%%%%%%%%%%%%%%%%%%%%%%

Y(extent,2)=0;% array with 4 coloumns, column 1 : velocity 1,  
% column 2 : position 1, column 3 : velocity 2,  column 4 : position 2.

%T(1:extent)=0; % stores time values at every integration point.
T=ones(1,extent); %initialize time array
T=deltat*double(find(T));%initialize time values at every integration point
lcoeff=linear_coeff*(omegax^2);% redefine coeff_of_linear
nlcoeff =non_linear_3*(omegax^2);% redefine 3rd-order nonlinear term

Force(1:extent)=0; % stores force values at every integration point.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% integration loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for index=1:extent   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RK4 algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%

Force(index) = (F0(index))*cos(omegax1*T(index)) ; % update force array

temp1=y1;%velocity
temp2=y2;%position
       
        % 1st step 
        
         F = (F0(index))*cos(omegax1*T(index)) ;
         d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1 + F;%cos(omegax1*t);
         d2 = y1;
         
         ak1=d1*deltat;
         bk1=d2*deltat;
                  
         y1=temp1+0.5*ak1;
         y2=temp2+0.5*bk1;
                 
        %2nd step
         
         F = (F0(index))*cos(omegax1*(T(index)+deltat/2));
         d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1+F;%cos(omegax1*t);
         d2 = y1;
               
         ak2=d1*deltat;
         bk2=d2*deltat;
                
         y1=temp1+0.5*ak2;
         y2=temp2+0.5*bk2;
                 
         %3rd step
         
         F = (F0(index))*cos(omegax1*(T(index)+deltat/2)) ;
         d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1+F;%cos(omegax1*t);
         d2 = y1;
                
         ak3=d1*deltat;
         bk3=d2*deltat;
                 
         y1=temp1+1*ak3;
         y2=temp2+1*bk3;
                
         %4th step
         F = (F0(index))*cos(omegax1*(T(index)+deltat)) ;
         d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1+F;%cos(omegax1*t);
         d2 = y1;
                 
         ak4=d1*deltat;
         bk4=d2*deltat;
                 
         y1 = temp1+(ak1+2.0*ak2+2.0*ak3+ak4)/6.0;
         y2 = temp2+(bk1+2.0*bk2+2.0*bk3+bk4)/6.0;
                 
         % store position and velocity values.
         Y(index,1)=y1;
         Y(index,2)=y2;        
end
        
save Y.mat Y
save T.mat T
%Plotting the columns of the returned array Y versus T shows the solution 
ttl=(['lin=',num2str(linear_coeff),', nonlin=',...
    num2str(non_linear_3),', \it{f}\rm_0 =',num2str(omegax/1e6, '%7.2f')...
    ,' MHz, \it{f}\rm_F =',num2str(omegax1/1e6, '%7.2f')...
    ' MHz, \gamma =',num2str(gamma/1e6, '%7.4f'),' MHz']);
figure(1)
plot(T(ss:ss+round(0.005*extent))*1e6/(2*pi),...
    (Y(ss:ss+round(0.005*extent),2)),'k','LineWidth',LW);
xlabel('Time, \it{t}\rm (\mus)','FontSize',FS);
ylabel('Position,  \it{x}\rm (\it{t}\rm)','FontSize',FS);
axis([T(ss)*1e6/(2*pi),T(ss+round(0.005*extent))*1e6/(2*pi),...
    1.1*min(Y(:,2)),1.1*max(Y(:,2))]);
%title(ttl,'FontSize',FS);
set(gca,'LineWidth',LW)
set(gca,'TickLength',[0.02 0.025])

figure(2)
plot(Y(ss:end,2)/1000,(Y(ss:end,1)/1000),'k','LineWidth',LW);
axis square;
xlabel('Position, \it{x}\rm(\it{t}\rm)','FontSize',FS);
ylabel('Velocity, d \it{x}\rm/d\it{t}\rm','FontSize',FS);
%title(ttl,'FontSize',FS);
set(gca,'LineWidth',LW)
set(gca,'TickLength',[0.02 0.045])
axis ([-0.03 0.03 -1e5 1e5]);