function [countRate,countRatePhotonEnergy,countRateRawData,doseRate,doseRatePhotonEnergy,doseRateRawData]=runModel(sourcePosition,sourceActivity,photonEnergy,detector,numberOfElements)
%% [countRate,countRateFromDifferentTypesOfSources,countRateRawData]=runModel(sourcePosition,sourceActivity,detector,numOfElements)
% Source position expects position of point sources in a form of a vector
% [x1,x2,...xn;y1,y2,...yn;z1,z2,...zn], where n is total number of sources
% Source activity expects matrix of activity. Each column represents
% activity of one point source. Each row represents activity of different
% nuclides.
% photonEnergy expects vector of energy of measured photons in eV
% detector expects struct which was created via function translate2struct
% numberOfElements expects single number which define total number of
% 'integrated' elements. Total number of elements will be
% numberOfElements^2
%
% countRate returns count rate of each detector.
% countRateFromDifferentTypesOfSources represents matrix where i-th page
% represents number of counts generated by i-th activity.
% countRateRawData is in form of 3D matrix. j-th line of k-th column of 
% m-th page represents contribution of counts from j-th source to k-th
% detector from m-th activity 
% doseRate returns dose rate of each detector.
% doseRateFromDifferentTypesOfSources represents matrix where i-th page
% represents energy absorbtion per second generated by i-th activity.
% doseRateRawData is in form of 3D matrix. j-th line of k-th column of 
% m-th page represents contribution of dose rate from j-th source to k-th
% detector from m-th activity 

%% testing some of the input data
if size(sourcePosition,2)~=size(sourceActivity,2)
    error('Different number of sourceActiviy and sourcePosition.')
end
if length(photonEnergy)~=length(detector.attenuation)
    error('Different number of photons energy and detector attenuations.')
end
    
%% preparation of arrays
density=detector.density;

countFunction=@(area,cubeL,r0,cubeAtt)area.*(1-exp(-cubeAtt*cubeL))./(4*pi*r0^2);
doseRateFunction=@(density,cubeL,r0,cubeAtt) 3*(1-exp(-cubeAtt*cubeL))./(4*pi*((r0+cubeL)^3-r0^3)*density);
numActivity=size(sourceActivity,1);
numDetector=size(detector.position,2);
numSource=size(sourcePosition,2);
numberOfCalcLengths=numberOfElements-1;
countRateRawData=zeros(numSource,numDetector,numActivity);
countRate=zeros(numDetector,1);

doseRateRawData=zeros(numSource,numDetector,numActivity);
doseRate=zeros(numDetector,1);

lengthInDetector=zeros(numberOfCalcLengths);

photonEnergy=photonEnergy*1.6E-19; % conversion from eV to J
for detectorIter=1:numDetector
    for sourceIter=1:numSource
        lengthInDetector=lengthInDetector*0;
        %% recalculation of r0 and d
        rs=sourcePosition(:,sourceIter);
        detectorSize=detector.size(:,detectorIter);
        cube=[pi^0.5*detectorSize(1)/2,0,0;0,pi^0.5*detectorSize(1)/2,0;0,0,detectorSize(2)];
        cubePosition=detector.position(:,detectorIter)-[cube(1,1);cube(1,1);0]/2;
        [newr0,newd]=recalculateR0andD(rs,cubePosition,cube);
        %% creating points where the length will be calculated
        dd=newd/(numberOfCalcLengths);
        nasobky=0:numberOfCalcLengths;
        
        x=newr0(1)+dd(1,1)*nasobky+dd(1,1)/2;
        y=newr0(2)+dd(2,1)*nasobky+dd(2,1)/2;
        z1=newr0(3)+dd(3,1)*nasobky+dd(3,1)/2;
        r1=[x(1:end-1);y(1:end-1);z1(1:end-1)];
        y=newr0(1)+dd(1,2)*nasobky+dd(1,2)/2;
        z=newr0(2)+dd(2,2)*nasobky+dd(2,2)/2;
        z2=newr0(3)+dd(3,2)*nasobky+dd(3,2)/2;
        r2=[y(1:end-1);z(1:end-1);z2(1:end-1)];        
        %% calculation of length
        for n=1:numberOfCalcLengths
            for m=1:numberOfCalcLengths
                lengthInDetector(n,m)=calculateLength(newr0,newd,rs,r1(:,n)+r2(:,m)-rs-newr0);
            end
        end
        %% Prepare for counting areas
        r1area=zeros(3,numberOfCalcLengths+1);
        r2area=zeros(3,numberOfCalcLengths+1);
        
        %% Create vectors which are between each vectors r1,r2 respectively
        for n=1:3
            % r1area(n,:)=rs(n)+(-r1(n,:)+dd(n,1)/2);
            r1area(n,1:end-1)=r1(n,:)-dd(n,1)/2;
            r2area(n,1:end-1)=r2(n,:)-dd(n,2)/2;
        end
        r1area(:,end)=r1(:,end)+dd(:,1)/2;
        r2area(:,end)=r2(:,end)+dd(:,2)/2;        
        %% counting probability of capture 
        p1=zeros(size(lengthInDetector,1),size(lengthInDetector,2),numActivity);
        p2=zeros(size(lengthInDetector,1),size(lengthInDetector,2),numActivity);
        x=r1area(:,2:end)-r1area(:,1:end-1);
        y=r2area(:,2:end)-r2area(:,1:end-1);
        for activityIter=1:numActivity
            for n=1:numberOfCalcLengths
                for m=1:numberOfCalcLengths
                    R1=r1(:,n)+r2(:,m)-rs-newr0;
                    p1(n,m,activityIter)=countFunction(dot(R1,cross(x(:,n),y(:,m)))/norm(R1),lengthInDetector(n,m),norm(R1),detector.attenuation(activityIter));
                    p2(n,m,activityIter)=doseRateFunction(density,lengthInDetector(n,m),norm(R1),detector.attenuation(activityIter));
                end
            end
        end
        p1(isnan(p1))=0;
        p2(isnan(p1))=0;
        countRateRawData(sourceIter,detectorIter,:)=sum(sum(abs(p1)));
        doseRateRawData(sourceIter,detectorIter,:)=sum(sum(abs(p2)));
        for n=1:numActivity
            countRateRawData(sourceIter,detectorIter,n)=countRateRawData(sourceIter,detectorIter,n).*sourceActivity(n,sourceIter);
            doseRateRawData(sourceIter,detectorIter,n)=doseRateRawData(sourceIter,detectorIter,n).*sourceActivity(n,sourceIter).*photonEnergy(n);
        end
    end
    countRate(detectorIter)=sum(sum(countRateRawData(:,detectorIter,:)));
    doseRate(detectorIter)=sum(sum(doseRateRawData(:,detectorIter,:)));
end
countRatePhotonEnergy=sum(countRateRawData,1);
doseRatePhotonEnergy=sum(doseRateRawData,1);
end

function l=calculateLength(r0,d,rs,v)
% r0 is point, which is the vertex of vectors d1,d2,d3
% d is matrix of vectors edges of detector d=[d1 d2 d3]
% rs is position of source
% v is vector defining direction of gamma ray

% n_i are vectors of planes.
% when asuming orthogonality of vectors d1,d2,d3 then
% n1=d1/abs(d1)
% n2=d2/abs(d2)
% n3=d3/abs(d3)
epsilon=1e-8;
d1=d(:,1);
d2=d(:,2);
d3=d(:,3);
n1=cross(d2,d3);
n2=cross(d3,d1);
n3=cross(d1,d2);

n1=n1/norm(n1);
n2=n2/norm(n2);
n3=n3/norm(n3);
n=[n1,n2,n3];

%ha is vector of offsets of planes with normal vectors n1,n2,n3 which are
%connected in vertex r0
%ta: A=rs+ta*v
ha=dot(n,[r0,r0,r0]);
ta=(ha-dot(n,[rs,rs,rs]))./dot([v,v,v],n);
A=[rs,rs,rs]+[ta;ta;ta].*[v,v,v];
%hb is vector of offsets of planes with normal vectors n1,n2,n3 which are
%connected in the opposite vertex of r0
%tb: A=rs+tb*v
hb=dot(-n,[r0,r0,r0]+d);
tb=(hb-dot(-n,[rs,rs,rs]))./dot([v,v,v],-n);
% A represents which points are taken into consideration for entering and
% exiting the cube
A(:,4:6)=[rs,rs,rs]+[tb;tb;tb].*[v,v,v];

for k=1:6
    % erases points which are not on the cube
    A(:,k)=inOut(A(:,k),d(:,1),d(:,2),d(:,3),r0);
end
relevantPoints=~isnan(sum(A));
points=A(:,relevantPoints);
if size(points,2)==2
    % particle enters the cube through one side and exit through another
    l=norm((points(:,1)-points(:,2)));
elseif size(points,2)<2
    % particle comes through corner
    l=0;
else
    % particle comes in through corner of cube
    l=0;
    for k=1:size(points,2)
        if sum(abs(points(:,1)-points(:,k)))> epsilon
            l=norm((points(:,1)-points(:,k)));
            break;
        end
    end
end
end
function res=inOut(A,d1,d2,d3,r0)
% testign if point A lays inside or outside of the cube
reserve=1e-8;
ra=r0-A;
denominator1=d1(3)*d2(2)*d3(1)-d1(2)*d2(3)*d3(1)-d1(3)*d2(1)*d3(2)+d1(1)*d2(3)*d3(2)+d1(2)*d2(1)*d3(3)-d1(1)*d2(2)*d3(3);
K1=(-d2(3)*d3(2)*ra(1)+d2(2)*d3(3)*ra(1)+d2(3)*d3(1)*ra(2)-d2(1)*d3(3)*ra(2)-d2(2)*d3(1)*ra(3)+d2(1)*d3(2)*ra(3))/denominator1;
K2=(d1(3)*d3(2)*ra(1)-d1(2)*d3(3)*ra(1)-d1(3)*d3(1)*ra(2)+d1(1)*d3(3)*ra(2)+d1(2)*d3(1)*ra(3)-d1(1)*d3(2)*ra(3))/denominator1;
K3=(-d1(3)*d2(2)*ra(1)+d1(2)*d2(3)*ra(1)+d1(3)*d2(1)*ra(2)-d1(1)*d2(3)*ra(2)-d1(2)*d2(1)*ra(3)+d1(1)*d2(2)*ra(3))/denominator1;

if (K1<0-reserve||K2<0-reserve||K3<0-reserve||K1>1+reserve||K2>1+reserve||K3>1+reserve)
    res=nan*A;
    
else
    res=A;
end

end
function [newR0,newD]=recalculateR0andD(rs,r0,d)
% finds nearest corner of the cube and recalculates r0 and d accordingly
distRSandClosestCorner(1)=norm(-rs+r0);
distRSandClosestCorner(2)=norm(-rs+r0+d(:,1));
distRSandClosestCorner(3)=norm(-rs+r0+d(:,2));
distRSandClosestCorner(4)=norm(-rs+r0+d(:,3));
distRSandClosestCorner(5)=norm(-rs+r0+d(:,1)+d(:,2));
distRSandClosestCorner(6)=norm(-rs+r0+d(:,1)+d(:,3));
distRSandClosestCorner(7)=norm(-rs+r0+d(:,2)+d(:,3));
distRSandClosestCorner(8)=norm(-rs+r0+d(:,1)+d(:,2)+d(:,3));

closestPoint=find(distRSandClosestCorner==min(distRSandClosestCorner));

switch closestPoint(1)
    case 1
        newR0=r0;
        newD=d;
    case 2
        newR0=r0+d(:,1);
        newD=[-d(:,1),d(:,2),d(:,3)];
    case 3
        newR0=r0+d(:,2);
        newD=[d(:,1),-d(:,2),d(:,3)];
    case 4
        newR0=r0+d(:,3);
        newD=[d(:,1),d(:,2),-d(:,3)];
    case 5
        newR0=r0+d(:,1)+d(:,2);
        newD=[-d(:,1),-d(:,2),d(:,3)];
    case 6
        newR0=r0+d(:,1)+d(:,3);
        newD=[-d(:,1),d(:,2),-d(:,3)];
    case 7
        newR0=r0+d(:,2)+d(:,3);
        newD=[d(:,1),-d(:,2),-d(:,3)];
    case 8
        newR0=r0+d(:,1)+d(:,2)+d(:,3);
        newD=[-d(:,1),-d(:,2),-d(:,3)];
        
end
end



