## Copyright (C) 2015 Loriane Weber
## 
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## PhantTibo_Analytical

## Author: Loriane Weber <lweber@rnice7-0102>
## Created: 2015-12-10

%%%%%% INPUT PARAMETERS
%% ps : is the physical pixel size, in um
%% RedAttFact : reduced attenuation factor
%% angles : angles of projections
%% save : 1 or 0 (1 if you want to save in edf, 0 sinon

function [MuP_tot, Phi_tot, B_tot ] = PhantTibo_Analytical (angles, ps, save)

disp('params');

RedAttFact=1; 

mu_Al=10.7765 * 100/RedAttFact; % *100 to be in m-1 
mu_Mg= 5.56697 * 100/RedAttFact; 
mu_PET=0.89889 * 100/RedAttFact; 

delta_Al=15.03*10^-7/RedAttFact;
delta_Mg=9.917*10^-7/RedAttFact;
delta_PET=8.265*10^-7/RedAttFact;

beta_Al=5.59*10^-9/RedAttFact;
beta_Mg=2.89*10^-9/RedAttFact;
beta_PET=0.46*10^-9/RedAttFact;

R_Al=73/2;
R_Mg=37/2;
R_PET=58/2;

cen_Al=[607-600; 585-600]; % 600 to put the center at the center of the image (1200*1200)
cen_Mg=[679-600; 640-600];
cen_PET=[730-600; 630-600];

du= 1; % echantillonage de la projection

angles_rad=(angles*pi/180) + pi/2; % in radians % +pi/2 to be coherent with radon transform

dim=1701; % taille de la proj, wrt radon function
u=[-(dim-1)/2:du:+(dim-1)/2]; 
disp('proj')


for i=1:length(angles)

	u0_Al= cos(angles_rad(i))*cen_Al(1) + sin(angles_rad(i))*cen_Al(2);
	u0_Mg= cos(angles_rad(i))*cen_Mg(1) + sin(angles_rad(i))*cen_Mg(2);
	u0_PET= cos(angles_rad(i))*cen_PET(1) + sin(angles_rad(i))*cen_PET(2);

## Proj attenuation
	MuP_Al=Analytical_1D(mu_Al, u, R_Al, u0_Al);
	MuP_Mg=Analytical_1D(mu_Mg, u, R_Mg, u0_Mg);
	MuP_PET=Analytical_1D(mu_PET, u, R_PET, u0_PET);

	MuP_tot(:, i) = ( MuP_Al+ MuP_Mg + MuP_PET ) * ps;

## Proj Phi
	Phi_Al=Analytical_1D(delta_Al, u, R_Al, u0_Al);
	Phi_Mg=Analytical_1D(delta_Mg, u, R_Mg, u0_Mg);
	Phi_PET=Analytical_1D(delta_PET, u, R_PET, u0_PET);

	Phi_tot(:, i) = ( Phi_Al+ Phi_Mg + Phi_PET ) * ps;

## Proj Beta
	B_Al=Analytical_1D(beta_Al, u, R_Al, u0_Al);
	B_Mg=Analytical_1D(beta_Mg, u, R_Mg, u0_Mg);
	B_PET=Analytical_1D(beta_PET, u, R_PET, u0_PET);

	B_tot(:, i) = ( B_Al+ B_Mg + B_PET ) *ps;

end

if save
	edfwrite('Sino_Analytical_Mu.edf', MuP_tot, 'float32')
	edfwrite('Sino_Analytical_Phi.edf', Phi_tot, 'float32')
	edfwrite('Sino_Analytical_B.edf', B_tot, 'float32')
end


endfunction
