{
Utilities for BCAM Device Calibration and Measurement Transformation
Copyright (C) 2004-2008 Kevan Hashemi, hashemi@brandeis.edu, Brandeis University

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
}

unit bcam;
{
	bcam contains routines that transform between BCAM coordintes, global coordinates,
	and image coordinates. It contains routines that analyze calibration measurements
	and calculate the bcam calibration constants. BCAM measurements we obtain separately
	using the routines in our spot unit. Consequently, this unit needs only the
	routines in utils to function, because it does not deal with images or displays at
	all.
}
interface

uses
	utils;
	
const
	num_roll_cage_pairs=3*2*1; {number of combinations of two from a set of four}
	bcam_z_angle=0.5596/2;{z-axis to slot-cone line in radians}
	bcam_mid_z=-40;{mm approx for discrimination}
	bcam_ccd_center_x=1.720;{along image coordinate x-axis (mm)}
	bcam_ccd_center_y=1.220;{along image coordinate y-axis (mm)}
	bcam_front_source_z=0.36;{front laser facet bcam-coord z-position (mm)}
	bcam_rear_source_z=-82.85;{rear laser facet bcam-coord z-position (mm)}
	num_mounts_per_pair=2;
	num_mounts_per_quad=4;
	num_sources_per_range=4;
	num_sources_per_pair=2;
	num_ranges_per_mount=2;
	bcam_camera_pivot_xy_spread=0.080;{mm}
	bcam_camera_pivot_z_spread=4.0;{mm}
	bcam_camera_axis_xy_spread=0.0001;{rad}
	bcam_camera_ccd_to_pivot_spread=0.300;{mm}
	bcam_camera_ccd_rotation_spread=0.001;{rad}
	bcam_source_spread=0.040;{mm}

type
{
	bcam_camera_type gives the pivot, ccd_to_pivot distance, axis, and
	ccd_rotation of the camera in 'bcam coordinates'. The bcam
	coordinate z-axis runs along the nominal optical axis of the
	camera. It points in the direction of the arrow made by the
	mounting balls. The x-axis is in the plane of the mounting balls.
	It points to the left when you are looking in the z-direction with
	the balls below. The y-axis completes a right-handed system.
}
	bcam_camera_type=record
		pivot:xyz_point_type;{bcam coordinates of pivot point (mm)}
		axis:xyz_point_type;{cosine vectors of camera axis}
		ccd_to_pivot:real;{from ccd to pivot along camera axis (mm)}
		ccd_rotation:real;{rotation of ccd about camera axis (rad)}
		id:short_string;{identifier}
	end;
{
	bcam_source_pair_type gives the centers of two light sources in bcam coordinates.
}
	bcam_source_pair_type=record
		sources:array [1..num_sources_per_pair] of xyz_point_type;
			{bcam coordinates of pivot point (mm)}
		id:short_string;
			{identifier}
	end;
{
	source_locations_type gives the coordinates of the sources in x and y, where x and y
	are local coordinates in the source plate, parallel to the bcam local x and y axese.
}
	source_locations_type=array [1..num_sources_per_range] of xy_point_type;
{
	calib_datum_type gives the position of a light spot on the camera ccd in ccd 
	coordinates, and the z-coordinate of the light source in global coordinates. The routine
	is for calibration stands in which the global z-coordinate is parallel to the bcam
	z-coordinate.
}
	calib_datum_type=record
		spot_center:xy_point_type;{light spot center image coordinates (um)}
		source_range:real;{z-coordinate of source in global coordinates (mm)}
	end;
{
	bcam_camera_calib_input_type contains all the information needed to determine the
	calibration parameters of a bcam camera.
}
	bcam_camera_calib_input_type=record
		mounts:array [1..num_mounts_per_quad] of kinematic_mount_type;
		measurements:array [1..num_mounts_per_quad,1..num_ranges_per_mount,1..num_sources_per_range] 
			of calib_datum_type;
		source_locations:source_locations_type;
		id,time:short_string;
		axis_direction:real;
	end;
{
	bcam_camera_calib_output_type contains the six bcam_camera_type records produced
	by the six pairs of orientations available from a roll-cage calibration.
}
	bcam_camera_calib_output_type=record
		pairs:array [1..num_roll_cage_pairs] of bcam_camera_type;
		average:bcam_camera_type;
		spread:bcam_camera_type;
		id,time:short_string;
	end;
{
	bcam_source_pair_calib_input_type contains all the information needed to determine the
	positions of a pair of bcam sources.
}
	bcam_source_pair_calib_input_type=record
		mounts:array [1..num_mounts_per_quad] of kinematic_mount_type;
		viewing_x_direction:real;
		viewing_scale:real;
		measurements:array [1..num_mounts_per_quad,1..num_sources_per_pair] 
			of calib_datum_type;
		id,time:short_string;
	end;
{
	bcam_source_pair_calib_output_type contains the six bcam_source_pair_type
	records produced by the six pairs of orientations available from a roll-cage calibration.
}
	bcam_source_pair_calib_output_type=record
		pairs:array [1..num_roll_cage_pairs] of bcam_source_pair_type;
		average:bcam_source_pair_type;
		spread:bcam_source_pair_type;
		id,time:short_string;
	end;
	
procedure bcam_check_camera_spread(var f:string;spread:bcam_camera_type);
procedure bcam_check_source_pair_spread(var f:string;spread:bcam_source_pair_type);
function bcam_coordinates_from_camera_mount(mount:kinematic_mount_type):coordinates_type;
function bcam_from_image_point(p:xy_point_type;camera:bcam_camera_type):xyz_point_type;
function bcam_from_global_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;
function bcam_from_global_point(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;
function bcam_from_global_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type;
function bcam_from_global_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type;
function image_from_bcam_point(p:xyz_point_type;camera:bcam_camera_type):xy_point_type;
function global_from_bcam_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;
function global_from_bcam_point(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;
function global_from_bcam_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type;
function global_from_bcam_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type;
function nominal_bcam_camera(axis_direction:real):bcam_camera_type;
function bcam_source_bearing(spot_center:xy_point_type;
	camera:bcam_camera_type):xyz_line_type;
function bcam_source_position(spot_center:xy_point_type;range:real;
	camera:bcam_camera_type):xyz_point_type;
function read_bcam_camera(var f:string):bcam_camera_type;
function bcam_camera_from_string(s:short_string):bcam_camera_type;
procedure write_bcam_camera_calib(var f:string;
	calib_output:bcam_camera_calib_output_type;
	verbose,check_spread:boolean);
procedure write_bcam_camera(var f:string;camera:bcam_camera_type);
procedure write_bcam_source_pair_calib(var f:string;
	calib_output:bcam_source_pair_calib_output_type;
	verbose,check_spread:boolean);
function bcam_camera_calib(calib_input:bcam_camera_calib_input_type):
	bcam_camera_calib_output_type;
function bcam_source_pair_calib(calib_input:bcam_source_pair_calib_input_type):
	bcam_source_pair_calib_output_type;

implementation

const
	n=3;{three-dimensional space}

type
{
	pair_data_type contains all the information the bcam_pair_calib
	routine needs to calibrate a bcam camera.
}
	pair_data_type=record
		mounts:array [1..num_mounts_per_pair] of kinematic_mount_type;
		measurements:
			array [1..num_mounts_per_pair,1..num_ranges_per_mount,1..num_sources_per_range] 
			of calib_datum_type;
		source_locations:source_locations_type;
		id:short_string;{identifier}
		axis_direction:real;
	end;

{
	camera_limit returns a bcam_camera_type with fields set equal to the 
	largest acceptable spreads for bcam camera paramters.
}
function camera_limit:bcam_camera_type;

var
	limit:bcam_camera_type;
	
begin
	with limit do begin
		pivot.x:=bcam_camera_pivot_xy_spread;
		pivot.y:=bcam_camera_pivot_xy_spread;
		pivot.z:=bcam_camera_pivot_z_spread;
		axis.x:=bcam_camera_axis_xy_spread;
		axis.y:=bcam_camera_axis_xy_spread;
		axis.z:=0;
		ccd_to_pivot:=bcam_camera_ccd_to_pivot_spread;
		ccd_rotation:=bcam_camera_ccd_rotation_spread;
		id:='limit';
	end;
	camera_limit:=limit;
end;

{
	source_pair_limit returns a source_pair_type with fields set equal to the 
	largest acceptable spreads for bcam source_pair paramters.
}
function source_pair_limit:bcam_source_pair_type;

var
	limit:bcam_source_pair_type;
	source_num:integer;
	
begin
	with limit do begin
		for source_num:=1 to num_sources_per_pair do begin
			with sources[source_num] do begin
				x:=bcam_source_spread;
				y:=bcam_source_spread;
				z:=bcam_source_spread;
			end;
		end;
		id:='limit';
	end;
	source_pair_limit:=limit;
end;

{
	bcam_check_source_pair_spread examines the fields of a camera_type that represents
	the spread of values obtained from the six pairs of positions available from a
	camera calibration. The routine appends spread excesses to a string.
}
procedure bcam_check_source_pair_spread(var f:string;spread:bcam_source_pair_type);

const
	fsr=6;
	fsd=3;
	fsi=1;

var
	limit:bcam_source_pair_type;
	source_num:integer;
	
begin
	limit:=source_pair_limit;
	for source_num:=1 to num_sources_per_pair do begin
		with spread.sources[source_num] do begin
			if x>limit.sources[source_num].x then 
				writestr(f,f,eol,'WARNING: source ',source_num:fsi,' x spread exceeds ',
					limit.sources[source_num].x:fsr:fsd,' mm');
			if y>limit.sources[source_num].y then 
				writestr(f,f,eol,'WARNING: source ',source_num:fsi,' y spread exceeds ',
					limit.sources[source_num].y:fsr:fsd,' mm');
			if z>limit.sources[source_num].z then 
				writestr(f,f,eol,'WARNING: source ',source_num:fsi,' z spread exceeds ',
					limit.sources[source_num].z:fsr:fsd,' mm');
		end;
	end;
end;

{
	bcam_check_camera_spread examines the fields of a camera_type that represents
	the spread of values obtained from the six pairs of positions available from a
	camera calibration. The routine appends spread excesses to a string.
}
procedure bcam_check_camera_spread(var f:string;spread:bcam_camera_type);

const
	fsr=6;
	fsd=3;

var
	limit:bcam_camera_type;
	
begin
	limit:=camera_limit;
	with limit do begin
		if spread.pivot.x>pivot.x then 
			writestr(f,f,eol,'WARNING: pivot.x spread exceeds ',
				pivot.x:fsr:fsd,' mm.');
		if spread.pivot.y>pivot.y then 
			writestr(f,f,eol,'WARNING: pivot.y spread exceeds ',
				pivot.y:fsr:fsd,' mm.');
		if spread.pivot.z>pivot.z then 
			writestr(f,f,eol,'WARNING: pivot.z spread exceeds ',
				pivot.z:fsr:fsd,' mm.');
		if spread.axis.x>axis.x then 
			writestr(f,f,eol,'WARNING: axis.x spread exceeds ',
				mrad_per_rad*axis.x:fsr:fsd,' mrad.');
		if spread.axis.y>axis.y then 
			writestr(f,f,eol,'WARNING: axis.y spread exceeds ',
				mrad_per_rad*axis.y:fsr:fsd,' mrad.');
		if spread.axis.z>axis.z then 
			writestr(f,f,eol,'WARNING: axis.z spread exceeds ',
				axis.z:fsr:fsd,' mm.');
		if spread.ccd_to_pivot>ccd_to_pivot then 
			writestr(f,f,eol,'WARNING: ccd_to_pivot spread exceeds ',
				ccd_to_pivot,' mm.');
		if spread.ccd_rotation>ccd_rotation then 
			writestr(f,f,eol,'WARNING: ccd_rotation spread exceeds ',
				mrad_per_rad*ccd_rotation:fsr:fsd,' mrad.');
	end;
end;

{
	read_bcam_camera reads a camera type from a string.
}
function read_bcam_camera(var f:string):bcam_camera_type;

var 
	i:integer;
	camera:bcam_camera_type;

begin
	with camera do begin
		id:=read_word(f);
		pivot:=read_xyz(f);
		axis:=read_xyz(f);
		with axis do begin
			x:=x/mrad_per_rad;
			y:=y/mrad_per_rad;
		end;
		ccd_to_pivot:=read_real(f);
		ccd_rotation:=read_real(f);
		ccd_rotation:=ccd_rotation/mrad_per_rad;
	end;
	read_bcam_camera:=camera;
end;

{
	bcam_camera_from_string converts a string into a bcam_camera_type;
}
function bcam_camera_from_string(s:short_string):bcam_camera_type;
begin
	bcam_camera_from_string:=read_bcam_camera(s);
end;

{
	write_bcam_camera appends a camera type to a string.
}
procedure write_bcam_camera(var f:string;camera:bcam_camera_type);

const 
	fsr=10;fsd=3;fsdr=3;fsdd=4;	

begin
	with camera do begin
		writestr(f,f,'id: ',id,eol);
		with pivot do 
			writestr(f,f,'pivot (mm):  ',x:fsr:fsd,y:fsr:fsd,z:fsr:fsd,eol);
		with axis do 
			writestr(f,f,'axis (mrad):   ',
				x*mrad_per_rad:fsr:fsdr,y*mrad_per_rad:fsr:fsdr,eol);
		writestr(f,f,'axis z-direction: ',axis.z:fsr:fsdr,eol);
		writestr(f,f,'ccd_to_pivot (mm):   ',ccd_to_pivot:fsr:fsd,eol);
		writestr(f,f,'ccd_rotation (mrad):   ',ccd_rotation*mrad_per_rad:fsr:fsdr,eol);
	end;
end;

{
	string_from_bcam_source_pair appends a camera type to a string, using only one line.
}
function string_from_bcam_source_pair(source_pair:bcam_source_pair_type)
	:short_string;

const 
	fsr=8;fsd=3;fsdr=3;fsid=10;fss=4;

var 
	source_num:integer;
	f:short_string='';

begin
	with source_pair do begin
		writestr(f,f,id:fsid);
		for source_num:=1 to num_sources_per_pair do
			with sources[source_num] do
				writestr(f,f,x:fsr:fsd,y:fsr:fsd);
		writestr(f,f,sources[1].z:fsr:fsd);
	end;
	string_from_bcam_source_pair:=f;
end;

{
	string_from_bcam_camera appends a camera type to a string, using only one line.
}
function string_from_bcam_camera(camera:bcam_camera_type):short_string;
	
const 
	fsr=8;fsd=3;fsdr=3;fsid=10;fss=4;

var 
	f:short_string='';

begin
	with camera do begin
		writestr(f,f,id:fsid);
		with pivot do 
			writestr(f,f,x:fsr:fsd,y:fsr:fsd,z:fsr:fsd);
		with axis do 
			writestr(f,f,x*mrad_per_rad:fsr:fsdr,y*mrad_per_rad:fsr:fsdr,z:fss:0);
		writestr(f,f,ccd_to_pivot:fsr:fsd);
		writestr(f,f,' ',ccd_rotation*mrad_per_rad:fsr:fsdr);
	end;
	string_from_bcam_camera:=f;
end;

{
	write_bcam_source_pair_calib appends all six source pair calibrations to a string 
	and the average and spreads of the parameters as well.
}
procedure write_bcam_source_pair_calib(var f:string;
	calib_output:bcam_source_pair_calib_output_type;
	verbose,check_spread:boolean);

var 
	source_num,pair_num:integer;
	s,w:short_string='';
	
begin
	if verbose then with calib_output do begin
		writestr(f,f,eol);
		if average.sources[1].z>bcam_mid_z then 
			writestr(f,f,'Calibration Constants for Front Sources on Device ',id,':',eol)
		else
			writestr(f,f,'Calibration Constants for Rear Sources on Device ',id,':',eol);
		writestr(f,f,' -------------------------------------------------',eol);
		writestr(f,f,'          | Inner Source  | Outer Source  |         ',eol);
		writestr(f,f,'          |-----------------------|-------|         ',eol);
		writestr(f,f,'          |   x   |   y   |   x   |   y   |   z     ',eol);
		writestr(f,f,'    Pair  |  (mm) |  (mm) |  (mm) |  (mm) |  (mm)   ',eol);
		writestr(f,f,' -------------------------------------------------',eol);
		for pair_num:=1 to num_roll_cage_pairs do
			f:=f+string_from_bcam_source_pair(pairs[pair_num])+eol;
		writestr(f,f,' -------------------------------------------------',eol);
		f:=f+string_from_bcam_source_pair(average)+eol;
		f:=f+string_from_bcam_source_pair(spread)+eol;
		f:=f+string_from_bcam_source_pair(source_pair_limit)+eol;
		writestr(f,f,' -------------------------------------------------',eol);
		writestr(f,f,'Calibration performed at time ',time,eol);
	end else begin
		s:=string_from_bcam_source_pair(calib_output.average);
		w:=read_word(s);
		f:=f+s;
	end;
	
	if check_spread then bcam_check_source_pair_spread(f,calib_output.spread);
end;

{
	If verbose is set, then write_bcam_camera_calib appends all six camera calibrations 
	to a string and the average and spreads of the parameters as well. If verbose is not
	set, then write_bcam_camera_calib appends the average calibration to the string. If
	check_range, then the routine adds warnings to the string if the calibrtion parameter
	ranges exceed the maxima defined above.
}
procedure write_bcam_camera_calib(var f:string;
	calib_output:bcam_camera_calib_output_type;
	verbose,check_spread:boolean);
	
var 
	pair_num:integer;
	s,w:short_string='';
	
begin
	if verbose then with calib_output do begin
		writestr(f,f,eol);
		if (average.axis.z>0) then 
			writestr(f,f,'Calibration Constants for Front-Facing Camera on Device ',id,':',eol)
		else
			writestr(f,f,'Calibration Constants for Rear-Facing Camera on Device ',id,':',eol);
		writestr(f,f,' -----------------------------------------------------------------------',eol);
		writestr(f,f,'          |      Pivot Position   |   Axis Direction  |  CCD  |  CCD    ',eol);
		writestr(f,f,'          |-----------------------|-------------------| -to-  |  Rot-   ',eol);
		writestr(f,f,'          |    x  |   y   |   z   |    x  |    y  |   | Pivot |  ation  ',eol);
		writestr(f,f,'    Pair  |  (mm) |  (mm) |  (mm) | (mrad)| (mrad)| z |  (mm) |  (mrad) ',eol);
		writestr(f,f,' -----------------------------------------------------------------------',eol);
		for pair_num:=1 to num_roll_cage_pairs do
			f:=f+string_from_bcam_camera(pairs[pair_num])+eol;
		writestr(f,f,' -----------------------------------------------------------------------',eol);
		f:=f+string_from_bcam_camera(average)+eol;
		f:=f+string_from_bcam_camera(spread)+eol;
		f:=f+string_from_bcam_camera(camera_limit)+eol;
		writestr(f,f,' -----------------------------------------------------------------------',eol);
		writestr(f,f,'Calibration performed at time ',time,eol);
	end else begin
		s:=string_from_bcam_camera(calib_output.average);
		w:=read_word(s);
		f:=f+s;
	end;
	
	if check_spread then bcam_check_camera_spread(f,calib_output.spread);
end;

{
	bcam_camera_average calculates the average calibration constants from a sequence
	of bcam_camera_type records.
}
function bcam_camera_average(cp:pointer;num_calibs:integer):bcam_camera_type;

var
	calib_num:integer;
	calib_ptr:^bcam_camera_type;
	sum_calib:bcam_camera_type;
	
begin
	with sum_calib do begin
		pivot:=xyz_origin;
		axis:=xyz_origin;
		ccd_to_pivot:=0;
		ccd_rotation:=0;
	end;
	for calib_num:=1 to num_calibs do begin
		calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_camera_type));
		with sum_calib do begin
			pivot.x:=pivot.x+calib_ptr^.pivot.x/num_calibs;
			pivot.y:=pivot.y+calib_ptr^.pivot.y/num_calibs;
			pivot.z:=pivot.z+calib_ptr^.pivot.z/num_calibs;
			axis.x:=axis.x+calib_ptr^.axis.x/num_calibs;
			axis.y:=axis.y+calib_ptr^.axis.y/num_calibs;
			axis.z:=axis.z+calib_ptr^.axis.z/num_calibs;
			ccd_to_pivot:=ccd_to_pivot+calib_ptr^.ccd_to_pivot/num_calibs;
			ccd_rotation:=ccd_rotation+calib_ptr^.ccd_rotation/num_calibs;
		end;
	end;
	sum_calib.id:='average';
	bcam_camera_average:=sum_calib;
end;

{
	bcam_camera_spread calculates the spread of each calibration constant from a sequence
	of bcam_camera_type records.
}
function bcam_camera_spread(cp:pointer;num_calibs:integer):bcam_camera_type;

var
	calib_num:integer;
	calib_ptr:^bcam_camera_type;
	max,min,spread:bcam_camera_type;
	
begin
	for calib_num:=1 to num_calibs do begin
		calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_camera_type));
		if calib_num=1 then begin
			max:=calib_ptr^;
			min:=calib_ptr^;
		end
		else begin
			with calib_ptr^ do begin
				if max.pivot.x<pivot.x then max.pivot.x:=pivot.x;
				if max.pivot.y<pivot.y then max.pivot.y:=pivot.y;
				if max.pivot.z<pivot.z then max.pivot.z:=pivot.z;
				if max.axis.x<axis.x then max.axis.x:=axis.x;
				if max.axis.y<axis.y then max.axis.y:=axis.y;
				if max.axis.z<axis.z then max.axis.z:=axis.z;
				if max.ccd_to_pivot<ccd_to_pivot then max.ccd_to_pivot:=ccd_to_pivot;
				if max.ccd_rotation<ccd_rotation then max.ccd_rotation:=ccd_rotation;
				if min.pivot.x>pivot.x then min.pivot.x:=pivot.x;
				if min.pivot.y>pivot.y then min.pivot.y:=pivot.y;
				if min.pivot.z>pivot.z then min.pivot.z:=pivot.z;
				if min.axis.x>axis.x then min.axis.x:=axis.x;
				if min.axis.y>axis.y then min.axis.y:=axis.y;
				if min.axis.z>axis.z then min.axis.z:=axis.z;
				if min.ccd_to_pivot>ccd_to_pivot then min.ccd_to_pivot:=ccd_to_pivot;
				if min.ccd_rotation>ccd_rotation then min.ccd_rotation:=ccd_rotation;
			end;
		end;
	end;
	with spread do begin
		pivot:=xyz_difference(max.pivot,min.pivot);
		axis:=xyz_difference(max.axis,min.axis);
		ccd_to_pivot:=max.ccd_to_pivot-min.ccd_to_pivot;
		ccd_rotation:=max.ccd_rotation-min.ccd_rotation;
	end;
	spread.id:='spread';
	bcam_camera_spread:=spread;
end;

{
	bcam_source_pair_average calculates the average calibration constants from a sequence
	of bacm_source_pair_type records.
}
function bcam_source_pair_average(cp:pointer;num_calibs:integer):bcam_source_pair_type;

var
	calib_num,source_num:integer;
	calib_ptr:^bcam_source_pair_type;
	sum_calib:bcam_source_pair_type;
	
begin
	for source_num:=1 to num_sources_per_pair do begin
		sum_calib.sources[source_num]:=xyz_origin;
	end;
	for calib_num:=1 to num_calibs do begin
		calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_source_pair_type));
		for source_num:=1 to num_sources_per_pair do begin
			sum_calib.sources[source_num]:=
				xyz_sum(sum_calib.sources[source_num],
					xyz_scale(calib_ptr^.sources[source_num],1/num_calibs));
		end;
	end;
	sum_calib.id:='average';
	bcam_source_pair_average:=sum_calib;
end;

{
	bcam_source_pair_spread calculates the spread of each calibration constant from a sequence
	of bcam_source_pair_type records.
}
function bcam_source_pair_spread(cp:pointer;num_calibs:integer):bcam_source_pair_type;

var
	calib_num,source_num:integer;
	calib_ptr:^bcam_source_pair_type;
	max,min,spread:bcam_source_pair_type;
	
begin
	for calib_num:=1 to num_calibs do begin
		calib_ptr:=pointer(integer(cp)+(calib_num-1)*sizeof(bcam_source_pair_type));
		if calib_num=1 then begin
			max:=calib_ptr^;
			min:=calib_ptr^;
		end
		else begin
			for source_num:=1 to num_sources_per_pair do begin
				with calib_ptr^ .sources[source_num] do begin
					if max.sources[source_num].x<x then max.sources[source_num].x:=x;
					if max.sources[source_num].y<y then max.sources[source_num].y:=y;
					if max.sources[source_num].z<z then max.sources[source_num].z:=z;
					if min.sources[source_num].x>x then min.sources[source_num].x:=x;
					if min.sources[source_num].y>y then min.sources[source_num].y:=y;
					if min.sources[source_num].z>z then min.sources[source_num].z:=z;
				end;
			end;
		end;
	end;
	with spread do begin
		for source_num:=1 to num_sources_per_pair do begin
			spread.sources[source_num]:=
				xyz_difference(max.sources[source_num],min.sources[source_num]);
		end;
	end;
	spread.id:='spread';
	bcam_source_pair_spread:=spread;
end;

{
	bcam_origin returns the origin of the bcam coordinates for the specified mounting balls.
}
function bcam_origin(mount:kinematic_mount_type):xyz_point_type;

begin
	bcam_origin:=mount.cone;
end;

{
	bcam_coordinates_from_camera_mount takes the global coordintes of the camera mounting balls 
	and calculates the origin and axis unit vectors of the bcam coordinate system expressed in global 
	coordinates.
}
function bcam_coordinates_from_camera_mount(mount:kinematic_mount_type):coordinates_type;
	
var
	bcam:coordinates_type;
	cs,cp,cs_normal:xyz_point_type;
	
begin
	with bcam,mount do begin
{
	The unit vectors cs and cp define a plane. The normal to this plane
	is the y-axis of the bcam coordinate system, as defined by the cross
	product of cp with cs.
}
		cs:=xyz_unit_vector(xyz_difference(slot,cone));
		cp:=xyz_unit_vector(xyz_difference(plane,cone));
		y_axis:=xyz_unit_vector(xyz_cross_product(cp,cs));
{
	The orientation of the bcam camera around its y-axis is set by the slot 
	and cone. We create the z-axis of the bcam system by rotating sc
	by z_angle about the y-axis. The resulting z-axis lies parallel to the 
	nominal optical axis of the camera. To perform the rotation, we define 
	cs_normal using cs and y_axis.
}
		cs_normal:=xyz_unit_vector(xyz_cross_product(cs,y_axis));
		z_axis:=
			xyz_unit_vector(
				xyz_sum(
					xyz_scale(cs,-cos(bcam_z_angle)),
					xyz_scale(cs_normal,-sin(bcam_z_angle))));
{
	The x-axis is the cross product of the y and z axes respectively.
}
		x_axis:=xyz_unit_vector(xyz_cross_product(y_axis,z_axis));
{
	We place the origin with the bcam_origin routine.
}
		origin:=bcam_origin(mount);
	end;
	bcam_coordinates_from_camera_mount:=bcam;
end;

{
	bcam_from_global_vector converts a direction in global coordinates into a 
	direction in bcam coordinates.
}
function bcam_from_global_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;

var
	M:xyz_matrix_type;
	bcam:coordinates_type;
	
begin
	bcam:=bcam_coordinates_from_camera_mount(mount);
	M:=xyz_matrix_from_points(bcam.x_axis,bcam.y_axis,bcam.z_axis);
	bcam_from_global_vector:=xyz_transform(M,p);
end;

{
	bcam_from_global_point converts a point in global coordinates into a point
	in bcam coordinates.
}
function bcam_from_global_point (p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;

begin
	bcam_from_global_point:=bcam_from_global_vector(xyz_difference(p,bcam_origin(mount)),mount);
end;

{
	bcam_from_global_z calculates the z-position of a source in bcam coordinates 
	given its z-position in global coordinates. The routine assumes that the source is on 
	the global z-axis. 
}
function bcam_from_global_z (z:real;mount:kinematic_mount_type):real;

var
	p,q:xyz_point_type;
	
begin
	p.x:=0;
	p.y:=0;
	p.z:=z;
	q:=bcam_from_global_point(p,mount);
	bcam_from_global_z:=q.z;
end;

{
	global_from_bcam_vector converts a direction in bcam coordinates into a 
	direction in global coordinates.
}
function global_from_bcam_vector(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;
var bc:coordinates_type;	
begin
	bc:=bcam_coordinates_from_camera_mount(mount);
	global_from_bcam_vector:=
		xyz_transform(
			xyz_matrix_inverse(
				xyz_matrix_from_points(bc.x_axis,bc.y_axis,bc.z_axis)),
			p);
end;

{
	global_from_bcam_point converts a point in bcam coordinates into a point
	in global coordinates.
}
function global_from_bcam_point(p:xyz_point_type;mount:kinematic_mount_type):xyz_point_type;

begin
	global_from_bcam_point:=xyz_sum(bcam_origin(mount),global_from_bcam_vector(p,mount));
end;

{
	global_from_bcam_line converts a bearing (point and direction) in bcam coordinates into
	a bearing in global coordinates.
}
function global_from_bcam_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type;

var
	gb:xyz_line_type;
	
begin
	gb.point:=global_from_bcam_point(b.point,mount);
	gb.direction:=global_from_bcam_vector(b.direction,mount);
	global_from_bcam_line:=gb;
end;

{
	bcam_from_global_line does the opposite of global_from_bcam_line
}
function bcam_from_global_line(b:xyz_line_type;mount:kinematic_mount_type):xyz_line_type;

var
	bb:xyz_line_type;
	
begin
	bb.point:=bcam_from_global_point(b.point,mount);
	bb.direction:=bcam_from_global_vector(b.direction,mount);
	bcam_from_global_line:=bb;
end;

{
	global_from_bcam_plane converts a bearing (point and direction) in bcam coordinates into
	a bearing in global coordinates.
}
function global_from_bcam_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type;

var
	gp:xyz_plane_type;
	
begin
	gp.point:=global_from_bcam_point(p.point,mount);
	gp.normal:=global_from_bcam_vector(p.normal,mount);
	global_from_bcam_plane:=gp;
end;

{
	bcam_from_global_plane does the opposite of global_from_bcam_plane
}
function bcam_from_global_plane(p:xyz_plane_type;mount:kinematic_mount_type):xyz_plane_type;

var
	bp:xyz_plane_type;
	
begin
	bp.point:=bcam_from_global_point(p.point,mount);
	bp.normal:=bcam_from_global_vector(p.normal,mount);
	bcam_from_global_plane:=bp;
end;

{
	Refer to Volume 8 page 107 for diagram. Suppose we have a bearing nearly parallel
	to the z-axis in bcam coordinates. We put the camera on two different mounts such
	that the bcam z-axis in both mounts is nearly parallel to the global z-axis. When we move
	the bcam from one mount to the other, the global direction of the bcam bearing changes. 
	Our calibration procedure measures the divergence of the two global bearings, and
	we call this divergance the 'link bearing'. Knowing the coordinates of the mounting balls 
	we can deduce the bearing in bcam coordinates.
}
function bcam_from_link_bearing(link_bearing:xyz_point_type;
	mount_1,mount_2:kinematic_mount_type):xyz_point_type;

var
	p:xyz_point_type;
	M_1,M_2,M:xyz_matrix_type;
	bc_1,bc_2:coordinates_type;
	
begin
{
	Calculate the bcam coordinates for the two mounting positions.
}
	bc_1:=bcam_coordinates_from_camera_mount(mount_1);
	bc_2:=bcam_coordinates_from_camera_mount(mount_2);
{
	Form the inverse of the coupling matrix.
}
	with bc_1 do M_1:=xyz_matrix_from_points(x_axis,y_axis,z_axis);
	with bc_2 do M_2:=xyz_matrix_from_points(x_axis,y_axis,z_axis);
	M:=xyz_matrix_difference(
			xyz_matrix_inverse(M_2),
			xyz_matrix_inverse(M_1));
{
	M is not necessarily invertible. The z-axis diagonal element tends to be zero.
	We want to invert it to get the coupling matrix as it applies to the other two
	axes. We can use the 3x3 matrix inversion if we set element [3,3] to one. We
	set the link vector z-component to 1 to match. Then we transform the link vector
	with our coupling matrix.
}
	M[num_xyz_dimensions,num_xyz_dimensions]:=1;
	link_bearing.z:=1;
	p:=xyz_transform(xyz_matrix_inverse(M),link_bearing);
{
	To complete the use of 3x3 matrices to perform a 2.5-dimensional calculation,
	we set the z-component of p to zero.
}
	p.z:=0;
	bcam_from_link_bearing:=p;
end;

{
	Refer to Volume 8 page 107 for diagram. Suppose we have a point in bcam
	coordinates. We put the camera on two different mounts such that the bcam z-axis
	in both mounts is nearly parallel to the global z-axis, and the bcam z=0 planes
	are nearly coincident. When we move the bcam from one mount to the other, the
	global position of our point changes. Our calibration procedure measures the
	global vector between the two locations of the point, in a z=constant plane that
	nearly contains both points. Knowing the coordinates of the mounting balls we
	can deduce the point's bcam coordinates. The 'link vector' is the vector in
	global coordinates that joins the two locations of the point. To implement
	bcam_from_link_vector, we use the bcam_from_link_bearing routine. Because the
	transformation from link vector to bcam point operates in a single z-plane, we
	must project all points involved in the calculation into this z-plane, including
	the origins of both bcam coordinates. We choose to project the bcam origin along
	the bcam z-axis so that the x- and y-components of our point, when calculated
	using the projected coordinate origins, will be equal to the x- and
	y-coordinates in the normal bcam coordinates. The z_plane parameter specifies
	the z-coordinate of the projection plane in bcam coordinates.
}
function bcam_from_link_vector(link_vector:xyz_point_type;
	mount_1,mount_2:kinematic_mount_type;
	z_plane:real)
	:xyz_point_type;

var
	origin_projected_into_pivot_plane:xyz_point_type;
	
begin
	with origin_projected_into_pivot_plane do begin
		x:=0;
		y:=0;
		z:=z_plane;
	end;
	bcam_from_link_vector:=
		bcam_from_link_bearing(
			xyz_difference(
				link_vector,
				xyz_difference(
					global_from_bcam_point(origin_projected_into_pivot_plane,mount_2),
					global_from_bcam_point(origin_projected_into_pivot_plane,mount_1))),
			mount_1,
			mount_2);			
end;

{
	ccd_center calculates the bcam coordinates of the ccd center.
}
function ccd_center(camera:bcam_camera_type):xyz_point_type;

begin
	with camera do 
		ccd_center:=xyz_sum(pivot,xyz_scale(axis,-ccd_to_pivot));
end;

{
	bcam_from_image_point converts a point on the ccd into a point in bcam coordinates. 
	The calculation takes account of the orientation of the ccd in the camera.
}
function bcam_from_image_point(p:xy_point_type;camera:bcam_camera_type):xyz_point_type;

var
	q:xy_point_type;
	r:xyz_point_type;
	cc:xyz_point_type;
	
begin
	q.x:=p.x-bcam_ccd_center_x;
	q.y:=p.y-bcam_ccd_center_y;
	cc:=ccd_center(camera);
	with camera do begin
		if axis.z>0 then begin
			r.x:=cc.x+q.x*cos(ccd_rotation)-q.y*sin(ccd_rotation);
			r.y:=cc.y+q.y*cos(ccd_rotation)+q.x*sin(ccd_rotation);
		end
		else begin
			r.x:=cc.x-q.x*cos(ccd_rotation)+q.y*sin(ccd_rotation);
			r.y:=cc.y+q.y*cos(ccd_rotation)+q.x*sin(ccd_rotation);
		end;
		r.z:=cc.z;
	end;
	bcam_from_image_point:=r;
end;

{
	image_from_bcam_point does the opposite of bcam_from_image_point. It assumes that 
	the z-coordinate bcam point is in the ccd.
}
function image_from_bcam_point(p:xyz_point_type;camera:bcam_camera_type):xy_point_type;

var
	q:xy_point_type;
	r:xy_point_type;
	cc:xyz_point_type;
	
begin
	cc:=ccd_center(camera);
	with camera do begin
		r.x:=p.x-cc.x;
		r.y:=p.y-cc.y;
		if axis.z>0 then begin
			q.x:=bcam_ccd_center_x+r.x*cos(ccd_rotation)+r.y*sin(ccd_rotation);
			q.y:=bcam_ccd_center_y+r.y*cos(ccd_rotation)-r.x*sin(ccd_rotation);
		end
		else begin
			q.x:=bcam_ccd_center_x-r.x*cos(ccd_rotation)+r.y*sin(ccd_rotation);
			q.y:=bcam_ccd_center_y+r.y*cos(ccd_rotation)+r.x*sin(ccd_rotation);
		end;
	end;
	image_from_bcam_point:=q;
end;

{
	bcam_source_bearing takes the position of a light spot on the ccd, and returns the pivot point
	and the direction of the source from the pivot point in bcam coordinates.
}
function bcam_source_bearing (spot_center:xy_point_type;camera:bcam_camera_type):xyz_line_type;

var
	bearing:xyz_line_type;
	image:xyz_point_type;
	
begin
	image:=bcam_from_image_point(spot_center,camera);
	bearing.point:=camera.pivot;
	bearing.direction:=xyz_unit_vector(xyz_difference(camera.pivot,image));
	bcam_source_bearing:=bearing;
end;

{
	bcam_source_position returns the bcam coordinates of a source, 
	given its bcam z-coordinate, its image center, and the camera
	calibration constants.
}
function bcam_source_position (spot_center:xy_point_type;range:real;camera:bcam_camera_type):xyz_point_type;

begin
	bcam_source_position:=
		xyz_line_plane_intersection(
			bcam_source_bearing(spot_center,camera),
			xyz_z_plane(range));
end;

{
	global_from_calib_datum takes a calib_datum_type and calculates the position of
	the source in global coordinates.
}
function global_from_calib_datum (p:calib_datum_type;
	camera:bcam_camera_type;
	mount:kinematic_mount_type):xyz_point_type;

begin
	global_from_calib_datum:=
		xyz_line_plane_intersection(
			global_from_bcam_line(bcam_source_bearing(p.spot_center,camera),mount),
			xyz_z_plane(p.source_range));
end;

{
	ccd_image_of_bcam_source takes a point in bcam coordinates and calculates 
	the theoretical position of its image on the ccd in image coordinates.
}
function ccd_image_of_bcam_source (p:xyz_point_type;
	camera:bcam_camera_type):xy_point_type;

var
	line:xyz_line_type;
	plane:xyz_plane_type;
	cc:xyz_point_type;

begin
	line.point:=camera.pivot;
	line.direction:=xyz_difference(p,line.point);
	cc:=ccd_center(camera);
	plane:=xyz_z_plane(cc.z);
	ccd_image_of_bcam_source:=
		image_from_bcam_point(
			xyz_line_plane_intersection(line,plane),
			camera);
end;

{
	global_offset_vector returns the vector in global coordinates that 
	joins the bcam axis to the source specified by the calib_datum_type.  
}
function global_offset_vector (p:calib_datum_type;
	camera:bcam_camera_type;
	mount:kinematic_mount_type):xyz_point_type;

var
	source_point,axis_point:xyz_point_type;
	axis_data_point:calib_datum_type;

begin
	source_point:=global_from_calib_datum(p,camera,mount);
	with axis_data_point do begin
		spot_center.x:=bcam_ccd_center_x;
		spot_center.y:=bcam_ccd_center_y;
		source_range:=p.source_range;
	end;
	axis_point:=global_from_calib_datum(axis_data_point,camera,mount);
	global_offset_vector:=xyz_difference(source_point,axis_point);
end;

{
	nominal_bcam_camera returns the nominal bcam_camera_type.
}
function nominal_bcam_camera(axis_direction:real):bcam_camera_type;

var
	camera:bcam_camera_type;
	
begin
	with camera do begin
		pivot:=xyz_origin;
		axis:=xyz_origin;
		if (axis_direction>=0) then begin
			axis.z:=+1;
			id:='nominal_fc';
		end
		else begin
			axis.z:=-1;
			id:='nominal_rc';
		end;
		ccd_to_pivot:=1;
		ccd_rotation:=0;
	end;{item}
	nominal_bcam_camera:=camera;
end;

{
}
function bcam_pair_calib(calib_data:pair_data_type):bcam_camera_type;

const
	close_range_num=1;
	far_range_num=2;
	first_mount_num=1;
	second_mount_num=2;
	show_details=false;
	fsn=20;
	fsd=6;
	
var
	offsets:array [1..num_mounts_per_pair,1..num_ranges_per_mount,1..num_sources_per_range] 
		of xyz_point_type;
	mount_num,source_num,second_source_num,range_num,count:integer;
	calibration:bcam_camera_type;
	divergance_global,divergance_bcam,camera_axis:xyz_point_type;
	r,s:xy_point_type;
	p,q,pivot_sum,axis_sum,xyz_datum:xyz_point_type;
	line:xyz_line_type;
	plane:xyz_plane_type;
	far_range,close_range,sum,datum,cos_sum,sin_sum,ns,fs:real;
	
begin
	with calib_data do begin
{
	Set the bcam calibration to its nominal values.
}
		calibration:=nominal_bcam_camera(axis_direction);
		calibration.id:=calib_data.id;
{
	Calculate the z-position of the pivot point.
}
		with calibration do begin
			if show_details then gui_writeln('pivot.z');
			sum:=0;
			count:=0;
			for mount_num:=1 to num_mounts_per_pair do begin
				for source_num:=1 to num_sources_per_range-1 do begin
					for second_source_num:=source_num+1 to num_sources_per_range do begin
{stop indentation}
ns:=xy_separation(
	measurements[mount_num,close_range_num,source_num].spot_center,
	measurements[mount_num,close_range_num,second_source_num].spot_center);
fs:=xy_separation(
	measurements[mount_num,far_range_num,source_num].spot_center,
	measurements[mount_num,far_range_num,second_source_num].spot_center);
datum:=
	bcam_from_global_z(
		measurements[mount_num,close_range_num,source_num].source_range-
			(measurements[mount_num,far_range_num,source_num].source_range
			-measurements[mount_num,close_range_num,source_num].source_range)
				/(ns/fs-1),
		mounts[mount_num]);
sum:=sum+datum;
inc(count);
{resume indentation}
						if show_details then begin 
							writestr(debug_string,id,mount_num:fsd,source_num:fsd,
								second_source_num:fsd,
								ns:fsn:fsd,
								fs:fsn:fsd,
								datum:fsn:fsd);
							gui_writeln(debug_string);
						end;
					end;
				end;
			end;
			pivot.z:=sum/count;
			if show_details then begin
				writestr(debug_string,sum/count:fsn:fsd);
				gui_writeln(debug_string);
			end;
		end;
{
	Calculate ccd to pivot distance. We use the value of pivot.z we
	obtained in the previous stage, even though the value we obtained is
	less accurate than the value we obtain from our drawings. The reason
	we use the calculated pivot.z is because it appears to express the
	errors in the calibration images. Our calibration of the ccd-pivot
	distance will be accurate despite these errors so long as we use our
	calculated pivot.z. We found this out by measuring the ccd-pivot and
	pivot.z directly on a 2-m granite beam with a ruler and a pair of
	sources. The pivot.z calculated in the previous stage was incorrect by
	2 mm, but if we used the correct value when calculating ccd-pivot, we
	obtained the wrong value for ccd-pivot. When we used our calculated
	pivot.z, however, we obtained a ccd-pivot distance accurate to better
	than 0.1%.
}
		with calibration do begin
			if show_details then gui_writeln('ccd_to_pivot');
			sum:=0;
			count:=0;
			for mount_num:=1 to num_mounts_per_pair do begin
				for range_num:=1 to 1 {num_ranges_per_mount} do begin
					for source_num:=1 to num_sources_per_range-1 do begin
						for second_source_num:=source_num+1 to num_sources_per_range do begin
{stop indentation}
datum:=axis.z*
	(bcam_from_global_z(
		measurements[mount_num,range_num,source_num].source_range,
		mounts[mount_num])
	-pivot.z)
	*xy_separation(measurements[mount_num,range_num,source_num].spot_center,
		measurements[mount_num,range_num,second_source_num].spot_center)
	/xy_separation(source_locations[source_num],source_locations[second_source_num]);
sum:=sum+datum;
inc(count);
{resume indentation}
							if show_details then begin
								writestr(debug_string,id,mount_num:fsd,range_num:fsd,
									source_num:fsd,second_source_num:fsd,datum:fsn:fsd);
								gui_writeln(debug_string);
							end;
						end;
					end;
				end;
			end;
			ccd_to_pivot:=sum/count;
			if show_details then begin
				writestr(debug_string,sum/count:fsn:fsd);
				gui_writeln(debug_string);
			end;
		end;
		
{
	Calculate rotation of CCD
}
		with calibration do begin
			if show_details then gui_writeln('ccd_rotation');
			sin_sum:=0;
			cos_sum:=0;
			count:=0;
			for mount_num:=1 to num_mounts_per_pair do begin
				for range_num:=1 to num_ranges_per_mount do begin
					for source_num:=1 to num_sources_per_range-1 do begin
						for second_source_num:=source_num+1 to num_sources_per_range do begin
{stop indentation}
p:=xyz_unit_vector(
		xyz_difference(
			global_from_calib_datum(measurements[mount_num,range_num,second_source_num],
				calibration,mounts[mount_num]),
			global_from_calib_datum(measurements[mount_num,range_num,source_num],
				calibration,mounts[mount_num])));
r:=xy_unit_vector(
		xy_difference(
			source_locations[second_source_num],
			source_locations[source_num]));
datum:=calibration.ccd_rotation-full_arctan(p.y,p.x)+full_arctan(r.y,r.x);
cos_sum:=cos_sum+cos(datum);
sin_sum:=sin_sum+sin(datum);
inc(count);
{resume indentation}
							if show_details then begin
								writestr(debug_string,id,mount_num:fsd,range_num:fsd,
									source_num:fsd,second_source_num:fsd,datum:fsn:fsd);
								gui_writeln(debug_string);
							end;
						end;
					end;
				end;
			end;
			ccd_rotation:=full_arctan(sin_sum/count,cos_sum/count);
			if (ccd_rotation < -pi/2) then ccd_rotation:=ccd_rotation+2*pi;
			if show_details then begin
				writestr(debug_string,ccd_rotation:fsn:fsd);
				gui_writeln(debug_string);
			end;
		end;
{
	Calculate the offset vectors in global coordinates that join the bcam axis to
	each source at each range in each mount. When we calculate these offsets, we use
	the nominal values of the pivot point position and the axis bearing. We use then
	use these offsets to calculate the actual pivot and axis.
}
		for mount_num:=1 to num_mounts_per_pair do begin
			for range_num:=1 to num_ranges_per_mount do begin
				for source_num:=1 to num_sources_per_range do begin
					offsets[mount_num,range_num,source_num]:=
						global_offset_vector(
							measurements[mount_num,range_num,source_num],
							calibration,
							mounts[mount_num]);
				end;
			end;
		end;
{
	Calculate the pivot to source range of the two source positions.
}
		count:=0;
		axis_sum:=xyz_origin;
		pivot_sum:=xyz_origin;
		if show_details then gui_writeln('axis and pivot');
		for source_num:=1 to num_sources_per_range do begin
{
	We assume that the range of the first source is the same for both mounts.
}
			close_range:=
				bcam_from_global_z(
					measurements[first_mount_num,close_range_num,source_num].source_range,
					mounts[first_mount_num])
				-calibration.pivot.z;
			far_range:=
				bcam_from_global_z(
					measurements[first_mount_num,far_range_num,source_num].source_range,
					mounts[first_mount_num])
				-calibration.pivot.z;
{
	With the camera on two different mounts, the lens center is in two positions,
	and the axis is in two orientations. We calculate the divergance of these two
	orientations in global coordinates by examining the offsets we calculated above.
	We express the divergence as a vector whose z- component is 0. It gives the
	divergance per meter moved along the global z-axis.
}
			divergance_global:=
				xyz_scale(
					xyz_difference(
						xyz_difference(
							offsets[first_mount_num,far_range_num,source_num],
							offsets[second_mount_num,far_range_num,source_num]),
						xyz_difference(
							offsets[first_mount_num,close_range_num,source_num],
							offsets[second_mount_num,close_range_num,source_num])),
					1/(far_range-close_range));
			divergance_global.z:=0; 
{
	Calculate the vector in bcam coordinates that gives rise to divergance_global, and convert
	this into a camera-axis unit vector. We assume the camera axis is nearly parallel to
	the global z-axis, and that the bcam z-axis of the bcam coordinates in both mounts is
	also nearly parallel to the z-axis. When the global z-axis and the bcam z-axis are in 
	opposite directions, we must negate the direction cosines of the calculated camera axis.
}
			xyz_datum:=xyz_scale(
				bcam_from_link_bearing(divergance_global,
					mounts[first_mount_num],
					mounts[second_mount_num]),
				calibration.axis.z);
			axis_sum:=xyz_sum(axis_sum,xyz_datum);
			if show_details then begin
				writestr(debug_string,id,source_num:fsd,xyz_datum.x:fsn:fsd,xyz_datum.y:fsn:fsd);
				gui_writeln(debug_string);
			end;
{
	Determine the vector between the two pivot positions and convert it into bcam coordinates 
	to get the actual lens position.
}
			xyz_datum:=bcam_from_link_vector(
					xyz_difference(
						xyz_difference(
							offsets[first_mount_num,close_range_num,source_num],
							offsets[second_mount_num,close_range_num,source_num]),
						xyz_scale(divergance_global,close_range)),
					mounts[first_mount_num],
					mounts[second_mount_num],
					calibration.pivot.z);
			pivot_sum:=xyz_sum(pivot_sum,xyz_datum);
			if show_details then begin
				writestr(debug_string,xyz_datum.x:fsn:fsd,xyz_datum.y:fsn:fsd);
				gui_writeln(debug_string);
			end;
{
	Done with this run through the loop so increment the data counter.
}
			inc(count);
		end;
		with calibration do begin
			pivot.x:=pivot_sum.x/count;
			pivot.y:=pivot_sum.y/count;
			axis.x:=axis_sum.x/count;
			axis.y:=axis_sum.y/count;
			if show_details then begin
				debug_string:=string_from_xyz(pivot);gui_writeln(debug_string);
				debug_string:=string_from_xyz(axis);gui_writeln(debug_string);
			end;
		end;
	end;
{
	Done.
}
	bcam_pair_calib:=calibration;
end;

{
	bcam_camera_calib takes a bcam_camera_calib_input_type and returns a
	bcam_camera_calib_output_type.
}
function bcam_camera_calib(calib_input:bcam_camera_calib_input_type):
	bcam_camera_calib_output_type;

var
	pair_data:pair_data_type;
	calib_output:bcam_camera_calib_output_type;
	mount_num,second_mount_num,source_num,range_num,pair_num:integer;
	p:xyz_point_type;
	
begin
{
	Apply apparatus description and camera name to pair data.
}
	pair_data.source_locations:=calib_input.source_locations;
	pair_data.axis_direction:=calib_input.axis_direction;
{
	Construct all possible pair_data_types for bcam_pair_calib and
	calculate calibration. 
}
	pair_num:=1;
	for mount_num:=1 to num_mounts_per_quad-1 do begin
		for second_mount_num:=mount_num+1 to num_mounts_per_quad do begin
			writestr(pair_data.id,mount_num:1,'_',second_mount_num:1);
			pair_data.mounts[1]:=calib_input.mounts[mount_num];
			pair_data.mounts[2]:=calib_input.mounts[second_mount_num];
			for range_num:=1 to num_ranges_per_mount do begin
				for source_num:=1 to num_sources_per_range do begin
					pair_data.measurements[1,range_num,source_num]:=
						calib_input.measurements[mount_num,range_num,source_num];
				end;
			end;
			for range_num:=1 to num_ranges_per_mount do begin
				for source_num:=1 to num_sources_per_range do begin
					pair_data.measurements[2,range_num,source_num]:=
						calib_input.measurements[second_mount_num,range_num,source_num];
				end;
			end;
			calib_output.pairs[pair_num]:=bcam_pair_calib(pair_data);
			inc(pair_num);
		end;
	end;
{
	Calculate spread and average parameters.
}
	with calib_output do begin
		average:=bcam_camera_average(@calib_output.pairs,num_roll_cage_pairs);
		spread:=bcam_camera_spread(@calib_output.pairs,num_roll_cage_pairs);
		id:=calib_input.id;
		time:=calib_input.time;
	end;
	bcam_camera_calib:=calib_output;
end;

{
	bcam_source_calib takes a bcam_source_pair_input_type and returns a
	bcam_source_pair_output_type.
}
function bcam_source_pair_calib(calib_input:bcam_source_pair_calib_input_type):
	bcam_source_pair_calib_output_type;

var
	first_mount_num,second_mount_num,source_num,pair_num:integer;
	p,s,global_link:xyz_point_type;
	viewing_link,v:xy_point_type;
	calib_output:bcam_source_pair_calib_output_type;
	
begin
{
	Calculate source positions.
}
	with calib_input do begin
		pair_num:=1;
		calib_output.id:=id;
		for first_mount_num:=1 to num_mounts_per_quad-1 do begin
			for second_mount_num:=first_mount_num+1 to num_mounts_per_quad do begin
				for source_num:=1 to num_sources_per_pair do begin
					viewing_link:=xy_difference(
						measurements[second_mount_num,source_num].spot_center,
						measurements[first_mount_num,source_num].spot_center);
{
	Rotate viwing link to account for orientation of viewing camera wrt global x.
}
					v:=xy_rotate(viewing_link,-viewing_x_direction/mrad_per_rad);
{
	Negate y-component because positive y in image is negative y in global coords.
}
					v.y:=-v.y;
{
	Get actual translation from ccd translation.
}
					v:=xy_scale(v,viewing_scale);
{
	Initialize global_link.
}
					with global_link do begin x:=v.x;y:=v.y;z:=0; end;
{
	Convert link into the bcam-coords of source position.
}
					p:=bcam_from_link_vector(
							global_link,
							mounts[first_mount_num],
							mounts[second_mount_num],
							measurements[first_mount_num,source_num].source_range);
{
	Restore z-value of position.
}
					p.z:=measurements[first_mount_num,source_num].source_range;
					calib_output.pairs[pair_num].sources[source_num]:=p;
				end;{for}
				writestr(calib_output.pairs[pair_num].id,
					first_mount_num:1,'_',second_mount_num:1);
				inc(pair_num);
			end;{for}
		end;{for}
	end;{with}
{
	Calculate spread and average parameters.
}
	with calib_output do begin
		average:=bcam_source_pair_average(@calib_output.pairs,num_roll_cage_pairs);
		spread:=bcam_source_pair_spread(@calib_output.pairs,num_roll_cage_pairs);
		id:=calib_input.id;
		time:=calib_input.time;
	end;
	bcam_source_pair_calib:=calib_output;
end;


end.