function  [FrontUpLoc,FrontDwnLoc,DailyMeanBiomass,Crayfish_Alive,Param] = Crayfish_dispersal_v2(Param)
% Function name: Crayfish_dispersal_v2
% Description: Contains the crayfish dispersal model
% Usage: [FrontUpLoc,FrontDwnLoc,DailyMeanBiomass,Crayfish_Alive,Param] = Crayfish_dispersal_v2(Param)
% Inputs: Param
% Outputs: FrontUpLoc FrontDwnLoc DailyMeanBiomass Crayfish_Alive Param
% Author: Dr Jim Kerr
% Last Modified: May 2020

%% Predefines river length cells and mid-points
River.Length=Param.RiverLength;
binwidth=10; % calculates the bin width for the biomass and density calculations
River.Width(1:River.Length/binwidth)=2;
xbins=(0:binwidth:River.Length); % Calculates bins along river length
Location=xbins(1:end-1)+(diff(xbins/2)); % Calculate midpoint of bin cell for later interpolation
Reproduction_xbins=(0:100:River.Length); % Calculates bins along river length
TemperatureScalingfactor=((((sin(linspace(-0.5*pi,1.5*pi,365))+1)/2)/1.11111111)+0.1);

%% Setup model parameters
Startdate=Param.Startdate;
ModelDuration=Param.ModelDuration; 
NoCray=Param.NoCray; 
MeanAgeCray=Param.MeanAgeCray; 
StdAgeCray=Param.StdAgeCray; 
MeanBreedingPeriod=Param.MeanBreedingPeriod; 
RangeBreedingPeriod=Param.RangeBreedingPeriod; 
SexuallyMatureLength=Param.SexuallyMatureLength; 
Gestation_Time=Param.Gestation_Time; 
StartPosition=Param.StartPosition;
Plot=Param.Plot;
Obstruct=Param.Obstruct;
BarrierLoc=Param.BarrierLoc;

%% Setup the crayfish data
Cray.ID=1:NoCray;
Cray.Age=zeros(1,NoCray);
Cray.Length=zeros(1,NoCray);
Cray.Sex=round(rand(1,NoCray)); % 0 = males, 1 = females
Cray.Location=zeros(1,NoCray);
Cray.Movement=zeros(1,NoCray);
Cray.Pregnant=zeros(1,NoCray);
Cray.Gestation=zeros(1,NoCray);
Cray.Mortality=zeros(1,NoCray);
Cray.Weight=zeros(1,NoCray);
Cray.Density=zeros(1,NoCray);
Cray.Biomass=zeros(1,NoCray);
Cray.Location=randn(1,NoCray)*20+StartPosition;
Cray.UpDwn=randn(1,NoCray);
Cray.Passer=rand(1,NoCray);

FrontUpLoc(1:ModelDuration*365)=NaN;
FrontDwnLoc(1:ModelDuration*365)=NaN;
DailyMeanBiomass(1:ModelDuration*365)=NaN;
Crayfish_Alive(1:ModelDuration*365)=NaN;
No_Crayfish_Pregnant(1:ModelDuration*365)=NaN;

%% Allocate startup information to crayfish data
Cray.Age=(randn(1,NoCray))*(StdAgeCray*365)+(MeanAgeCray*365); 
Cray.Length=23.669*((Cray.Age/365).^0.614); % Length is calculated from the relationships presented in Guan and Wiles, 1999 (see Crayfish Info.xlsx).
Cray.Weight=0.00024*(Cray.Length).^3.05263; % Weight is calculated from the relationships presented in Guan and Wiles, 1999 (see Crayfish Info.xlsx).

%% Setup the figure if specified
if Plot==1
    figure(1)
    hold on
    xlabel('Longitudinal distance (m)')
    ylabel('Biomass (g m^2)')
    xlim([0 River.Length]);
    ylim([0 100])
    ax = gca;
    ax.XRuler.Exponent = 0;
    if Param.Video==1
        v = VideoWriter(strcat(Param.RunNumber,'.avi'));
        open(v);
    end
end

%% Iteratively run the model
for j=1:ModelDuration*365

    %% Increment all the time dependent factors
    JulianDay=rem(Startdate+j,365); % Calculate the julian day
    Year=floor((Startdate+j)/365);  % Calculate the year
    Cray.Gestation(Cray.Pregnant==1)=Cray.Gestation(Cray.Pregnant==1)+1; % Increase gestation progress for all pregnant crayfish
    Cray.Age=Cray.Age+1; % Increase the age of the Crayfish by one day
    Cray.Length=23.669*((Cray.Age/365).^0.614); % Calculate the size of the crayfish (CL mm)
    Cray.Weight=0.00024*(Cray.Length).^3.05299; % Calculate the weight of the crayfish (g)
    
    %% Output descriptive data
    Crayfish_Alive(j)=length(Cray.ID); % Calculate the number of crayfish that are alive
    No_Crayfish_Pregnant(j)=sum(Cray.Pregnant); % Calculate the number of crayfish that are pregnant
   
    %% Calculates the density (number of individuals) of crayfish in same river section for each crayfish
    [Number,~,idx]=histcounts(Cray.Location,xbins);
    Cray.Density=interp1(Location,Number,Cray.Location,'nearest'); % Calcualte the density of crayfish surrounding each crayfish (for density dependent mortality)
    
    %% Calculates the biomass of crayfish in same river section for each crayfish
    Biomass=(accumarray(idx',Cray.Weight'))';
    Biomass(end+1:length(Location))=0; % extends the array with zeros to cover whole length of river
    NormBiomass=(Biomass./River.Width)/binwidth; % grams per metre squared of river.
    Cray.Biomass=interp1(Location,NormBiomass,Cray.Location,'nearest');
    
    %% Move the crayfish
    T=rand(1,Crayfish_Alive(j)); % preallocate T to be random numbers between 0 and 1
    Q=zeros(1,Crayfish_Alive(j)); % preallocate Q to be all zeros
    ScaleFactor=ones(1,Crayfish_Alive(j)); % preallocate the ScaleFactor variable
    Q(Cray.UpDwn< 0)=-rand(1,sum(Cray.UpDwn< 0)); % randomly allocate the downstream movers a value between -1 and 0
    Q(Cray.UpDwn>=0)= rand(1,sum(Cray.UpDwn>=0)); % randomly allocate the upstream movers a value between 0 and 1
    Q(T<=0.2)=(rand(1,sum(T<=0.2))-0.5)*2; % randomly assigns 20% of all movement to be between -1 and 1.
    Q(Q<0)=-(((exp((3.161-log(abs(Q(Q<0))))/1.236))-12.902787463944376)/2); % allocate all downstream movements a length according to the Bubb et al. (2004) movement distributions.
    Q(Q>0)=+(((exp((3.845-log(abs(Q(Q>0))))/1.464))-13.823445837555905)/2); % allocate all upstream movements a length according to the Bubb et al. (2004) movement distributions.
    Q(Q<-170)=rand(1,length(Q(Q<-170)))*-170; % limit the downstream movements according to the Bubb et al. (2004) maximum 
    Q(Q> +90)=rand(1,length(Q(Q> +90)))* +90; % limit the upstream movements according to the Bubb et al. (2004) maximum 
    ScaleFactor(Cray.Biomass<70.33)=Cray.Biomass(Cray.Biomass<70.33)/70.33; % calculate the density dependent scaling factor. Carrying capacity is set at 70.33 based on mean of values from Guan (2000) and Chadwick (2019).
    ScaleFactor=ScaleFactor*TemperatureScalingfactor(JulianDay+1); % scale the scaling factor by to be temperature dependent sinuosoidal variance from 0.1 at start of year 1 mid year back to 0.1 at end of year.
    Cray.Movement=Q.*ScaleFactor; % apply the density dependent scaling factor
    TempLoc=Cray.Location+Cray.Movement;
    
    if Obstruct==1 && max(Cray.Location)>min(BarrierLoc)-300
        for i=1:length(BarrierLoc)
            idx=(BarrierLoc(i)>Cray.Location & BarrierLoc(i)<=TempLoc & or(Cray.Passer>0.22 & Cray.Sex==0,Cray.Passer>0.12 & Cray.Sex==1));
            if sum(idx)>0
                TempLoc(idx)=BarrierLoc(i)-(rand*20);
            end
        end
    end
    
    Cray.Location=TempLoc; % move the crayfish
    
    %% Calculate the probablity of mortality for each crayfish
    AnnualMort=-0.26*log(Cray.Weight.*0.1651)+0.91; % General invertebrate annual mortality rates given by McCoy and Gillooly, 2008. Uses weight to dry weight conversion factor from Thompson et al. 2004.
    AnnualMort=AnnualMort/100; % making sure values are less than 1.
    DailyMort=1-(1-AnnualMort).^(1/365); % Calculating daily mortality rates
    DailyMort=DailyMort*100; % Correct for earlier correction to AnnualMort variable
    ScaleFactor=1+(0.0072*Cray.Biomass).^10;
    Cray.Mortality=1-((1-DailyMort)./ScaleFactor);

    %% Kill off the specified crayfish individuals
    Q=rand(1,Crayfish_Alive(j));
    idx = Q < Cray.Mortality | Cray.Location < 0  | Cray.Location > River.Length | Cray.Age>7*364;
    Cray.ID(idx)=[];
    Cray.Age(idx)=[];
    Cray.Length(idx)=[];
    Cray.Sex(idx)=[];
    Cray.Location(idx)=[];
    Cray.Movement(idx)=[];
    Cray.Pregnant(idx)=[];
    Cray.Gestation(idx)=[];
    Cray.Mortality(idx)=[];
    Cray.Weight(idx)=[];
    Cray.Density(idx)=[];
    Cray.Biomass(idx)=[];
    Cray.UpDwn(idx)=[];
    Cray.Passer(idx)=[];

    %% Make female crayfish pregnant if they are large enough, its the right time of year and they are near a male
    if JulianDay>MeanBreedingPeriod-RangeBreedingPeriod && JulianDay<MeanBreedingPeriod+RangeBreedingPeriod % if within the breeding period
        idx = Cray.Sex==0 & Cray.Length>SexuallyMatureLength; % finds males of breeding size
        [Number_males,Location_males]=hist(Cray.Location(idx),Reproduction_xbins); % finds number of males in river 100m section
        idx = Cray.Sex==1 & Cray.Length>SexuallyMatureLength & interp1(Location_males,Number_males,Cray.Location,'nearest')>0 & Cray.Pregnant~=1; % finds females of breeding size that are in a river section with at least one male and which are not already pregnant
        Cray.Pregnant(idx)=1;
    end

    %% Make female crayfish give birth at end of gestation time
    idxloc = find(Cray.Gestation>Gestation_Time);
    if ~isempty(idxloc)
        Fecundity=round(0.0436*Cray.Length(idxloc).^2.2233); % from Capurra et al., 2015. relationship between potential fecundity and size.
        IDX=zeros(1,sum(Fecundity));
        count=1;
        for i=1:length(idxloc)
            IDX(count:Fecundity(i)+count)=idxloc(i);
            count=count+Fecundity(i);
        end
        Cray.ID(end+1:end+length(IDX))=NoCray+1:NoCray+length(IDX);
        Cray.Age(end+1:end+length(IDX))=1;
        Cray.Length(end+1:end+length(IDX))=23.669*(1/365)^0.614;
        Cray.Weight(end+1:end+length(IDX))=0.00024*(0.623)^3.05299;
        Cray.Sex(end+1:end+length(IDX))=round(rand(1,length(IDX)));
        Cray.Location(end+1:end+length(IDX))=Cray.Location(IDX);
        Cray.Pregnant(end+1:end+length(IDX))=zeros(1,length(IDX));
        Cray.Gestation(end+1:end+length(IDX))=zeros(1,length(IDX));
        Cray.Mortality(end+1:end+length(IDX))=NaN;
        Cray.Density(end+1:end+length(IDX))=NaN;
        Cray.Biomass(end+1:end+length(IDX))=NaN;
        Cray.UpDwn(end+1:end+length(IDX))=(randn(1,length(IDX)))/4;
        Cray.Movement(end+1:end+length(IDX))=NaN;
        NoCray=NoCray+length(IDX);
        Cray.Pregnant(idxloc)=0;
        Cray.Gestation(idxloc)=0;
        Cray.Passer(end+1:end+length(IDX))=rand(1,length(IDX));
    end

    %% Plot figure if specified
    if Plot==1 && rem(j,20)==0
        h=bar(Location,NormBiomass');
        hold on
        title(['Year: ',num2str(Year),', Julian Day: ',num2str(JulianDay),', CA: ',num2str(Crayfish_Alive(j)),', CP: ',num2str(No_Crayfish_Pregnant(j))]);
        if Obstruct==1
            vline(BarrierLoc)
        end
        drawnow
        if Param.Video==1
            frame = getframe(gcf);
            writeVideo(v,frame);
        end
        delete(h)
    end
 
    FrontUpLoc(j)=Location(find(Biomass>max(Biomass)/4,1,'last'));
    FrontDwnLoc(j)=Location(find(Biomass>max(Biomass)/4,1,'first'));
    DailyMeanBiomass(j)=mean(NormBiomass(round(FrontDwnLoc(j)/10):round(FrontUpLoc(j)/10)));
end

if Param.Video==1
    close(v);
end

%% Plot results after model has completed
% figure(1)
% yyaxis left
% bar(Location,NormBiomass)
% title(strcat('Year:',num2str(Year),'Julian Day:',num2str(JulianDay)));
% xlim([0 River.Length]);
% yyaxis right
% plot(Location,River.Width)
% 
% figure(2)
% X=Startdate:j+Startdate-1;
% plot(X,FrontUpLoc)
% hold on
% plot(X,FrontDwnLoc)
% xlim([Startdate j+Startdate])
% 
% times=10:365:length(X)-10;
% for i=1:length(times)-1
%     UpRate(i)=mean(FrontUpLoc(times(i+1)-9:times(i+1)+9))-mean(FrontUpLoc(times(i)-9:times(i)+9));
%     DwnRate(i)=mean(FrontDwnLoc(times(i)-9:times(i)+9))-mean(FrontDwnLoc(times(i+1)-9:times(i+1)+9));
% end
% 
% figure(3)
% plot(0:length(times)-2,UpRate)
% hold on
% plot(0:length(times)-2,DwnRate)
% 
% figure(4)
% histogram(Cray.Movement)
% 
% figure(5)
% histogram(Cray.Age/365,200)
% 
% figure(6)
% plot(No_Crayfish_Pregnant)
% 
% figure(7)
% plot(Crayfish_Alive)