Program Craig;	

{Craig Engine Simulation Program }
{This simulation employes a Craig Engine approach to simulate the active sampling of volatile}
{organic compounds onto a polymeric adsorbent trap.}

Uses DOS,Crt;

Type
	statarray = array[1..2,1..30] of real;			{Define stationary phase array}
	statptr = ^statarray;
Type
	mobarray = array[1..2,1..30] of real;			{Define mobile phase array}
	mobptr = ^mobarray;

Var 	
	stat							:statptr;
	mob 							:mobptr;
	dist,mobfrac,statfrac,sampcon,
	samvol,	t_ret,phasecap,mobtot				:array[1..4] of real;

	col_len,compno,startnum,plates				:integer;

	deadvol,statsite,molwt,elevol,
	statsat,break_vol,BT,mass_trap				:real;

	plotname						:string[12];
	filename						:string[6];
	n,times,stage						:word;
	key							:char;

Const ava = 6e23;										{Avagadro's No.}

{------------------------------------------------------------------------------------------------------------------------}
Procedure Initialise; 	{set filename for elution data}
Var
name							       : String[6];

Begin
        Writeln('Enter filename, 6 letters');			{no. of components}
        Readln(name);
	filename:=name;
	startnum:=1;
End;
{------------------------------------------------------------------------------------------------------------------------}
Procedure Reset_Array; 	{ sets mobile and stat array values to zero}

Var loop1,loop2											:integer;

Begin
	for loop1:=1 to plates do				
	begin
		for loop2:=1 to 2 do
		begin
			stat^[loop2,loop1]:=0;mob^[loop2,loop1]:=0;
		end;
	end;
End;
{------------------------------------------------------------------------------------------------------------------------}
Procedure Input;

Var	plotfix							:string[2];
        number	                				:real;
	q,r							:word;

Begin
	Writeln('number of components?');			{no. of components}
        Readln(r);
	if ( r > 4) or ( r < 1 ) then r := 1;
	compno := r;

	For  q  := 1 to compno do
	begin													
		writeln('Partition parameter component ',q,' ?');{K component }
			readln(number);
		if ( number < 1) or ( number > 1e9 ) then number := 1;
		dist[q] := number;

		writeln('Concentration g/cm3',q,' ?');		{Concentration component 1}
		readln(number);
		if  number > 1E-3  then number := 1E-9;
		sampcon[q] := number;
        end;

	molwt := 78;						{set molecular weight of analyte}

	str(times,plotfix);
	plotname := filename+ plotfix + '.dat';			{sets up filename}

	BT := sampcon[1]*0.98					{sets breakthough vol for expt}
End;
{------------------------------------------------------------------------------------------------------------------------}
Procedure Process_Input;	{Use data from input procedure to obtain necessary variables for program,i.e.
 capacity of stationary phase element}

Begin
	mass_trap := 0;
	break_vol := 0;
	statsat := molwt*statsite/(ava*plates);			{capacity stat. element w.r.t mass of analyte}
End;
{------------------------------------------------------------------------------------------------------------------------}
Procedure Sample_to_Mob(num : integer); 			{place sample conc. in 1st element of array}

Begin
	mob^[num,1] := sampcon[num];					
End;
{------------------------------------------------------------------------------------------------------------------------}
Procedure Multcomp;	{Performs actual distribution between both phases}

Var wxprime,divisor,accuracy					:real;
	subfrac							:double;
	mobconc,statmass					:array[1..4] of real;
	dist_end						:boolean;
	q							:integer;

Begin
for q:=1 to compno do
begin
	mobconc[q] := mobtot[q];				{all component in mobile, stat empty}
	statmass[q] := 0;
end;
subfrac := -mobtot[1]*0.5;
dist_end := false;
Repeat
	divisor := 0;
	mobconc[1] := mobconc[1] + subfrac;
	statmass[1] := statmass[1] - (subfrac*elevol);
	if compno > 1 then					{distribute other components relative to first}
	begin
		for q := 2 to compno do
			begin
													
				if ((dist[1]*mobconc[1])+(statmass[1]*dist[q]/elevol)) = 0 then
				begin
					statmass[q] := mobtot[q]*elevol 
				end
				else
				begin
					statmass[q] := mobtot[q]*statmass[1]*dist[q]*elevol/
					((dist[1]*mobconc[1]*elevol)+(statmass[1]*dist[q]));
				end;
				mobconc[q] := mobtot[q]-(statmass[q]/elevol);
			end;
	end;
	for q := 1 to compno do					{langmuir eqn. based on current concs}
	begin
		divisor:=divisor+(dist[q]*mobconc[q]);
	end;
	wxprime := statsat*dist[1]*mobconc[1]/(1+divisor);
	subfrac := abs(subfrac);
	if wxprime < statmass[1] then 
	begin
		subfrac := (subfrac*0.5) 							{compares mass at surface cf. prediction}
	end
	else
	begin
		subfrac := -(subfrac*0.5);
	end;
	accuracy := statmass[1]*1e-3;				{accuracy of difference before calc. stops}
	if (abs(wxprime-statmass[1]) < accuracy)
	or (wxprime = 0)
	or (abs((statmass[1]/wxprime)-1) < 0.001) then
	begin
		dist_end := true;
	end;
Until dist_end;

for q := 1 to compno do
begin
	mobfrac[q] := mobconc[q];
	statfrac[q] := statmass[q];
end;
End;
{------------------------------------------------------------------------------------------------------------------------}
Procedure Distribute(var curr:word);	{Langmuir distribution procedure}

Var q													:integer;

Begin
	for q := 1 to compno do
	begin
		mobfrac[q] := mob^[q,curr];
		statfrac[q] := stat^[q,curr];
		mobtot[q] := mobfrac[q]+(statfrac[q]/elevol);
	end;
	Multcomp;						{distribution calculation}
	for q := 1 to compno do
	begin
		mob^[q,curr] := mobfrac[q];
		stat^[q,curr] := statfrac[q];
	end;
End;
{------------------------------------------------------------------------------------------------------------------------}
Procedure Craig_Engine;	{perform Craig Engine routine}

Var  statposn,current,mobcurr1,mobcurr2				:word;

	t													:text;
	q,no_eluent,no_cycles					:integer;
	elu_mass,elu_vol					:real;
	stop							:boolean;

Begin
	elevol := deadvol/plates;
	statposn := 1;
	no_cycles := 0;
	for q := 1 to compno do
	begin
		Sample_to_Mob(q)								{put sample in first element}
	end;
	Writeln('At work');
	Assign(t,'c:\temp\'+plotname);				{opens file to save data}
	Rewrite(t);
	Writeln(t,'Single comp.');
	writeln(t,'K = ',dist[1]);								{saves component characteristics}
	writeln(t,'[c] = ',sampcon[1]);
	no_eluent := 1;										{counts integer volume for saving}
	elu_vol := elevol;
	elu_mass := 0;										{counts actual volume passed} 
	stop := false;

	Repeat
	begin
		for current := 1 to statposn do
		begin
			Distribute(current);							{call distribute routine}
		end;
		if statposn = plates then							{writes effluent conc. every x ml}
		begin
			if elu_vol > no_eluent then
			begin
				writeln(t,mob^[1,statposn],mob^[2,statposn]);	
				writeln(times,mob^[1,statposn]);				
				no_eluent := no_eluent + 1;	
				no_cycles:= no_cycles+1;
			end;
			elu_vol := elu_vol + elevol;
			elu_mass := elu_mass + (sampcon[1] - mob^[1,statposn]);
		end;
		if mob^[1,plates] > BT then
		begin
			mass_trap := elu_mass;
			break_vol := elu_vol;
			stop := true;
		end;
		for q := 1 to compno do				{shift mob. phase along one element}
		begin
			for current := statposn downto 1 do
			begin
				mob^[q,current+1] := mob^[q,current];
			end;
			mob^[q,1] := sampcon[q];
		end;
		if statposn<plates then statposn := statposn+1;
	end;
	Until stop;
	mass_trap := mass_trap * elevol * 1000;
	writeln(t,'Mass trapped = ',mass_trap);
	writeln(t,'BTV = ',break_vol);
	Close(t);
End;
{------------------------------------------------------------------------------------------------------------------------}
Begin	{The main loop of the program}
	Writeln;
	Writeln('Multicomponent Tenax tube simulation');
	Writeln;
	Writeln('Optimised Version 3.0');

	col_len := 4;										{column length - 4 cm}
	deadvol := 0.312;									{dead volume ml}
	statsite := 1.5e18;									{no. adsorption sites in the stat phase}

	plates := 6;											{no. elements in stat. array}

	New(stat); New(mob);					{set up mobile & stat phase arrays}

	Initialise;											{calls analyte parameter input}

	For times := startnum to 1 do
	begin
		Reset_Array; 									{Clear array contents}

		Input;											{Set up Compound Characteristics}

		Process_Input;

		Craig_Engine;									{Main Calculation step}
	end;
End.