% Non-line-of-sight multiple-scatter propagation model
% Author: Hongwei Yin(National University of Defense Technology,China.
% E-mail:yinhongwei@nudt.edu.cn)
% Please cite: H. Yin, S. Chang, H. Jia, J. Yang, and J. Yang,¡°Non-line-of-sight
% multiscatter propagation model¡±, J. Opt. Soc. Am. A 26, 2466~2469 (2009).

% Hats off to Craig F. Bohren.

% The name of every variable or functions have an added "_" to aviod
% homonymy with the built-in variable or function of matlab.
clear 
%The transmitter and receiver are structure, they have five same
%parameters: theta,fai,positioninoriginal,delta and normal. The unit of the
%angles is radian, and that of the range is meter.
transmitter_=struct('theta',pi/2-45/180*pi,'fai',pi/2,'positioninoriginal',...
    [1 1 1]','delta',15*pi/180);
receiver_=struct('theta',pi/2-45/180*pi,'fai',3*pi/2,'positioninoriginal',...
    [1 500+1 1]','delta',15*pi/180);
%Calculating the normals of the transmitter and the receiver 
[transmitter_,receiver_]=angle_normal(transmitter_,receiver_);
%ke_ is the extinction coefficient, and ks_ the scattering coefficient (unit: /m)
ke_=1*10^(-3);ks_=0.5*10^(-3);
%g_ is the g-factor of Henyey-Greenstein scattering phase function
g_=0;
%photonoriginalweight_ is the original weight of each photon,
%totalphotonnumber_ is the number of the photons in simulation,
%maxscattercount_ is the max scattering count considered for each photon. 
photonoriginalweight_=10^10; totalphotonnumber_=10^4; maxscattercount_=4;
%range resolution of the scattered irradiance (unit: m)
rangereso_=10;
%result_.rangedis saves the scattered irradiance versus the
%range,result_.scattercountdis saves the scattered irradiance versus the
%scattering count
result_.rangedis(maxscattercount_,1:round(10^5/rangereso_))=0;
result_.scattercountdis(1:maxscattercount_)=0;

%Simulation start
for photonnumber_=1:totalphotonnumber_
    %birth of one single photon, 
    [photon_]=transmitter_emit(transmitter_,photonoriginalweight_);
    %scattering start
    for scattercount_=1:maxscattercount_
        %sampling the free path of the photon, updating the photon
        [photon_]=photon_path(photon_,ke_,ks_);
        %z=0 is set as the ground, z<0 means that the photon hit the ground
        if (photon_.positioninoriginal(3)<=0)
            photon_.weight=0;
        end
        [pointto_]=pointto_receiver(photon_.positioninoriginal,receiver_,g_);
        [result_,photon_]=receiver_receive(photon_,pointto_,result_,ke_,rangereso_);
        [photon_]=photon_scatter(photon_);
    end
    photonnumber_
end

time_=(1:length(result_.rangedis(1,:)))/3/10^8*rangereso_;
result_.rangedis=result_.rangedis/totalphotonnumber_/photonoriginalweight_;
result_.scattercountdis=result_.scattercountdis/totalphotonnumber_/photonoriginalweight_;

subplot(2,1,1)
plot(time_(1:200),result_.rangedis(1,1:200)*3*10^8/rangereso_)
xlabel('time (s)');ylabel('single-scatter irradiance (W/m2)')
subplot(2,1,2)
bar(1:maxscattercount_,result_.scattercountdis)
xlabel('scattering count');ylabel('scattered irradiance (W/m2)')

