#!/bin/awk -f

BEGIN{
	a=320;
	b=2.5;
	c=4.54
	n[0]=c-1;
	p=3984.45092096;
	srand();

	deadcount=bound(1,ARGV[2],ARGV[3]);
	power=bound(0,3500,ARGV[1]);
	assemblies=bound(32,256,ARGV[2]);
	trials=1;
	trial=0;
	max_generation=0
	neutron_count[1]=0;
	neutron_history[1]=0;
	
	if(power == 0)
	{
		print "Power can not be zero!";
		exit;
	};
	
	for(x=1;x<=a;x++)
	{
		n[x]=s=(1.5*(x/a)^(-3/2))*b^(-0.08/(x/a)^2)+c-(2*exp(-x/(1.8745*a))-1)^-1;
	};
	while(trial<trials)
	{
		temp=decay(a,p,assemblies,deadcount,power);
		print trial"\t"temp > "lifetime";
		max_generation=max_generation+temp;
		trial++;
	};
	print "Average lifetime: "max_generation/trials;
	for(nh=1;nh<=length(neutron_history);nh++)
	{
			printf "%s\t%s\n",nh,neutron_count[nh]/neutron_history[nh] > "neutron_profile";
	};
}

function decay(a,p,assemblies,deadcount,power)
{
	dead=0;
	generation=0;
	count=0;
	dcount=0;
	for(r=1;r<=assemblies;r++)
	{
		reactor[r]=0";"n[0];
	};
	while(1)
	{
		partial=0;
		while(dead<deadcount)
		{
			count++;
			drate=0;
			neutron_potential=count_neutron(reactor,assemblies);
			maximum_power=int(p-p/(1+(8.575*neutron_potential/1000)^(9/5)));
			actual_power=minimum(power,maximum_power);
			decay_rate=bound_rate(exp((log(actual_power)-log(p-actual_power)+1.8*log(1000))/1.8)/neutron_potential);
			
			#print "Recording neutron count: "neutron_potential" on count: "count;
			neutron_count[count]=neutron_count[count]+neutron_potential;
			neutron_history[count]++;
			
			if(decay_rate <= 1)
			{
				if(dcount == 0)
				{
					dcount=round(count+(60/decay_rate));
					#print "count: "count" modified-count: "dcount" decay rate: "decay_rate;
				}
			}
			else
			{
				dcount=0;
			};
			if(dcount == 0 || count >= dcount)
			{
				dcount=0;
				#print "Count: "count" Power setting: "power" Maximum power: "maximum_power" neutron_potential :"neutron_potential" Actual power: "actual_power" Decay rate: "decay_rate
				partial=int(0.1*(length(reactor)-dead));
				check_partial=0;
				while(drate<decay_rate)
				{
					drate++;
					while(check_partial<partial)
					{
						check_partial++;
						r=1+int(assemblies*rand());
						split(reactor[r],fuelinfo,";");
						if(fuelinfo[1]<a)
						{
							fuelinfo[1]=fuelinfo[1]+1;
							fuelinfo[2]=n[fuelinfo[1]];
							reactor[r]=fuelinfo[1]";"fuelinfo[2];
							#print "Aged slot: "r"; Age: "fuelinfo[1];
						}
						else
						{
							if (fuelinfo[3]=="")
							{
								reactor[r]=fuelinfo[1]";"fuelinfo[2]";""dead";
								dead++;
								print "Marking assembly no: "r" as dead.  Dead count: "dead". Partial: "partial;
								partial=int(0.1*(length(reactor)-dead));
							};
						};
					};
				};
			};
		};
		if(dead == deadcount)
		{
			generation = count;
			break;
		};
		if(partial == 0)
		{
			generation = count;
			break;
		};
	};
	return generation;
};
function count_neutron(source,assemblies)
{
	nc=0;
	for(s=1;s<=assemblies;s++)
	{		
		split(source[s],fuelinfo,";");
		nc=nc+fuelinfo[2]
	};
	return nc;
};
function minimum(a,b)
{
	if(a < b)
		return a;
	if(a == b)
		return a;
	return b;
};
function bound(a,b,c)
{
	if(b > a)
	{
		if(c >= b)
			return b;
		if(c <= a)
			return a;
		return c;
	}
	if(a > b)
	{
		if(c >= a)
			return a;
		if(c <= b)
			return b;
		return c;
	}
	if (a == b)
	{
		print "Power bound error!"
		exit;
	}
};
function bound_rate(a)
{
	if(a > 8.575)
		return 8.575;
	if(a < 0.05)
		return 0.05;
	return a;
};
function round(x)
{
   ival = int(x)    # integer part, int() truncates

   # see if fractional part
   if (ival == x)   # no fraction
      return ival   # ensure no decimals

   if (x < 0) {
      aval = -x     # absolute value
      ival = int(aval)
      fraction = aval - ival
      if (fraction >= .5)
         return int(x) - 1   # -2.5 --> -3
      else
         return int(x)       # -2.3 --> -2
   } else {
      fraction = x - ival
      if (fraction >= .5)
         return ival + 1
      else
         return ival
   }
};