Index: /experiments/frams/foraminifera/data/scripts/foraminifera.expdef
===================================================================
--- /experiments/frams/foraminifera/data/scripts/foraminifera.expdef	(revision 478)
+++ /experiments/frams/foraminifera/data/scripts/foraminifera.expdef	(revision 479)
@@ -23,13 +23,12 @@
 code:~
 
-
+global chambers;
+global colors;
+global dir_change;
+global max_chamber_volume;
+global movePerStep;
 global nutrientenergywaiting;
+global o;
 global reprocounter;
-global colors;
-global chambers;
-global o;
-global max_chamber_energ;
-global dir_change;
-global movePerStep;
 
 @include "foraminifera.inc"
@@ -83,5 +82,5 @@
 
 	//time
-	ExpParams.secPerStep = 60;
+	ExpParams.secPerStep = 480;
 	ExpParams.foramSpeedMmPerMin = 0.05;
 	movePerStep = getMovePerStep();
@@ -92,9 +91,13 @@
 
 	//reproduction
-	ExpParams.foramPop = 10;	
+	ExpParams.foramPop = 4;	
 	ExpParams.crossprob = 0;
 	ExpParams.mutationprob = 0;
-	ExpParams.repro_time = 20;
+	ExpParams.repro_time = 720;
+	ExpParams.gametoPeriod = 43200;
+	ExpParams.divisionCost = 15.6;
 	reprocounter = 0;
+
+	init_chambers();
 
 	//morphology
@@ -102,15 +105,17 @@
 	ExpParams.zone1_range = micronsToFrams(1000);
 	ExpParams.zone2_range = micronsToFrams(3000);
-	init_chambers();
+	ExpParams.chamber_growth_time = 720;
+	ExpParams.chamberCostPerSec = 0.000001;
 	ExpParams.chamber_proculus_haplo = micronsToFrams(50);
 	ExpParams.chamber_difference_haplo = 0.0;
 	ExpParams.chamber_proculus_diplo = micronsToFrams(20);
 	ExpParams.chamber_difference_diplo = 0.2;
-	max_chamber_energ = [Vector.new(), Vector.new()];
+
+	max_chamber_volume = [Vector.new(), Vector.new()];
 	for (var j = 0; j < 2; j++)
 	{
 		for (var i = 0; i < chambers[0].size; i++)
 		{
-			max_chamber_energ[j].add(((energyFromVolume(getProperty(j, "chamber_proculus")) + energyFromVolume(getProperty(j, "chamber_proculus") + (i) * getProperty(j, "chamber_difference")))*(i+1))/2);
+			max_chamber_volume[j].add(((volumeInMicrons(getProperty(j, "chamber_proculus")) + volumeInMicrons(getProperty(j, "chamber_proculus") + (i) * getProperty(j, "chamber_difference")))*(i+1))/2);	
 		}                                  
 	}
@@ -120,19 +125,23 @@
 	ExpParams.min_repro_energ_diplo = 6;
 
-	ExpParams.e_meta = 0.00005; 
-	ExpParams.energy_hib = 0.000025;
-	ExpParams.energy_move = 0.00005;
-
-	ExpParams.energies0_haplo = max_chamber_energ[0][0] - 0.001*max_chamber_energ[0][0];
-	ExpParams.energies0_diplo = max_chamber_energ[1][0] - 0.001*max_chamber_energ[1][0];
-	ExpParams.feedtrans = 0.125;
-	ExpParams.e_repro_cost_haplo = 0.7; 
-	ExpParams.e_repro_cost_diplo = 0.3;
+	ExpParams.e_meta = 0.0000005; 
+	ExpParams.energy_hib = 0.0000001;
+	ExpParams.energy_move = 0.0000005;
+
+	ExpParams.energies0_haplo = energyFromVolume(ExpParams.chamber_proculus_haplo,1);
+	ExpParams.energies0_diplo = energyFromVolume(micronsToFrams(1.25),1);
+	ExpParams.feedtrans = 0.001;
+	ExpParams.e_repro_cost_haplo = 0.3; 
+	ExpParams.e_repro_cost_diplo = 0.2;
+
+	ExpParams.e_death_level_haplo = 0.5;
+	ExpParams.e_death_level_diplo = 0.5;
 
 	//nutrients
-	ExpParams.nutrientsize = micronsToFrams(10);
-	ExpParams.energy_nut = 100 * energyFromVolume(ExpParams.nutrientsize);
-	ExpParams.nutrientPop = 1;
-	ExpParams.feedrate = 100;
+	ExpParams.foodperiod = 480;
+	ExpParams.nutrientradius = micronsToFrams(10);
+	ExpParams.energy_nut = energyFromVolume(ExpParams.nutrientradius,1);
+	ExpParams.feedrate = 1;
+	ExpParams.ingestion = 0.25;
 	nutrientenergywaiting = 0;
 	ExpState.totaltestedcr = 0;
@@ -142,7 +151,19 @@
 @include "standard_placement.inc"
 
-function energyFromVolume(radius)
-{
-	return 4.0/3.0*Math.pi*Math.pow(radius,3);
+function volumeInMicrons(radiusInFrams)
+{
+	return 4.0/3.0*Math.pi*Math.pow(framsToMicrons(radiusInFrams),3);
+}
+
+function energyFromVolume(base, isRadiusInFrams)
+{
+	if (isRadiusInFrams == 1) //radius in frams
+	{
+		return ExpParams.picoCarbonPerMikro*volumeInMicrons(base);
+	}
+	else //volume in microns
+	{
+		return ExpParams.picoCarbonPerMikro * base;
+	}
 }
 
@@ -169,13 +190,29 @@
 }
 
-function addInitialForam(species, i)
-{
-	var geno = createForamGenotype(0, species, 0);
+function addForam(species, iter, chambernum, ploid)
+{
+	var geno = createForamGenotype(ploid, species, chambernum);
 	var cr = Populations[0].add(geno);
-	cr.name = "Initial creature" + species + "_" + i;
+	cr.name = "Initial creature" + species + "_" + iter;
 	placeCreatureRandomly(cr, 0, 0);
-	cr.energy0 = getProperty(0, "energies0");
+	cr.energy0 = energyFromVolume(max_chamber_volume[ploid][chambernum],0);
 	cr.energy = cr.energy0;
 	setGenotype({"opt" : 0, "cr" : cr, "species" : species});
+	if (ploid == 1)
+	{
+		cr.data->lifeparams->gen = 1;
+		cr.data->genes = [cr.data->genes, cr.data->genes]; //TODO two different genes sets
+	}
+}
+
+function addInitialForam(species, iter)
+{
+	var ploid = 0;
+	if (Math.rnd01 > 0.5)
+	{
+		ploid = 1;
+	}	
+	//add new foram with random energy bewtween starting energy and reproduction treshold
+	addForam(species, iter, Math.rndUni(0,getProperty(ploid, "min_repro_energ")),ploid);
 }
 
@@ -261,7 +298,7 @@
 // -------------------------------- foram begin -----------------------------------
 
-function setForamMeta(cr, gen)
-{
-	//cr.idleen = (ExpParams.e_meta * max_chamber_energ[gen][Math.min(lastChamberNum(cr), max_chamber_energ[gen].size-1)])*ExpParams.secPerStep;
+function setForamMeta(cr)
+{
+	//percent of current energy
 	cr.idleen = (ExpParams.e_meta * cr.energy)*ExpParams.secPerStep;
 }
@@ -272,10 +309,15 @@
 }
 
+function getZoneRange(cr, zone_num)
+{
+	return ExpParams.["zone"+zone_num+"_range"];
+}
+
 function onForamsBorn(cr)
 {
-	setForamMeta(cr, 1);
+	setForamMeta(cr);
 	if (ExpParams.visualize == 1)
 	{
-		var ret = Populations[2].add("//0\np:sh=3,sx=0.01,sy="+ExpParams.zone1_range+",sz="+ExpParams.zone1_range+",ry=1.57,vr=1.0,1.0,1.0");
+		var ret = Populations[2].add("//0\np:sh=3,sx=0.01,sy="+getZoneRange(cr,1)+",sz="+getZoneRange(cr,1)+",ry=1.57,vr=1.0,1.0,1.0");
 		cr.data->reticulopodiacreature = ret;
 	}
@@ -322,5 +364,5 @@
 		setGenotype({"cr" : cr2, "parent_genes" : cr.data->genes, "parent_lifeparams" : cr.data->lifeparams, "opt" : 2});
 		cr2.moveAbs(cr.center_x - cr2.size_x / 2, cr.center_y - cr2.size_y / 2, cr.pos_z);
-		setForamMeta(cr2, cr2.data->lifeparams->gen);
+		setForamMeta(cr2);
 
 		if (visualization(cr))
@@ -335,5 +377,5 @@
 {
 	var p = cr.getMechPart(0);
-	var n = cr.signals.receiveSet("nutrient", ExpParams.zone2_range);
+	var n = cr.signals.receiveSet("nutrient", getZoneRange(cr,2));
 
 	//if signals are received find the source of the nearest
@@ -354,5 +396,5 @@
 			distvec.sub(p.pos);
 			dist = distvec.length;
-			if (dist < ExpParams.zone1_range)
+			if (dist < getZoneRange(cr,1))
 			{
 				if (n[i].value != null)
@@ -390,6 +432,6 @@
 	if (cr.data->lifeparams->hibernated == 0)
 	{
-		//cr.energy_m += (ExpParams.energy_move * max_chamber_energ[cr.data->lifeparams->gen][Math.min(lastChamberNum(cr), (max_chamber_energ[cr.data->lifeparams->gen].size)-1)])*ExpParams.secPerStep;
-		cr.energy_m += (ExpParams.energy_move * cr.energy)*ExpParams.secPerStep;
+		//percent of maximal energy
+		cr.energy_m += (ExpParams.energy_move * cr.data->lifeparams->max_energy_level)*ExpParams.secPerStep;
 	}
 }
@@ -422,5 +464,5 @@
 	}
 	//random move
-	else if (cr.lifespan%(dir_change/ExpParams.secPerStep) == 0)
+	else if (cr.lifespan%int(dir_change/ExpParams.secPerStep) == 0)
 	{
 		cr.data->lifeparams->dir = randomDir();
@@ -445,9 +487,10 @@
 {
 	cr1.localDrive = XYZ.new(0,0,0);
-	var e = ExpParams.feedtrans*ExpParams.secPerStep; //TODO efficiency dependent on age
-	e = Math.min(cr2.energy, e) + 0.0000001;
-	//Simulator.print("transferring "+e+" to "+cr1.name+" from "+cr2.name+" ("+cr2.energy+")");
-	cr2.energy_m = cr2.energy_m + e;
-	cr1.energy_p = cr1.energy_p + e;
+	var e =  ExpParams.feedtrans*cr1.energy;//ExpParams.feedtrans*cr1.energy*ExpParams.ingestion*ExpParams.secPerStep; //TODO efficiency dependent on age
+	e = Math.min(cr2.energy, e);
+	e *= ExpParams.secPerStep;
+	//Simulator.print("transferring "+e +"("+e*ExpParams.ingestion+")"+" to "+cr1.name +" ("+ cr1.energy+") " +" from "+cr2.name+" ("+cr2.energy+") "+ e/ExpParams.secPerStep+ " per sec");
+	cr2.energy_m = cr2.energy_m + e + 0.0000001;
+	cr1.energy_p = cr1.energy_p + e*ExpParams.ingestion;
 	if (cr1.data->lifeparams->hibernated == 1)
 	{
@@ -460,9 +503,8 @@
 	if (cr.data->lifeparams->hibernated == 1)
 	{
-		setForamMeta(cr, cr.data->lifeparams->gen); //unhibernate
+		setForamMeta(cr); //unhibernate
 	} 
 	else
 	{
-		//cr.idleen = (ExpParams.energy_hib * max_chamber_energ[cr.data->lifeparams->gen][Math.min(lastChamberNum(cr), (max_chamber_energ[cr.data->lifeparams->gen].size)-1)])*ExpParams.secPerStep; //hibernate
 		cr.idleen = (ExpParams.energy_hib * cr.energy)*ExpParams.secPerStep; //hibernate
 	}
@@ -470,33 +512,71 @@
 }
 
+
 function onForamsStep(cr)
 {
-	cr.getMechPart(0).orient.set(o);
-	if (visualization(cr))
-	{
-		cr.data->reticulopodiacreature.moveAbs(cr.center_x-ExpParams.zone1_range, cr.center_y-ExpParams.zone1_range, cr.center_z-ExpParams.zone1_range-getProperty(cr.data->lifeparams->gen, "chamber_proculus"));
-	}
-
-	if (deathConditions(cr) == 1)
-	{
-		Populations[0].kill(cr);
-		return;
-	}
-
-	foramMove(cr);
-
-	var repro = foramReproduce(cr);
-	if (repro == 1)
-	{
-		return;
-	}
-
-	cr.data->lifeparams->max_energy_level = Math.max(cr.energy, cr.data->lifeparams->max_energy_level);
-	if  (lastChamberNum(cr) <= chambers[0].size-1)
-	{
-		if ((cr.data->lifeparams->max_energy_level >= max_chamber_energ[cr.data->lifeparams->gen][lastChamberNum(cr)])) 	
-		{
-			foramGrow(cr, lastChamberNum(cr)); 
-		}	
+	//checking for gametogenesis process
+	if (cr.data->lifeparams->division_time > 0)
+	{
+		cr.data->lifeparams->division_time = Math.max(cr.data->lifeparams->division_time-1,0);
+	}
+	//checking for end of gametogenesis
+	else if (cr.data->lifeparams->division_time == 0)
+	{
+		//waiting for gamets fusion
+	}
+	//checking for chamber growth process
+	else if (cr.data->lifeparams->chamber_growth > 0)
+	{
+		cr.data->lifeparams->chamber_growth = Math.max(cr.data->lifeparams->chamber_growth-1,0);
+		//Simulator.print("chamber growing, time left = " + cr.data->lifeparams->chamber_growth*ExpParams.secPerStep);
+		cr.energy_m += ExpParams.chamberCostPerSec * cr.energy * ExpParams.secPerStep;
+
+	 	//Simulator.print("energy " + cr2.energy + " subtracting " + growth_cost);
+	}
+	//checking for end of chamber growth process
+	else if (cr.data->lifeparams->chamber_growth == 0)
+	{
+		foramGrow(cr, lastChamberNum(cr)); 
+		cr.data->lifeparams->chamber_growth = -1;
+		//Simulator.print("chamber "+ (lastChamberNum(cr) + 1) +" complete");
+	}
+	else
+	{
+		//update of metabolism rate
+		if (cr.data->lifeparams->hibernated == 0)
+		{
+			setForamMeta(cr);
+		}
+
+		cr.getMechPart(0).orient.set(o);
+		if (visualization(cr))
+		{
+			cr.data->reticulopodiacreature.moveAbs(cr.center_x-getZoneRange(cr,1), cr.center_y-getZoneRange(cr,1), cr.center_z-getZoneRange(cr,1)-getProperty(cr.data->lifeparams->gen, "chamber_proculus"));
+		}
+
+		if (deathConditions(cr) == 1)
+		{
+			Populations[0].kill(cr);
+			return;
+		}
+
+		foramMove(cr);
+
+		var repro = foramReproduce(cr);
+		if (repro == 1)
+		{
+			return;
+		}
+
+		cr.data->lifeparams->max_energy_level = Math.max(cr.energy, cr.data->lifeparams->max_energy_level);
+
+		//cheking conditions of chamber growth process start
+		if  (lastChamberNum(cr) != chambers[0].size-1)
+		{
+			if ((cr.data->lifeparams->max_energy_level >= energyFromVolume(max_chamber_volume[cr.data->lifeparams->gen][lastChamberNum(cr)],0))) 	
+			{
+				cr.data->lifeparams->chamber_growth = int(ExpParams.chamber_growth_time/ExpParams.secPerStep);
+			}	
+		}
 	}	
 }
@@ -504,6 +584,8 @@
 function deathConditions(cr)
 {
-	if ((cr.energy <= getProperty(cr.data->lifeparams->species,"e_death_level")) || (Math.rnd01 < ExpParams.hunted_prob))
+	if ((cr.energy <= getProperty(cr.data->lifeparams->gen,"e_death_level")*cr.data->lifeparams->max_energy_level) || (Math.rnd01 < ExpParams.hunted_prob))
+	{
 		return 1;
+	}
 	else
 		return 0;
@@ -528,7 +610,7 @@
 // -------------------------------- nutrient begin --------------------------------
 
-function createNutrientGenotype(nutrientsize, zone1_range)
-{
-	return "//0\np:sh=3,sx="+nutrientsize+",sy="+nutrientsize+",sz="+nutrientsize+",ry=1.5,vr=0.0,1.0,0.0";
+function createNutrientGenotype(nutrientradius)
+{
+	return "//0\np:sh=3,sx="+nutrientradius+",sy="+nutrientradius+",sz="+nutrientradius+",ry=1.5,vr=0.0,1.0,0.0";
 }
 
@@ -540,5 +622,5 @@
 function addNutrient()
 {
-	var cr = Populations[1].add(createNutrientGenotype(ExpParams.nutrientsize, ExpParams.zone1_range));
+	var cr = Populations[1].add(createNutrientGenotype(ExpParams.nutrientradius));
 
 	cr.name = "Nutrients";
@@ -553,5 +635,5 @@
 	if (ExpParams.visualize == 1)
 	{
-		var nutsize = ExpParams.nutrientsize*10;
+		var nutsize = ExpParams.nutrientradius*10;
 		var nut = Populations[2].add("//0\np:sh=2,sx="+nutsize+",sy="+nutsize+",sz="+nutsize+",ry=1.5,vr=0.0,1.0,0.0");
 		cr.data->reticulopodiacreature = nut;
@@ -571,7 +653,7 @@
 {
 	nutrientenergywaiting = nutrientenergywaiting + 1;
-	if (nutrientenergywaiting > ExpParams.feedrate/ExpParams.secPerStep)
-	{
-		for (var i = 0; i < ExpParams.nutrientPop; i++)
+	if (nutrientenergywaiting*ExpParams.secPerStep >= ExpParams.foodperiod)
+	{
+		for (var i = 0; i < ExpParams.feedrate*(ExpParams.secPerStep/ExpParams.foodperiod); i++)
 		{    
 			addNutrient();
@@ -668,4 +750,5 @@
 
 // -------------------------------- step end --------------------------------
+//TODO default params values in frams instead of microns/seconds
 
 @include "standard_events.inc"
@@ -692,4 +775,22 @@
 
 prop:
+id:gametSuccessRate
+name:Ratio of successful gamets
+type:f 0.001
+group:Foraminifera
+
+prop:
+id:gametoPeriod
+name:Time of gametogenesis
+type:f 720
+group:Foraminifera
+
+prop:
+id:picoCarbonPerMikro
+name:Picograms of carbon in cubed micrometer
+type:f 0.13
+group:Foraminifera
+
+prop:
 id:secPerStep
 name:Seconds per simulation step
@@ -705,7 +806,25 @@
 
 prop:
+id:divisionCost
+name:Cost of division in pG
+type:f
+group:Foraminifera
+
+prop:
 id:e_repro_cost_diplo
 name:Cost of reproduction
 type:f 0.1 0.9 0.3
+group:Foraminifera
+
+prop:
+id:chamber_growth_time
+name:Time of the chamber growth in seconds
+type:f
+group:Foraminifera
+
+prop:
+id:chamberCostPerSec
+name:Cost of growning chamber per second
+type:f
 group:Foraminifera
 
@@ -791,5 +910,5 @@
 id:e_death_level_haplo
 name:Minimal level of energy to sustain life of haploid
-type:f 0 20 4
+type:f 0 20 0.2
 group:Foraminifera
 
@@ -797,5 +916,5 @@
 id:e_death_level_diplo
 name:Minimal level of energy to sustain life of diploid
-type:f 0 20 0.25
+type:f 0 20 0.2
 group:Foraminifera
 
@@ -857,26 +976,32 @@
 
 prop:
+id:ingestion
+name:Ingestion rate
+type:f
+group:Energy
+
+prop:
 id:energy_nut
 name:Nutrient energy
-type:f 0 1000
+type:f 0 100000
 group:Energy
 
 prop:
 id:feedtrans
-name:Ingestion multiplier
-type:f 0 100
+name:Energy transfer per second
+type:f 0 100000
 group:Energy
 
 prop:
-id:nutrientsize
+id:foodperiod
+name:Time between food occurrences
+type:f 0 100000
+group:Energy
+
+prop:
+id:nutrientradius
 name:Nutrient size
 type:f 0.1 0.9 0.1
 group:Energy
-
-prop:
-id:nutrientPop
-name:Nutrient population size
-group:Energy
-type:d 1 1000 10
 
 prop:
Index: /experiments/frams/foraminifera/data/scripts/foraminifera.inc
===================================================================
--- /experiments/frams/foraminifera/data/scripts/foraminifera.inc	(revision 478)
+++ /experiments/frams/foraminifera/data/scripts/foraminifera.inc	(revision 479)
@@ -55,10 +55,10 @@
 	if (mode->opt == 0) //initial
 	{
-		mode->cr.data->genes = {"min_repro_energies" : [max_chamber_energ[0][getProperty(0, "min_repro_energ")], max_chamber_energ[1][getProperty(1, "min_repro_energ")]], "hibernation" : mode->species};
-		mode->cr.data->lifeparams = {"max_energy_level" : getProperty(0,"energies0"), "gen" : 0,  "hibernated" : 0, "species" : mode->species, "reproduce" : 0, "dir" : randomDir()};
+		mode->cr.data->genes = {"min_repro_energies" : [energyFromVolume(max_chamber_volume[0][getProperty(0, "min_repro_energ")],0), energyFromVolume(max_chamber_volume[1][getProperty(1, "min_repro_energ")],0)], "hibernation" : mode->species};
+		mode->cr.data->lifeparams = {"max_energy_level" : getProperty(0,"energies0"), "gen" : 0,  "hibernated" : 0, "species" : mode->species, "reproduce" : 0, "dir" : randomDir(), "chamber_growth" : -1, "division_time" : -1};
 	}
 	else if (mode->opt  == 1) //child
 	{
-		mode->cr.data->lifeparams = {"max_energy_level" : getProperty(1 - mode->parent_lifeparams->gen,"energies0"), "gen" : 1 - mode->parent_lifeparams->gen,  "hibernated" : 0, "species" : mode->parent_lifeparams->species, "reproduce" : 0, "dir" : randomDir()};
+		mode->cr.data->lifeparams = {"max_energy_level" : getProperty(1 - mode->parent_lifeparams->gen,"energies0"), "gen" : 1 - mode->parent_lifeparams->gen,  "hibernated" : 0, "species" : mode->parent_lifeparams->species, "reproduce" : 0, "dir" : randomDir(), "chamber_growth" : -1, "division_time" : -1};
 		mode->cr.data->genes = mode->parent_genes;
 	}
@@ -68,4 +68,17 @@
 		mode->cr.data->lifeparams = mode->parent_lifeparams;
 	}
+}
+
+function gametsDivision(parent_energy, energy0)
+{
+	var number = 1;
+	var result = parent_energy;
+	while ((result-ExpParams.divisionCost) >= energy0)
+	{
+		result = (result-ExpParams.divisionCost)/2;
+		number *= 2; 
+	}
+	//Simulator.print("parent: " + parent_energy + " result: " + result + " number " + number);
+	return {"energy" : result, "number" : number};
 }
 
@@ -75,6 +88,7 @@
 	if (clone == 1)
 	{
-		energy0 = getProperty(0,"energies0");
-		number = (( 1 - getProperty(1, "e_repro_cost")) * parent.energy) / energy0;
+		var offspring = gametsDivision(parent.energy,getProperty(0,"energies0"));
+		energy0 = offspring->energy;
+		number = offspring->number;
 		new_genes = parent.data->genes;
 		parent.data->lifeparams->gen = 1 - parent.data->lifeparams->gen; //because of reversal of "gen" in createOffspring function
@@ -83,6 +97,8 @@
 	else
 	{
-		energy0 = getProperty(1,"energies0"); 
-		number = (((1 - getProperty(parent.data->lifeparams->gen, "e_repro_cost")) * parent.energy) + ((1 -(getProperty(parent.data->lifeparams->gen, "e_repro_cost"))) * parent2.energy)) / energy0;
+		var offspring1 = gametsDivision(parent.energy,getProperty(1,"energies0"));
+		var offspring2 = gametsDivision(parent2.energy,getProperty(1,"energies0"));
+		energy0 = (offspring1->energy+offspring2->energy); 
+		number = ExpParams.gametSuccessRate*(offspring1->number+offspring2->number)/2;
 		new_genes = [parent.data->genes, parent2.data->genes];
 		gen = 1 - parent.data->lifeparams->gen;
@@ -152,11 +168,27 @@
 				if (parent1 != null && parent2 != null)
 				{
-					reproduce_haploid(parent1, parent2, 0);
-					print_repro_info(parent1);
-					print_repro_info(parent2);
-					pop.kill(parent1);
-					pop.kill(parent2);
-					parent1 = null;
-					parent2 = null;
+					//when parents are ready for reproduction start gametogenesis
+					if (parent1.data->lifeparams->division_time == -1 && parent2.data->lifeparams->division_time == -1)
+					{
+						var time = int(ExpParams.gametoPeriod/ExpParams.secPerStep);
+						parent1.data->lifeparams->division_time = time;
+						parent2.data->lifeparams->division_time = time;
+						parent1.idleen = 0;
+						parent2.idleen = 0;
+						//Simulator.print("parents "+parent1.uid + " " + parent2.uid + " ready to repro: "+Simulator.stepNumber);
+					}
+					//when gametogenesis is finished fuse gamets 
+					else if (parent1.data->lifeparams->division_time == 0 && parent2.data->lifeparams->division_time == 0)
+					{
+						reproduce_haploid(parent1, parent2, 0);
+						print_repro_info(parent1);
+						print_repro_info(parent2);
+						//Simulator.print("parents "+parent1.uid + " " + parent2.uid + " reproduced: "+Simulator.stepNumber);
+						pop.kill(parent1);
+						pop.kill(parent2);
+						parent1 = null;
+						parent2 = null;
+						
+					}
 				}	
    			} 
