function [dosePrefix,nSlices] = getDose(targetFolder,path,fPrefix, bldTh, expT,res,machineparamobj, resinparamobj)
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here



d = 0:1e-6:1000e-6;
beers = 1*machineparamobj.flatF*exp(-d/resinparamobj.critD);  % Gives dose / height relationship
ind = find(beers<(0.03*beers(1)),1);  % Finds first array element where dose energy is 3% of surface
maxLayers = d(ind)/bldTh;  % Max number of layers to calculate

top = cd(path);
a = dir([path '/*.bmp']);
nLayers = size(a,1)-1;

mkdir(targetFolder);
for cLayer = 1:nLayers
    % Load the images of the layers
    pxFile = strcat(fPrefix,num2str(cLayer,'%03.0f'),'.bmp');
    % Expecting a filename of format filePrefix_01.bmp
    pxTemp = imbinarize(rgb2gray(imread(pxFile)));
    % Read current pxMap layer from file
    pxTemp = pxTemp == 0;  % Convert to logical
    space = 4;  % Add 4 pixels spacing in each side, should be 110 um
    pxMap = zeros(size(pxTemp,1)+2*space,size(pxTemp,2)+2*space);
    pxMap(space+1:end-space, space+1:end-space)=pxTemp;
    
    
    [xSize,ySize] = size(pxMap);
    xRes = xSize*machineparamobj.rhoX/res;
    yRes = ySize*machineparamobj.rhoY/res; % Size of intensity map
    [X,Y] = meshgrid(0:res:xRes*res,0:res:yRes*res);
    
    iMap = getLayerIntensity(pxMap, xSize, ySize, X, Y, machineparamobj.pxIntensity, machineparamobj.gaus2W, machineparamobj.rhoX, machineparamobj.rhoY);
    
    clear curedZ;
    
    if (cLayer - maxLayers) < 0
        curedZ = 0:res:cLayer*bldTh;
    else
        curedZ = (cLayer-maxLayers)*bldTh:res:cLayer*bldTh;
    end
    layerID = int16(curedZ(1)/res);
    cd(targetFolder);
    
    for z = curedZ
        thisLayer = getDoseLayer(iMap,z,expT, cLayer,bldTh,resinparamobj.critD);
        thisLayer((1-exp(-resinparamobj.K*thisLayer))<resinparamobj.phi_c) = 0;
        file = strcat('dose_',num2str(layerID,'%02.0f'),'.csv');
        if isfile(file)
            prevLayer = readmatrix(file);
            newLayer = prevLayer+thisLayer;
            writematrix(newLayer,file);
            clear thisLayer;
        else
            writematrix(thisLayer,file);
        end
        layerID=layerID+1;
    end
    cd(path)
    clear curedZ;
    
end


a = dir('*.csv');
dosePrefix = "dose_";
nSlices = size(a,1);

end

function iMap = getLayerIntensity(pxMap,xSize,ySize,X,Y,pxIntensity,gaus2W,rhoX,rhoY)
    pxProfile = zeros(size(X));
    iMap = zeros(size(X));
    for i=1:xSize
        for j=1:ySize
            % Sweep over every pixel in the array
            if pxMap(i,j)~=0
                xPxCentre = rhoX*(i-0.5);
                yPxCentre = rhoY*(j-0.5);
                distSq = (X-xPxCentre).^2+(Y-yPxCentre).^2;
                flag = abs((Y-yPxCentre)./(X-xPxCentre))<=1;
                alpha = Y-yPxCentre;  % Default assumes voxel is closer to centre in X
                alpha(flag) = X(flag)-xPxCentre;  % If pixel is to right or left use this value
                omega = gaus2W*sqrt(distSq)./alpha;
                pxProfile = pxIntensity*exp(-(2./(omega.^2)).*distSq);
                pxProfile(alpha==0)=pxIntensity;
                iMap = iMap + pxProfile;
                % Add current pixel to total for this layer
            end
        end
    end
end

function dose = getDoseLayer(iMap,curedZ,expT,cLayer,bldTh,critD)
    % For each pixel in iMap get the dose energy at the current z
    dose = expT*iMap*exp((curedZ-cLayer*bldTh)/critD);
end