#define BASES 3
#define QUADS 30
#define ENGLAND 0
#define SCOTLAND 1
#define WALES 2

continuous BtoQ[BASES][QUADS];
continuous TravelDist;
int aircraft[BASES];
int hi[QUADS], lo[QUADS];
int knots[BASES][QUADS];

data()
{
aircraft = { 41, 40, 48 };

lo = {3,2,4,4,1,1,1,1,1,2,6,3,2,5,2,2,3,4,3,5,4,5,2,1,3,5,3,4,3,4};

hi = {4,6,5,9,8,6,9,5,3,4,9,6,7,8,8,4,6,5,9,8,6,9,5,3,4,9,6,7,8,8};

knots =     { 475, 200, 575,
		   550, 225, 475,
		   650, 250, 400,
		   750, 300, 300,
		   350, 225, 700,
       		   450, 275, 600,
		   550, 275, 500,
		   650, 300, 425,
		   750, 350, 325,
		   850, 400, 275,
		   250, 375, 800,
		   350, 350, 725,
		   450, 325, 475,
		   550, 400, 525,
		   650, 400, 475,
		   750, 400, 400,
		   850, 475, 325,
		   950, 500, 300,
		   300, 425, 825,
		   375, 400, 750,
		   475, 400, 700,
		   575, 400, 600,
		   675, 400, 525,
		   750, 475, 475,
		   850, 500, 425,
		   950, 550, 400,
		   500, 475, 725,
		   600, 475, 650,
		   700, 475, 600,
		   800, 500, 525};
}

assign_all_aircraft()
{
	and(int i=0;i<BASES;i++)					// From each base
		sigma(int j=0;j<QUADS;j++) BtoQ[i][j]
		== aircraft[i];						// Send all aircraft
}

stay_within_limits()
{
	and(int j=0;j<QUADS;j++)
		sigma(int i=0;i<BASES;i++) BtoQ[i][j] >= lo[j];
	and(int j=0;j<QUADS;j++)
		sigma(int i=0;i<BASES;i++) BtoQ[i][j] <= hi[j];

	and(int j=0;j<QUADS;j++)					// Redundant bounding
		and(int i=0;i<BASES;i++)				// constraints to help
			BtoQ[i][j] <= hi[j];				// tighten model
}

set_travel_distance()
{
	TravelDist == 
		sigma(int i=0;i<BASES;i++)
			sigma(int j=0;j<QUADS;j++)
				knots[i][j]*BtoQ[i][j];
}

branch_on_bases(int j)
{
	or(int i=0;i<BASES;i++) {				// Select base i
		and(int t=0;t<BASES;t++)
			if t != i;							// Other bases t cannot send 
			then BtoQ[t][j] == 0;					// planes to quadrant j
		up_then_down(BtoQ[i][j],lo[j],hi[j]);		// Branch on planes
	}
}

determine_quad(int & quad)
{
	and(int j=0;j<QUADS;j++) 
		if not ok(j);
		then {
			quad = j; 			// Record injured quadrant
			return;			// Return with success
		}
	0==1;			// Return with failure
}

ok(int j)
{
	c_or(int i=0;i<BASES;i++){				// Is there a base i sending 
		integral(BtoQ[i][j]);					// an integral number of 
										// planes to quadrant j  
		and(int t=0;t<BASES;t++)				// such that all others
			if t != i;						// send none 
			then wp(BtoQ[t][j]) == 0;  
	} 
}

assign_bases_to_quadrants()
{
int quad;			// Storeback variable

	// A sequence of successful calls to 
	// branch_on_bases(quad0),...,branch_on_bases(quadk-1) 
	// where quad0,...,quadk-1 are distinct quadrants
	and(;;) {			// Loop variable plays no role
		min: TravelDist;
		if determine_quad(quad);			// Injured quadrant ?
		then branch_on_bases(quad);		// Repair injury
		else break; 						// Eureka effect
	}
}

up_then_down(continuous X, int lo,hi)
{
int start;
	start = floor(X);
	either {
		//X==start,...,X==k-1 do not lead to a solution
		or(int k=start+1;k<=hi;k++){// If k>hi, sidetrack
			c_either X == k;
			or break;			// Sidetrack to lower loop
		}
	}
	or {
		// X==start,...,X==k+1 do not lead to a solution
		or(int k=start;k>=lo;k--){	// If k<lo, backtrack
			c_either X == k;
			or break;			// Backtrack
		}
	}
}

output()
{
	printf("Travel distance is %.0f \n", TravelDist);
}

2lp_main()
{
	data();
	stay_within_limits();
	assign_all_aircraft();
	set_travel_distance();
	
	find_min: TravelDist; 
	subject_to
		assign_bases_to_quadrants();
	
	output(); 
}
 