%YL, 08/06/2017

filename = 'example';
np = 1e6;
layers_matrix = [1, 0.1 100 0.9 50];
n_above = 1;
n_below = 1;
dz = 0.005;
dr = 0.01;
nz = 200;
nr = 200;
na = 1;
h = create_MCML_input_file(filename,np,layers_matrix,n_above,n_below,...
   dz,dr,nz,nr,na);
disp([h.input_filename,' file created']);
dos(['Mcml.exe ' h.input_filename '&']);
while isempty(dir(h.output_filename)) %terminate upon simulation done
   pause(1);
end
disp('Simulation done!');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Parameter Space

%read output
fid = fopen(h.output_filename);
for i=1:14                      % Tar bort onödiga rader
   tline = fgetl(fid);
end
step_size = fscanf(fid,'%g'); %[dz, dr]
tline = fgetl(fid);
step_num = fscanf(fid,'%g'); %[nz,nr,na]
%RAT
rad='abc';                      
crit='RAT'; %reflectance, absorption and transmittance
while rad~=crit
   tline = fgetl(fid);
   if length(tline)>=3
       rad=tline(1:3);
   end
end
spec_refl = fscanf(fid,'%g');
tline = fgetl(fid);              
diff_refl = fscanf(fid,'%g');    
tline = fgetl(fid);              
abs_frac = fscanf(fid,'%g');    
tline = fgetl(fid);              
trans_frac = fscanf(fid,'%g');
%A_l
for i=1:3
   tline = fgetl(fid);          
end
abs_layer = fscanf(fid,'%g');  
%A_z
tline = fgetl(fid);  
abs_z = fscanf(fid,'%g');
%Rd_r
tline = fgetl(fid);
refl_r = fscanf(fid,'%g');
%Rd_a
tline = fgetl(fid);
refl_a = fscanf(fid,'%g');
%trans_r
tline = fgetl(fid);
trans_r = fscanf(fid,'%g');
%trans_a
tline = fgetl(fid);
trans_a = fscanf(fid,'%g');
%A_rz
for i=1:6
   tline = fgetl(fid);          
end
abs_rz = fscanf(fid,'%g');
abs_rz = reshape(abs_rz,step_num(1),step_num(2));
%refl_ra
for i=1:6
   tline = fgetl(fid);          
end
refl_ra = fscanf(fid,'%g');
refl_ra = reshape(refl_ra,step_num(2),step_num(3));
%trans_ra
for i=1:6
   tline = fgetl(fid);          
end
trans_ra = fscanf(fid,'%g');
trans = reshape(refl_ra,step_num(2),step_num(3));
%compute fluence rate: A(r,z)/mua
for i=1:step_num(1)
   for j=1:step_num(2)
       flod(i,j) = abs_rz(i,j)/layers_matrix(2);
   end
end
f_rz=flod;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Comparsion with analytical solutions under structured illumination

mua = 0.2;
mus = 4;
Fx = linspace(0.2,3.6,18);
Rd = Rd_fx(mua,mus,Fx);

layers_matrix = [1.35, 0.2, 10*mus, 0.9, 5];
dr = 0.01;
dz = 0.01;
nr = 500;
nz = 500;
s = MCML('SI_test',5e6,layers_matrix,1,1,dz,dr,nz,nr,1);

refl = s.refl_r;

refl_MC = hankelTrans(refl,Fx,dr);

figure
plot(Fx, Rd, '-rd','linewidth',1.5); hold on;
plot(Fx, refl_MC, '-cs','linewidth',1.5);
xlabel('Spatial frequency (mm^{-1})');
ylabel('Reflectance');
legend('Analytical','Monte Carlo');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Comparsion with Finite Element