% fdrz.m
clear
warning off

%%%%%%%%%%%%%%%%%
% 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   % diffusivity [cm^2/s]
alpha = a;   % diffusivity [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

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


% 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) = Co;   % concentration of light scattered in first time step
               % laser pulse duration matches time step
Uo     = Co*pi*dr^2*dz;

%Cp     = Co*exp(-mutp*c*dt); % photon energy concentration remaining in laser pulse
disp(sprintf('Uo = %5.3f [J]', Uo))

time = 0;
drawCrz
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;
drawCrz

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

% Scatter light from pulse
if 0 %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 %%%%%%%%%%%%%%%%%
