	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%			                      %
	%   FILE NAME: forcedHO           %
	%			                      %
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% This program accompanies the book: "Introduction to Quantum Mechanics:
% A Time Dependent Perspective", by David Tannor (University Science Books, 2003).



% This program was written by Erez Boukobza in winter 2002.


% This program caculates and plots the phase space for the forced harmonic oscilator.
% The calculation is based on the proposed solution for differential equations
% as appears on exercise 9b on chapter 3 in the book (with a correction for t0).
% The plotting in the second part (corresponds to an interruption in the forcing)
% may be commented out in order to run the program with one forcing.
% Written by Erez Boukobza winter 2002.
% Variables and parameters for the calculation

syms x t1 t2; %symbolic variables, where x is the time variable. 
phi=pi; % The phase of the interrupting forcing, can be changed.
tmf=8*pi; % The interruption time (the time when the forcing changes),can be changed. 
m=1; % mass 
w=1; % basic frequency of the oscilator.
w1=1; % frequncy for a forcing of a form sin(w1*t), should be changed.
x0=0; % initial position
p0=0; % initial momentum
t0=0; % initial time
A=0.05; % constant for a forcing of the form A or A*t, can be changed.

% The forcing can be changed in the line below to: A, A*t, sin(w1*t), or anything else.
force1=A*x; 



% First part: Phase space calculation for a forced harmonic oscilator
z0=p0+i*m*w*x0;

% In the line below the trajectory is calculated symbolically according to the
% Complex quantity method:
% Z(t)= Pt + imwXt = e(iwt)*int(F(t')*e(-iwt'),t0,t) + e(iw(t-t0))*Z(t0)
z1=exp(i*t1*w)*int(exp(-i*x*w)*force1,t0,t1)+exp(-i*t0*w)*exp(i*t1*w)*z0;

% Substitution into the symbolic quantity and plotting. 
xx=[0:pi/100:tmf];
Sz1=subs(z1,t1,xx);
plot(Sz1);
xlabel('pt');
ylabel('mwXt');
title('Phase space for the forced harmonic oscilator');
hold on;



% Second part: Change of forcing (or phase of forcing sin(w1*t+phi))

% The plot command line can be commented out to run the program only with one forcing
% (no interruption).

% In the three lines below  Z(t) is calculated at the interruption time and after
% the forcing has been changed according to the complex quantity method.
z1t1f=exp(i*tmf)*int(exp(-i*x)*force1,t0,tmf)+exp(-i*t0*w)*exp(i*tmf*w)*z0; 
force2=sin(w1*x+phi);
z2=exp(i*t2*w)*int(exp(-i*x)*force2,tmf,t2)+exp(-i*tmf*w)*exp(i*t2*w)*z1t1f;

% Substitution into the symbolic quantity and plotting. 
xx=[0:pi/100:tmf];
yy=[tmf:pi/100:16*pi];
Sz2=subs(z2,t2,yy);

%plot(Sz2,'r'); % This line can be commented out to run one forcing only
% The red line corresponds to the trajectory after the change of phase
