clf
clear
tcount=1;

%Radiation count initialisation
msvcount=0;
eta=0.001/7; % mSv (arbitrary per particle)
RADMEASURE=0; % in mSv

%%%Coastline definition%%%
%Load coast
close all
I=imread('Snip20141109_2.png');
imshow(I)

level=graythresh(I);
bw=im2bw(I,level);
figure,imshow(bw)
bw=double(bw);


for i=1:size(bw,1)
    for j=1:3
        bw(i,j)=1;
    end
    
end


for i=1:3
    for j=1:size(bw,2)
        bw(i,j)=1;
    end
    
end

bw=flipud(bw);

for i=1:size(bw,1)
    for j=1:size(bw,2)
        if bw(i,j)==0
            bw(i,j)=1;
        else bw(i,j)=0;
        end
    end
end

%GRID
maxX=size(bw,1);
maxY=size(bw,2);
nParticles=7000;

NNPPF=500;

%Variable initialisation
x=0;
y=0;
stuck=0;
particle=0;
escape=0;
die=0;
r=0;
%Movie tracker/counter
f=1;

while (( particle <nParticles) + (escape~=0))~=0
    stuck = 0; % important
    die = 0;% important
    
    %Particle iteration updater/counter
    particle=particle+1;
    
    %Time tracker
    tcount=tcount+1;
    tcountme=tcount*(1/140);
    
    %Generate point source from one of a selection of used/disused nuclear
    %power stations
    aran=rand();
    
    %Function that takes in the value of 'aran' and assigns a pointed
    %source (weighted in line with consultation of relevant literature)
    % See function file for more information.
    sourcegen = sourcefunk(aran,x,y);
    
    %Assign x and y values accordingly
    x=sourcegen(1);
    y=sourcegen(2);
    
    %value storage
    old_x=x;
    old_y=y;
    
    stuck=0;
    
    %Eoghan Totten added a 'particle' print command - relating to iteration
    fprintf('new: bw(%d,%d)=%f %d \n',x,y,bw(x,y),particle);
    while (stuck == 0)&&(die == 0)
        old_x=x;
        old_y=y;
        r=rand();
        
        %%%OCEAN CURRENT FUNCTION IMPLEMENTATION
        out = oceancurrent(r,x,y);
        
        %Once again: Assign x and y values accordingly
        x=out(1);
        y=out(2);
        
        %%% Remedy: POINT SOURCE CANNOT SIT ON COAST
        if (bw(x,y)==1)
            x=old_x;
            y=old_y;
            %Burkhard's print contributions
            fprintf('land: bw(%d,%d)=%f %d \n',x,y,bw(x,y),particle);
            continue;
        end
        
        %Neighbour definition
        xRight = x+1;
        xLeft = x-1;
        yUp = y+1;
        yDown = y-1;
        %%%%%%%%%%%%%%%%
        
        %Weighted sticking according to region on grid
        p=rand();
        PROB=0.7;
        
        if ((x==1) + (y==1) + (x==maxX) + (y==maxY))~=0 %Particle kill re: boundaries
            die=1;
        else
            stuck=0;
            if p<=PROB && (x>=491 && x<=690) && (y>=312 && y<470)
                if (bw(xRight,y) ==1 || bw(xLeft,y) ==1 || bw(x,yUp) ==1 || bw(x,yDown) ==1)
                    fprintf('stuck: %d %d %f %f %f %f %d \n',x,y,bw(xRight,y),bw(xLeft,y),bw(x,yUp),bw(x,yDown), particle);
                    bw(x:x+7,y:y+7)=2; % Do not make it 1 so that other particle can still diffuse through this site
                    %Note that square grid is larger in the Irish Sea
                    %region than other transient regions
                    stuck=1;
                    
                    %Draw commands
                    imagesc(bw);
                    set(gca,'YDir','normal');
                    colormap(1-gray);
                    title(['Time:' num2str(tcountme) ' ' 'years']);
                    axis equal;
                    axis tight;
                    drawnow
                end
            elseif p>PROB
                if (bw(xRight,y) ==1 || bw(xLeft,y) ==1 || bw(x,yUp) ==1 || bw(x,yDown) ==1)
                    fprintf('stuck: %d %d %f %f %f %f %d \n',x,y,bw(xRight,y),bw(xLeft,y),bw(x,yUp),bw(x,yDown), particle);
                    %Attempt to avoid the border of the North of Ireland
                    if ( ((y>=237 && y<=269) && (x>=624 && x<=667)) || ((y>=269 && y<=288) && (x>=624 && x<=645)) || ((y>=204 && y<=241) && (x>=673 && x<=708)) || ((y>=219 && y<=236) && (x>=627 && x<=649))  )
                        continue;
                    else
                        bw(x:x+5,y:y+5)=2; % Do not make it one:other particle can still diffuse through this site
                        
                    end
                    stuck=1;
                    
                    %Draw commands
                    imagesc(bw);
                    set(gca,'YDir','normal');
                    colormap(1-gray);
                    title(['Time:' num2str(tcountme) ' ' 'years']);
                    axis equal;
                    axis tight;
                    drawnow
                    
                end
                
                
            end
            
        end
        
    end
    
    %Burkhard's print contributions
    fprintf('stuck: bw(%d,%d)=%f  -> %d %d   m: %d %d Iteration %d \n',x,y,bw(x,y),stuck,die,maxX,maxY, particle);
    
    %Draw commands
    imagesc(bw);
    set(gca,'YDir','normal');
    colormap(1-gray);
    title(['Time:' num2str(tcountme) ' ' 'years']);
    axis equal;
    axis tight;
    drawnow
    
    %{
Movie slide tracker/counter
    if mod(particle+9,10) == 0
        f=particle;
        saveas(gcf, ['pollutant',sprintf('%f',f),'.png'])
    end
%}    
    
end

for i = 1:size(bw,1)
    for j = 1:size(bw,2)
        if bw(i,j)==2
            msvcount=msvcount+1;
            
        end
    end
end
RADMEASURE = msvcount*eta;
fprintf('Radiation count in Irish sea model box: %d mSv',RADMEASURE)