clear all
% consider all rotational levels between j=0 and j=jmax;
jmax=15;
j=[0:jmax];
% g is the rotational degeneracy
g=2*j+1;
% pick a rotational constant (in cm-1)
b=10;
% ej is the energy of each rotational level
ej=b*j.*(j+1);
% kmax is the upper bound to any individual state-to-state rate
kmax=5.e-11; 
alpha=2;
% T is the temperature in Kelvin
T=1000;
% rho is the number density of collision partners
rho=3.2e22;
% use an exponential gap model to estimate the lower triangle of the
% matrix of rotationally inelastic rate constants
k=zeros(jmax+1,jmax+1);
for ii=2:jmax+1
   for jj=1:ii-1
      k(ii,jj)=kmax*exp(-alpha*(ej(ii)-ej(jj))*1.4388/T);
   end
end
% now use detailed balance to determine the upper triangle of the matrix
% delta-E factor
deltaE=repmat(ej,jmax+1,1);
deltaE=deltaE-deltaE';
edeltaE=exp(deltaE*1.4388/T);

% degeneracy factor
degen=repmat(g,jmax+1,1);
degen=degen'./degen;

k=k+tril(k,1)'.*degen.*edeltaE;

% determine diagonal elements to conserve total number
vec=sum(k);
k=k-diag(vec);

%pboltz is the vector of boltzmann populations
pboltz=g.*exp(-ej*1.4388/T);
pboltz=pboltz/sum(pboltz);
% pi is the diagonal matrix of boltzmann populations
pi=diag(pboltz);
sqpi=diag(sqrt(pboltz));
% sqpii is the inverse of the sqpi
sqpii=diag(1./sqrt(pboltz));
% transform rate matrix into symmetric form

ksym=sqpii*k*sqpi;

% determine eigenvectors and eigenvalues of this matrix
[D Lambda]=eig(ksym);

% consider initial population in one level 
j0 = 6;
if j0>jmax
   disp(['j0 = ' int2str(j0) ', which is larger than jmax = ' int2str(jmax) '; abort'])
   return
end

n0=zeros(jmax+1,1);n0(j0+1)=1;


% determine initial ntilde vector
ntilde0= sqpii*n0;
% determine initial N-script vector

scriptn0=D'*ntilde0;

ev=rho*diag(Lambda);
% dt is the time step
dt=1.e-14;
for i=1:201
   t(i)=dt*(i-1);
   % scriptn is the population in the n-tilde basis at time t
   scriptn=scriptn0.*exp(ev*t(i));
   % n is the population in the n basis at time t
   n=sqpi*D*scriptn;
   erotav(i)=sum(ej*n)/sum(n);
   figure(1)
   bar(j,n);xlabel('j');ylabel('p(j)');title(['T = ', int2str(T) ' K; t = ' num2str(t(i)) ' s'])
   axis([-0.5 jmax+0.5 0 max(n)*1.2])
   set(gca,'fontsize',14)
   %pause(1)
end
eboltz=pboltz*ej';
figure(3);plot(t,erotav,[0 max(t)],eboltz*[1 1]) 
set(gca,'fontsize',14);
xlabel('time / s');ylabel('< Erot >');title('Average rotational energy')

disp('hit any key to look at overall infill of an initially bleached level')

% n is the initial n-vector

n0=diag(pi);
n0(j0+1)=0;
n0=n0/sum(n0);
nmax=max(n0)*1.2;

% determine initial ntilde vector
ntilde0= sqpii*n0;
% determine initial N-script vector

scriptn0=D'*ntilde0;

ev=rho*diag(Lambda);
% dt is the time step
dt=1.e-14;
for i=1:200
   t(i)=dt*(i-1);
   % scriptn is the population in the n-tilde basis at time t
   scriptn=scriptn0.*exp(ev*t(i));
   % n is the population in the n basis at time t
   n=sqpi*D*scriptn;
   erotav(i)=sum(ej*n)/sum(n);
   figure(1)
   
   bar(j,n);xlabel('j');ylabel('p(j)');title(['T = ', int2str(T) ' K; t = ' num2str(t(i)) ' s'])
   axis([-0.5 jmax+0.5 0 max(n)*1.2])
   set(gca,'fontsize',14)
   pause(1)
end
figure(3);plot(t,erotav,[0 max(t)],eboltz*[1 1]) 
set(gca,'fontsize',14)
xlabel('time / s');ylabel('< Erot >');title('Average rotational energy')
