code_part1

//Delphi 2007
// optimum is fixed 
unit MainUnit;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls, math, TeeProcs, TeEngine, Chart, Series;

const
INI_IND=2000;
NSTER=50;
LAMBDA=3;
MAX_IND=10000000;
NEQUILIBRIUM=1000;
NSIM=2101;
mutationrate_plas=0.001;
mutationrate_disp=0.001;
mutationrate_spec=0;
inioptimum=50;      // equals half the width of the climate window
CWS=1.3;
Nreplic=3;
steekproefgrootte=10000;     
breedte=100;
max_x=1550;

Type
  TForm1 = class(TForm)
    CreateLandscapeButton: TButton;
    Image1: TImage;
    RunButton: TButton;
    CloseButton: TButton;
    procedure CreateLandscapeButtonClick(Sender: TObject);
    procedure RunButtonClick(Sender: TObject);
    procedure CloseButtonClick(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

berekening=record
 dens:integer;      //total number of individuals
 som:double;          // sum of plasticity within cell
 end;

individual=record
 posx,posy:integer;
 plas:double;
 dispersalsd:double;
 spec:double
 end;

 TPoint=record
 x,y:integer;
 end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

uses CreateLandscapeUnit;

var
  time:integer;
  occupied:  array[1..max_x, 1..breedte] of boolean;      //whether pixel is occupied
  population: array [1..MAX_IND] of individual;
  young: array [1..MAX_IND] of individual;    // important voor event reproduction
  totalpop: integer;
  eigenschap_cel: array [1..max_x, 1..breedte] of berekening;
  survivors:array [1..MAX_IND]of individual;  // important for event survival
  competitors:array[1..MAX_IND] of individual;// important for event competition
  output:textfile;
  Topt:double;    //position of global optimum within landscape
  steekproef: array [1..10000]of individual;
  habitat:array[1..max_x,1..breedte]of integer;
  factor: integer;

/////**********************************************************************/////
/////                   Procedure to prepare landscape                     /////
/////**********************************************************************/////

Procedure preparelandscape;
var
x,y : integer;
begin
for x:=1 to max_x do
  for y:= 1 to breedte do
      begin
        occupied[x,y]:=false;
        habitat[x,y]:=1;
         //      case habitat[x,y] of
   //      0:  form1.Image1.Canvas.Pixels[x,y]:=clwhite;
   //      1:  form1.Image1.Canvas.Pixels[x,y]:=clwhite;
      end;
end;




/////**********************************************************************/////
/////                      Procedure to draw landscape                     /////
/////**********************************************************************/////

Procedure drawlandscape;
var
x,y: integer;
begin
for x:=1 to max_x do
  for y:= 1 to breedte do
      begin
        case habitat[x,y] of
           0:  form1.Image1.Canvas.Pixels[x,y]:=clwhite;
           1:  form1.Image1.Canvas.Pixels[x,y]:=clblack;
        end;
      end;
end;

//****************************************************************************//
//                     Procedure place optimum                                //
//****************************************************************************//
Procedure placeoptimum;
begin
  Topt:=inioptimum;
end;

//****************************************************************************//
//                    Procedure move optimum                                  //
//****************************************************************************//
procedure moveoptimum;
begin
  if (time>NEQUILIBRIUM) and (Topt<(max_x-inioptimum)) then Topt:=(Topt + cws);

end;
//------------------------------------------------------------------------------
//- function Poisson: function (level 4) called by Reproduction procedure      -
//------------------------------------------------------------------------------

function Poisson(lambda:real):integer;  // calculates the number of offsprings
                                        // of a single female from a Poisson
                                        // distribution using the variable
                                        // "lambda" (fertility) as mean value
var
  a,b      :real;
  i        :integer;
begin                                   // start algorithm from:
if  lambda>88 then  lambda:=88;         // Law, A., Kelton, D., Simulation
a:=exp(- lambda);                       // Modeling And Analysis, McGraw-Hill
b:=1;                                   // Book Company, New York, 1982,
i:=-1;                                  // p. 258 - 260
repeat
  i:=i+1;
  b:=b*random;
until (b<a)or (i>199);
poisson:=i;
end;

/////**********************************************************************/////
/////                   Procedure competition                              /////
/////**********************************************************************/////

Procedure competition;
var
Nr,b,a,c,x,y,d:integer;
numcompetitors:integer;
chance_of_mort, comp_coef,e:double;
begin

  for x := 1 to max_x do
  begin
    for y := 1 to breedte do
    begin
      eigenschap_cel[x,y].dens:=0;
      eigenschap_cel[x,y].som:=0;
      end;
      end;

  for d := 1 to totalpop do
  begin
    inc(eigenschap_cel[(population[d].posx),(population[d].posy)].dens);
    e:= (eigenschap_cel[(population[d].posx),(population[d].posy)].som);
    if ((exp((-1)*
    (sqr(population[d].posx-(Topt + population[d].spec))/(100 *population[d].plas))))*10*(exp((-factor)*population[d].plas)))>0 then
     eigenschap_cel[(population[d].posx),population[d].posy].som:=(e +((exp((-1)*
    (sqr(population[d].posx-(Topt + population[d].spec))/(100 *population[d].plas))))*10*(exp((-factor)*population[d].plas))));
    end;

  numcompetitors:=0;

  for Nr := 1 to totalpop do
    begin
      comp_coef:= ((exp((-1)*
    (sqr(population[Nr].posx-(Topt + population[Nr].spec))
    /(100 *population[Nr].plas))))*10*(exp((-factor)*population[Nr].plas)));
      if comp_coef<=0 then chance_of_mort:=1 else
      begin
      chance_of_mort:= (1-(NSTER/(
      (eigenschap_cel[(population[Nr].posx),(population[Nr].posy)].som)/comp_coef)));
      if eigenschap_cel[(population[Nr].posx),(population[Nr].posy)].som< comp_coef then
      writeln (eigenschap_cel[(population[Nr].posx),(population[Nr].posy)].som,
      comp_coef, eigenschap_cel[(population[Nr].posx),(population[Nr].posy)].dens:5);
      end;
      if random> chance_of_mort then
      begin
      inc (numcompetitors);
      competitors[numcompetitors] := population[Nr];
    //  form1.image1.Canvas.pixels[young[i].posx,young[i].posy]:=clyellow;
      end;
      end;

  for a := 1 to totalpop do
    begin
    population[a].posx:=0;
    population[a].posy:=0;
    population[a].dispersalsd:=0;
    population[a].plas:=0;
    population[a].spec:=0;
    end;

    totalpop:=numcompetitors;

  for b := 1 to totalpop do
    begin
    population[b]:=competitors[b];
    end;

  for c := 1 to numcompetitors do
    begin
    competitors[c].posx:=0;
    competitors[c].posy:=0;
    competitors[c].dispersalsd:=0;
    competitors[c].plas:=0;
    competitors[c].spec:=0;
    end;


end;


/////**********************************************************************/////
/////                   Procedure reproduction                             /////
/////**********************************************************************/////

Procedure reproduction;
var
Nr,n,b:integer;
Newpop:integer;
offspring, populatie_ouders:integer; //offspring: absolute number of offspring per adult
mutation_spec, mutation_plas, mutation_disp:double;

begin
  newpop:=0;
  offspring:=0;
  populatie_ouders:=totalpop;

  for Nr := 1 to totalpop do
  begin
  offspring:= Poisson(LAMBDA);
  if offspring>0 then
  begin
  for n := 1 to offspring do

  begin
  young[newpop+n]:= population[Nr];
  if random<mutationrate_plas then mutation_plas := (random-0.5)else mutation_plas:=0;
  if ((young[newpop+n].plas) + mutation_plas) <=0.00000000001 then mutation_plas:= ((-1)*mutation_plas) ;
  young[newpop+n].plas:=population[Nr].plas + mutation_plas;

  if random<mutationrate_disp then mutation_disp:= (random-0.5)else mutation_disp:=0;
  if ((young[newpop+n].dispersalsd) + mutation_disp) <0 then mutation_disp:= ((-1)*mutation_disp);
  young[newpop+n].dispersalsd :=population[Nr].dispersalsd + mutation_disp;

  if random<mutationrate_spec then mutation_spec:= ((random*10)-5)else mutation_spec:=0;
  young[newpop+n].spec:=population[Nr].spec + mutation_spec;
  end;

  newpop:=newpop+offspring;
  end;
  end;

  totalpop:=newpop;

  for b := 1 to populatie_ouders do
  begin
  population[b].posx:=0;
  population[b].posy:=0;
  population[b].dispersalsd:=0;
  population[b].plas:=0;
  population[b].spec:=0;
  end;
  
end;

//continue working with array young//



/////**********************************************************************/////
/////                   Function for Gaussian dispersal                    /////
/////**********************************************************************/////
function DisperseGauss(location:TPoint;stratDist:double):TPoint;
var
dx,dy:extended;
begin
  dy:=randG(0,stratdist);
  dx:=randG(0,stratdist);
  location.x:=round(location.x+dx);
  location.y:=round(location.y+dy);
  DisperseGauss:=location;
end;
/////**********************************************************************/////
/////                   Procedure to move to other point                   /////
/////**********************************************************************/////

procedure move;
var
nr:integer;
location,newlocation:TPoint;
a:double;

begin
for Nr := 1 to totalpop do
begin
location.x:=young[Nr].posx;
location.y:=young[Nr].posy;
a:= young[Nr].dispersalsd;
newlocation:=DisperseGauss(location,a);
young[nr].posx:=newlocation.x;
young[nr].posy:=newlocation.y;
if young[nr].posy>breedte then young[nr].posy:= (young[nr].posy - breedte);
if young[nr].posy<=0 then young[nr].posy:= (young[nr].posy + breedte);
end;
    end;

/////**********************************************************************/////
/////                   Procedure to survive                               /////
/////**********************************************************************/////

procedure survive;

var
a,b,k,l:integer;
numsurvivors, populatie_zaden:integer;
begin

  numsurvivors:=0;
  populatie_zaden:=totalpop;

  for a := 1 to totalpop do
    begin
      if (habitat[young[a].posx,young[a].posy]=1)
      and (round((Topt-inioptimum+1))<=young[a].posx)
      and (young[a].posx<=round(Topt+inioptimum))
      and (1<=young[a].posy)
      and (young[a].posy<=breedte)
      then
      begin
      inc (numsurvivors);
      survivors [numsurvivors] := young [a];
    //  form1.image1.Canvas.pixels[young[i].posx,young[i].posy]:=clyellow;
      end;
      end;

      totalpop:=numsurvivors;

      for b := 1 to totalpop do
        begin
        population[b]:=survivors[b];
        end;

       for l :=1 to populatie_zaden do
       begin
        young[l].posx:=0;
        young[l].posy:=0;
        young[l].dispersalsd:=0;
        young[l].plas:=0;
        young[l].spec:=0;
       end;

      for k := 1 to numsurvivors do
        begin
          survivors[k].posx:=0;
          survivors[k].posy:=0;
          survivors[k].dispersalsd:=0;
          survivors[k].plas:=0;
          survivors[k].spec:=0;
        end;

end;


/////**********************************************************************/////
/////                   Procedure to interupt simulation                   /////
/////**********************************************************************/////

procedure TForm1.CloseButtonClick(Sender: TObject);
begin
   halt(0);
   close;  // exits the program
end;

/////**********************************************************************/////
/////                    Procedure to initialise                           /////
/////**********************************************************************/////

procedure initialise;
var x,y: integer;
    i: integer;

     begin
      totalpop:=0;

        for i := 1 to INI_IND do
      begin
        begin
        Repeat
          begin
          x:= round(random(100)+1);      // random number is drawn within interval[0,99]
          y:= round(random(breedte)+1);
          population[i].posx:=x;
          population[i].posy:=y;
          end;
        Until
         (habitat[x,y]=1)
          and (occupied[population[i].posx,population[i].posy]= false);
        end;
        occupied[population[i].posx,population[i].posy]:=true;
        if (factor=1) or (factor=5) then population[i].plas:=(random*6+0.00000000000001);
        if factor=10 then population[i].plas:=(random*2+0.00000000000001);
        population[i].dispersalsd:=(5*random);
        population[i].spec:=0;//((random*(2*inioptimum))-50);
  //    form1.image1.Canvas.pixels[population[i].posx,population[i].posy]:=clyellow;
        inc(totalpop);
        end;
         //densloc array aanmaken
end;
/////**********************************************************************/////
/////                    Procedure to create landscape                     /////
/////**********************************************************************/////

procedure TForm1.CreateLandscapeButtonClick(Sender: TObject);
begin
CreateLandscapeUnit.createlandscape;
//drawlandscape;
writeln('p habitat = ',p_habitat:3:2);
end;

/////**********************************************************************/////
/////                          Procedure nemensteekproef                   /////
/////**********************************************************************/////

Procedure nemensteekproef;
var
i,a:integer;
begin
if totalpop>=100 then
begin
for i := 1 to steekproefgrootte do
begin
  repeat a:=round(random*totalpop) until a>0 ;
  steekproef[i]:=population[a];
end;
end;
end;

/////**********************************************************************/////
/////                          Main programme                              /////
/////**********************************************************************/////

procedure TForm1.RunButtonClick(Sender: TObject);
var
i,N,a:integer;
begin
assignfile(output,'results.txt');
rewrite(output);
for N := 1 to Nreplic do
begin
  factor:=N;
  if factor=2 then factor:=5;
  if factor=3 then factor:=10;  
  for a := 1 to 5 do
   begin
  //CreateLandscapeUnit.createlandscape;
 preparelandscape;
 //drawlandscape;
 placeoptimum;
 initialise;
 for time := 1 to NSIM do
 begin
 for i := 1 to steekproefgrootte do
 begin
   steekproef[i].plas:=0;
   steekproef[i].posx:=0;
   steekproef[i].posy:=0;
   steekproef[i].dispersalsd:=0;
   steekproef[i].spec:=0;
 end;
 writeln('generation',time:6);
  competition;
  if (time = 1000) or (time=2000) or (time=2100) then
  begin
  nemensteekproef;
  for i := 1 to steekproefgrootte do
  begin
  writeln(output,time:6, factor:6, totalpop:9, steekproef[i].dispersalsd:10:6,steekproef[i].plas:12:8,steekproef[i].posx:6, steekproef[i].spec:8:4);
  end;
  end;
  reproduction;
  move;
  survive;
  moveoptimum;
//writeln(output, time:6, totalpop:6);
application.ProcessMessages;
 end;
end;
end;
closefile(output);
writeln('Finished!');

end;

end.