%
%
%Model of spreading of Fouth and Mouth disease and tool to simulate control
%IDE course of Jim Powel, April 2001
%Student: A.G.T. (Tom) Schut
%
%


np=128; 
xl=10; dx=2*xl/np;
yl=10; dy=2*yl/np;
x=linspace(-xl,xl-dx,np);
y=linspace(-yl,yl-dy,np)
[X,Y]=meshgrid(x,y);

b=1.0;																	%dispersion parameter
c=1.01;																	%dispersion parameter
fs=0.8;																	%fraction "seed"dispersal
Lfs=60;																	%Level of virus were diffusion starts

IFprob=0.8;																%maximum chance of infection
radius=1;																%radius of cleared farms
Vmin=.005;																%minimum amount of virus needed for survival

Rfarm=0.5;																%radius of farms
days=15;																   %number of hours before control measure is taken
rate=5.0; 																%logistic growth rate per day of virus;
IFinit=25;																%initial infection level
IFK=500;																	%Carrying capacity 
Control=100;															%level were there is control measure taken

%Define farm area's
Farms=( rand(np,np)>0.95 );										%select random points where farms are located

Kf=(sqrt(X.^2+Y.^2)<Rfarm);										%define farm size kernel with size Rfarm
Kf=Kf / trapz(trapz(dx*dy*Kf) );								   %normalize
Kf=dx*dy*fft2(Kf);													%take fourier transform

Farms=( abs(fftshift(ifft2(Kf.*fft2(Farms)))) >0.1 );		%Make farms of certain size
IFfarm=(X>0 & X<Rfarm) & (Y>0 & Y<Rfarm);						%make sure that infected farm is present
Farms=( (Farms+IFfarm)>0.5 );										%farm area's have value 1, rest 0

figure(3),surf(X,Y,Farms), colormap(jet), shading interp;
pause(.1);


%dispersion kernel
R=sqrt(X.^2+Y.^2);
																			%"seed" dispersal kernel
Ks=(1/2) * c *( (abs(R).^(c-1) )/b^c) .* exp(-(abs(R)/b).^c);    
Ks=Ks / trapz(trapz(dx*dy*Ks));									%normalize, as it is a probability function
Ks=Ks.*IFprob;															%maximize to max IF probability
FKs=dx*dy*fft2(Ks);													%take fourier transform

%clearing kernel, remove farms if above threshold level
Kcm=(sqrt(X.^2+Y.^2)<radius);										%define clearing kernel with size radius
Kcm=Kcm / trapz(trapz(dx*dy*Kcm) );								%normalize
FKcm=dx*dy*fft2(Kcm);												%take fourier transform


%define initial conditions
V0=IFinit *( (X>0 & X<0.5) & (Y>0 & Y<0.5) ) ;				%define initially infected area
 
CM=V0.*0;																%define last control measure map
TCM=V0.*0;																%define total control measure map

for j=1:days
   V=V0.*exp(rate*(1-V0/IFK));									%logistic growth of virus population
   V=V.*Farms;															%growth only on farms
    
   																		%disperse if above threshold Lfs
   V0=(1-fs)*V+fs*real(ifftshift(ifft2(fft2(V>=Lfs).*FKs)));   

	VM=(V0>Vmin);														%Define map where virusses are above minimum
   
	V0=V0.*VM;															%area's where virusses are above minimum
 
   CM=abs(fftshift(ifft2(FKcm.*fft2((V0>=Control)))));	%area where control measures are taken
   cmax=max(max(CM));
   TCM=(TCM>.5)+(CM>.5*cmax);										%Total area of control
   
   NC=(TCM<0.5);														%area where no control measure is taken
   
   V0=V0.*NC;															%remove CM area's from virus map   
 	figure(1),surf(X,Y,V0), colormap(jet), shading interp,legend('Virus infection');
   pause(.15);
   figure(2),surf(X,Y,TCM),colormap(jet), shading interp,legend('Total control measures');   
   pause(.15);
end

