% fdrz.m
clear
warning off
makec2

%%%%%%%%%%%%%%%%%
% USER CHOICES
%%%%%%%%%%%%%%%%%
Nr   = 20;	 % number of bins
Nz   = 2*Nr;
mua  = 0;  % [cm^-1]
musp = 10;   % [cm^-1]
b    = 0.1;  % stability factor
%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%

mutp = mua + musp; % [cm^-1]
D    = 1/(mua + musp)/3;
c    = 3e10; % [cm/s]
a    = c*D   % cm^2/s

dr   = a/b/c;% r bin size
dz   = dr;   % z bin size

%%%% Set up arrays
C    = zeros(Nz,Nr); % array for photon concentration
dC   = zeros(Nz,Nr); % array for incremental concentration per timestep
W    = zeros(Nz-2, 2*(Nr-2));

rmax = (Nr-2)*dr;
zmax = (Nz-2)*dz;
z((1:Nz)',1) = dz*((1:Nz)' - 1.5); % rows    at z
r((1:Nr)',1) = dr*((1:Nr)' - 1.5); % columns at r
rw   = (-rmax:dr:rmax);
zw   = (0:dz:zmax)';

dt   = b*dr^2/a; % time step [s]
disp(sprintf('dt = %5.3f ps', dt*1e12))

%%% draw grid
figure(1);clf
set(figure(1), 'position', [200 30 800 700])
plot([0 0],[0 zmax],'c-')
hold on
for ir=2:Nr-1
	for iz=1:Nz
		plot( [r(ir)+dr/2, r(ir)+dr/2], [-dz zmax+dz], 'c-')
		plot(-[r(ir)+dr/2, r(ir)+dr/2], [-dz zmax+dz], 'c-')
	end
end
for iz=2:Nz-1
	for ir=1:Nr
		plot([-rmax-dr rmax+dr], [z(iz)+dz/2, z(iz)+dz/2], 'c-')
		plot(r(ir), z(iz),'ro')
	end
end
plot([-rmax rmax rmax -rmax -rmax], [0 0 zmax zmax 0],'k-') % domain of computation
plot([0 0],[0 zmax],'k-')
axis([-rmax rmax -dz zmax]*1.1)
axis ij
set(gca,'fontsize',24)
xlabel('r [cm]')
ylabel('z [cm]')

% boundary conditions
for ir=2:Nr-1
	plot(r(ir),z(1), 'ro','markerfacecolor','r')
	plot(r(ir),z(Nz), 'ro','markerfacecolor','r')
end
for iz=2:Nz-1
	plot(r(1),z(iz), 'ro','markerfacecolor','r')
	plot(r(Nr),z(iz), 'ro','markerfacecolor','r')
end


%%%%%%%%%%% INITIALIZE %%%%%%%%%%%%%%%%%

C(:,:) = 0;
Co     = 1;    % [J/cm^3] photon energy concentration of laser pulse
Eo     = c*Co; % [W/cm^2] irradiance of laser
C(2,2) = musp*dt*Eo;  % concentration of light scattered in first time step
               % laser pulse duration matches time step
Cp     = Co*exp(-mutp*c*dt); % photon energy concentration remaining in laser pulse
disp(sprintf('Eo = %5.3e [W/cm2]', Eo))

time = 0;
drawC
plot([0 0],[-zmax*.7 0],'r-','linewidth',3)
text(rmax/10,-zmax/2,'Launch pulse of light','fontsize',24)
pause(1)

i = 0;
%%%%%%%%%%%%%%%%% CYCLE %%%%%%%%%%%%%%%%%
while (time <= 10e-9)
	i = i + 1;	
	time = time + dt;

% C(1,1), C(1,N), C(N,1), and C(N,N) not used

% boundary conditions
for ir=2:Nr-1
	C(1,ir) = C(2,ir);    % INSULATING top surface
	C(Nz,ir) = C(Nz-1, ir); % INSULATING bottom surface
end
for iz=2:Nz-1
	C(iz,1) = C(iz,2);    % INSULATING central z-axis
	C(iz,Nr) = C(iz, Nr-1); % INSULATING outer radial surface
end

% within volume
for ir=2:Nr-1
	for iz=2:Nz-1
		
		F1 = (C(iz,ir+1) - C(iz,ir-1))/r(ir)/2/dr;
		
		F2 = (C(iz,ir+1) - 2*C(iz,ir) + C(iz,ir-1))/dr^2;
		
		F3 = (C(iz+1,ir) - 2*C(iz,ir) + C(iz-1,ir))/dz^2;
		
		dC(iz,ir) = dt*a*(F1 + F2 + F3);
		
	end
end

% Update C
C = C + dC;
drawC

% Draw photon pulse
zs = time*c;
plot(0,zs,'ro','markerfacecolor','y', 'markersize',12,'linewidth',3)

% Scatter light from pulse
if zs<=zmax
	ii = round(zs/dz) + 2;
	C(ii,2) = C(ii,2) + musp*c*dt*Cp;
	Cp = Cp*exp(-mutp*c*dt);
end


drawnow
end % while
%%%%%%%%%%%%% end CYCLE %%%%%%%%%%%%%%%%%
