Write Up

Document Sample
Write Up Powered By Docstoc
					                            Patrick Gloster




     Summer Project 2008:
  Compact Straightness Monitor
Simulation/Multilateration Report
       By Patrick Gloster




                                    Page 1
                                                                        Patrick Gloster




                                      Contents



Section 1 – Accessing code and results………………………………………....Page 4

Section 2 – Introduction…………………………………………………………Page 5

Section 3 – The theory of least squares minimization…………………………..Page 5

Section 4 – Problem 1 - Finding the position of a target point using measurements
from 4 known observer points in 2D……………………………………….……Page 9

Section 5 – Problem 2 - Finding the position of a set of target points and the observer
points themselves using a set of measurements in 2D………………………….Page 12

Section 6 – Problem 3 – The bow-tie problem in 2D…………………………..Page 17

Section 7 – Problem 4 – The bow-tie problem in 3D…………………………..Page 28

Section 8 – Saving data………………………………………………………...Page 41

Section 9 – Current Matlab code…………………………………….…………Page 42




                                                                                Page 2
                                                                        Patrick Gloster


1)Accessing code and results

All of my files have been transferred from my H drive to the computer PPXP435.
They are all contained in a folder entitled ‘CSM Summer 2008’ on the C drive.

While writing this document I wrote instructions on how to access the results of
different programs. Since then I have reorganised the saved files and code. The
following should enable all data to be accessed:

      When Matlab starts up, choose the current directory (found centrally at the top
       of the screen) to be ‘C:\CSM Summer 2008\MATLAB\Code’ (Note – I’m
       not sure if you need the C at the beginning)
      To run a program type the name of the program followed by the input
       arguments in brackets, separated by commas. Each program is described in
       detail in this document, and each description tells you exactly how to start
       running the program, but they all have the same general format. For example,
       to run the program bowtie5.m with 100 plate 2 positions and save it to a file
       named ‘bowtie5_10_figure1.bmp’ type ‘bowtie5(100,10)’ in the Matlab
       command window and press return
      To look at the code click on the open folder icon, and then choose the Code
       folder, followed by whichever program you wish to see the code for.
      All programs with automatic saves will save files to the Matlab Code folder. If
       you haven’t moved the saved files from this folder after you generate results,
       they can be accessed as written in the text
      Previous results have been moved into separate folders within Matlab. To
       access the text files, first switch directory to ‘C:\CSM Summer
       2008\MATLAB’, then type ‘type FOLDER/FILENAME.TXT’ into the
       command window. For example, the results of bowtie 5 are in the bowtie5
       folder. To access the information from the 102nd run, make sure you are in the
       Matlab directory, then type ‘type bowtie5/bowtie5_102.txt’
      Picture files are saved as bmp files and can be accessed by clicking the open
       folder icon, followed by the folder which has the name of the program they
       were created by. If the folder looks empty it may be because files of type is set
       to ‘Matlab files’, in which case change it to ‘all files’
      The available folders containing text files are bowtie, bowtie2, bowtie3,
       bowtie4, bowtie5, bowtie6, bowtie7, bowtie8, bowtiehistogram, newbowtie3,
       newbowtie4 and newbowtiehistogram. Their names indicate what programs
       produced the results they contain
      The other folders contain figures. ‘2D bowtie varying angle’, 2D bowtie
       varying initial error’, 2D bowtie varying plate range’ and ‘2D bowtie varying
       platesize’ also contain workspace variables which can be imported into
       Matlab. These variables are vectors and a matrix containing the results of tests
       described towards the end of section 5 in this document
      Currently the results of the 3D programs are NOT saved automatically. I have
       been saving the histograms manually in designated folders. If you wish to save
       information automatically, follow the format of lines 204-255 in bowtie5. Also
       see the section on saving data in this document

The folder ‘CSM Summer 2008’ also contains a folder called ‘Excel data’ which has a
spreadsheet containing data generated when trying to find the optimum position of
plate 2 in the 2D bowtie set-up. This is referenced on page 24 of this write up.
                                                                                 Page 3
                                                                         Patrick Gloster



Finally there is a folder named ‘Write ups’, which contains various documents written
as the project progressed including odd things I sent to Patrick Brockill, a copy of this
write up, and some files which I used to write down some things as they happened.
All useful information has been included in this write up, but I have left these other
documents intact. Pages 12 and 13 of the file ‘locating observer coordinates’ contain
some information on allocating weights to different errors and may come in useful in
the future.

If the Matlab help section is not helping I suggest using the website:

http://www.mathworks.fr/matlabcentral/




                                                                                  Page 4
                                                                         Patrick Gloster



2. Introduction

The goal of this project was to create a program to calculate the positions of one
endplate of the compact straightness monitor with respect to the other. One endplate
has 4 launch pairs, arranged in a cross. The other has 4 retro-reflectors, also arranged
in a cross. The positions are calculated using the method of least squares. The
program also had to include a simulation of positions, so that we could compare
calculated values with true values of the positions.

The project was approached by writing several programs for different versions of the
problem, with gradually increasing complexity. For completeness I include a full
account of all of the code written, but many programs are no longer used at all. For
each section I will point out the programs which are the end products of that section.
The reason for so many discarded programs is the diverse approach taken for each
version of the problem. For the first section I worked almost exclusively with my own
least squares minimization function; however as the project progressed I
experimented with various minimization functions already present in Matlab, and
finally settled on a non-linear least squares minimization function called lsqnonlin.

3. The Theory of Least Squares Minimization

If we have a set of measurements in a vector, L , which fit a set of formulas involving
the position variables in the vector X , we can write our equation as A X  L where A
is some matrix (the design matrix), operating on the position vector X to give us our
measurements. Due to the limits on the precision of measuring instruments, taking
redundant physical measurements gives inconsistent results. This means no unique
solution exists; all we can do is make a unique estimate of the solution. The chosen
method is the method of least squares.

We can specify an initial approximation to the solution vector, X 0 , which we perturb
slightly to give the total solution vector X  X 0  X

Since all measurements are fundamentally inaccurate, however, the measured lengths
L will give inconsistent results. We therefore add a ‘residual’ vector V which will
cancel out the inconsistencies. Thus A X  L is satisfied where L  L  V

The method of least squares tells us that the best estimate for X is the estimate which
minimises the sum of the squares of the inconsistencies, i.e. minimises V T V




                                                                                  Page 5
                                                                                     Patrick Gloster



The weight matrix

Our measurements (which make up the elements of L) will not necessarily all have
the same precision. This is accounted for by assigning a known ‘weight’ matrix to
each measurement. P is the matrix whose elements are these weights. In this case the
least squares criterion becomes

V T PV = minimum


Linearization of the problem

Mathematical models expressing the relationship between the total observed and
solution vectors X and L have the general form:

F ( X , L)  0

Linearization of these models is accomplished by replacing the nonlinear functions F
by their Taylor’s series linear approximation, expanded about the point defined by the
initial approximation to the solution vector (X0), and the measured values of the
observation vector (L).

Remembering that X  X 0  X and L  L  V

linearization gives:

                              F ( X  X 0 , L  L)        F ( X  X 0 , L  L)
F ( X , L)  F ( X 0 , L)                            X                           V 0
                                      X                           L

Or W  AX  BV  0

                                   F ( X  X 0 , L  L)              F ( X  X 0 , L  L)
Where W  F ( X 0 , L) , A                                 and B 
                                            X                                 L

In our problem the observed quantities can be explicitly expressed as functions of the
parameters X (i.e. we can write our problem as F ( X )  L )

In this case, linearization gives us

                          F ( X  X 0 )
F(X )  L  F(X 0 )                       X  (L  V )  0
                                 X
Or W  AX  V  0
                                                                  L
Where W  F ( X 0 )  L




                                                                                              Page 6
                                                                           Patrick Gloster


                   F1        F1        
                                       
                   X 1       X 2       
   F ( X  X 0 )  F2                   
A               =                      
        X         X 1                  
                                      
                                       
                                         

Our problem is now:

V T PV = minimum
X’s subject to constraint: W  AX  V  0


To implement this constraint we use Lagrange multipliers.

Lagrange Multipliers

Suppose we want to minimize some function f1 ( x, y) subject to some constraint
 f 2 ( x, y)  0

The Lagrange method consists of 3 steps:

   1) Form the variation function

  f1 ( x, y)  kf 2 ( x, y) where k is an undetermined constant called the Lagrange
multiplier

   2) Set the derivatives of the variation function to zero

      
    0,    0
x      y


   3) Solve the three equations

                         
f 2 ( x, y)  0,       0,    0
                   x      y


In our case we have


 V T PV = minimum
 X’s subject to constraint:
 W  AX  V  0


So we form the variation function

                                                                                    Page 7
                                                                       Patrick Gloster



  V T PV  2 K T ( AX  V  W )

Here K is a column vector of Lagrange multipliers, and has been multiplied by a
factor of 2 for convenience. The derivatives are:


    2V T P  2K T  0
V

which becomes

PV  K  0

and


    2K T A  0
X

which becomes

AT K  0

So the three equations we want to solve are

W  AX  V  0
PV  K  0
AT K  0

These equations can be combined into a single hypermatrix equation (a hypermatrix
being a matrix whose elements are themselves matrices)




 P    BT    0  V   0 
              
 I   0     A   K   W   0
                    
 0    AT    0  X   0 
                 

This can be partitioned (see pages 114-116 of Wells and Krakiwsky) to give

X  ( AT PA) 1 AT PW
K  P( AX  W )
V  P 1 K  AX  W

The solution we are interested in here is X. Thus, if we can determine the matrices A,
P and W we can solve for X.



                                                                                Page 8
                                                                                     Patrick Gloster


        4. Problem 1 – Finding the position of a target point using measurements from 4
        known observer points in 2D

        The first problem was set up as follows:



Figure 1 – 4
observer points
with measured
distances to the
target point
(x,y)




        Here we are measuring the distance of a particular point from four ‘observer points’.
        We know the positions of these four observers. Our goal was to write a simulation
        which could generate a large number of target points and, for each of these points,
        generate measurements from the target point to the 4 observer coordinates, which
        would be inexact (as all physical measurements are). The method of least squares was
        then used to calculate an approximate position for the target point and the results were
        represented graphically.

        I only created one simulation for this problem, and used my own minimization
        function. I used the weight matrix

        A=
                 ( x0  x1 )                         ( y 0  y1 )               
                                                                                
         ( x0  x1 ) 2  ( y 0     y1 ) 2   ( x0  x1 ) 2  ( y 0  y1 ) 2     
                 ( x0  x 2 )                       ( y0  y2 )                 
                                                                                
         ( x0  x 2 ) 2  ( y 0    y2 ) 2   ( x0  x 2 ) 2  ( y 0  y 2 ) 2   
                                                                                
                 ( x 0  x3 )                       ( y0  y3 )                 
         (x  x )2  ( y 0         y3 ) 2   ( x 0  x3 ) 2  ( y 0  y 3 ) 2   
            0     3
                                                                                 
                 ( x0  x1 )                         ( y 0  y1 )               
        
         ( x0  x 4 ) 2  ( y 0                                                 
                                                                                 
                                   y4 ) 2   ( x0  x 4 ) 2  ( y 0  y 4 ) 2   

        x0 and y0 were approximations to the target position, generated in the simulation. You
        can choose how many points you want to generate and calculate positions of.

        The program calls several functions while it is running, so what follows is a
        description of each of these functions, followed by a description of the simulation
        itself.




                                                                                             Page 9
                                                                          Patrick Gloster


F.m

The F function takes two vectors describing two positions, and calculates the distance
between the two positions. The input must be in vector form. For example, to
calculate the distance between the points (2,4) and (0,0) you need to define two
vectors X1=[2;4] and X2=[0;0]. By typing ‘F(X1,X2)’ in the command window the
distance between these two points will be returned.

Locate.m

The locate.m function is used to find the location of an unknown point using two
measurements and three known locations. It has an input in the Matlab command
window of ‘locate(x1,y1,L1,x2,y2,L2,x3,y3)’. These are 3 sets of coordinates
describing the 3 known points, and measurements from two of these points to the
unknown point. Two possible positions for the target (unknown) point are calculated
using the cosine rule. Whichever of these two calculated points is closer to the third
known point is selected as the approximate location of the unknown point. As an
example, if you have 3 known points, (0,0), (5,3) and (6,8), and you want to know the
approximate location of a point which is 2 units from (0,0) and 5 units from (5,3) you
would type locate(0,0,2,5,3,5,6,8). One early problem with this function was what
happened if there was no real solution. This happened if both of the measurements
were less than the true distance. In this case the program returns the point halfway
between the two known points as the solution.

Blackcircle.m, bluecircle.m, greencircle.m, redcircle.m

These functions plot a circle of radius r centred at (x,y) in the colour included in the
name. So to plot a red circle centred on (4,2) with a radius of 1, type ‘redcircle(4,2,1)
in the Matlab command window.

Simulation.m

To begin the simulation, type ‘simulation(N)’ in the Matlab command window, where
N is the number of simulations you want to run (that is, how many target positions
you want to generate and calculate). This starts the program.

First a series of GUIs are generated prompting the user to choose the positions of their
known observer coordinates. The user is then given the option to jump to a simulation
number of their choice, and asked to input the errors on each measurement. This is so
that the weight matrix P can be created.

There follows a small loop in the code generating a series of random numbers. This is
done because the seed of the random number generator used in this program is set to
0. If you want to jump to a particular simulation, the correct number of random
numbers needs to be generated beforehand, so that the correct target point is
generated. This loop takes care of this by generating the same number of numbers as
would be done if you had run through all of the previous simulations rather than
jumping to that one.

For each simulation the program generates a target position somewhere in a 10 x 10
square. The true distances from each observer point to the target are then calculated,

                                                                                  Page 10
                                                                          Patrick Gloster


and measurements are generated by adding a different random number (normally
distributed about 0 with a standard deviation of 1) multiplied by the error for that
measurement (entered by the user earlier in the program) to each of the
measurements.

The simulation then generates an approximate position to start the minimization from.
This is done using the locate function described above. A figure is plotted showing
circles centred on the four observer coordinates, with radii corresponding to the
measurements.

The minimization proceeds by applying the general formula above;
 X  ( AT PA) 1 AT PW . This is done several times, converging to a solution. At the
beginning of the program I define a variable moddx. The iterations of the
minimization continue until the modulus of the perturbation to the solution vector
drops below moddx. There is a limit on iterations of 1000000 to prevent the program
crashing if it is not converging. The user is given the option to plot each successive
solution for the simulation, or to jump to the final solution.




                                                         Figure 2 – Results of the program
                                                         simulation.m. The user can zoom in
                                                         on successive iterations




                                                                                 Page 11
                                                                          Patrick Gloster


Gauss.m, gauss2.m

There were also two functions written to plot histograms, gauss.m and gauss2.m. The
first, gauss.m, is identical to the simulation but at the end it calculates the difference
between the true and calculated point for every simulation, and then plots a histogram
of the differences. To use this function type ‘gauss(N,simulationnumber)’ into the
command window, where N is the number of target points you want to generate and
simulationnumber is the simulation you want to jump to. The second is more useful,
and is run by typing ‘gauss2(N)’ into the command window, where N is again the
number of target points you want to generate. The program runs the simulation and
then calculates the distance from each observer point to the calculated point. It then
calculates the difference between these 4 distances and the 4 true distances. This is
done for each simulated target point. After all of the simulations 4 histograms are
plotted of the differences between true and calculated distances to each observer
point.




                                                                    Figure 3 – An example of 4
                                                                    histograms of the differences
                                                                    between calculated and true
                                                                    distances, produced by gauss2.m.




Notice the bottom two histograms having a wider spread than the top two. This is
interesting and something which I could never adequately explain. However, it was
felt that rather than dwelling on this problem we should move on to a new set-up.

5. Problem 2 – Finding the position of a set of target points and the observer
points themselves using a set of measurements in 2D

The next stage in the project was to introduce uncertainty on where the observer
points themselves were, and use many measurements to different target points to
calculate where these observer points were (and where the target points were). In
problem 1, we solved for each simulated target point separately, and had just 2
unknowns (the x and y coordinates of the target point). In this problem, however, we
had to solve for all of the simulated target points simultaneously, since it was only by
considering all of the sets of measurements at once that we could solve for the
positions of the observer points.

For example, we could generate 300 target positions; this would give us 600
unknowns (300 x coordinates, 300 y coordinates). We also have the 8 unknown
coordinates describing the 4 observer points. Since we can arbitrarily define an origin
                                                                                  Page 12
                                                                                    Patrick Gloster


        and a direction in this problem, we can reduce the number of unknown coordinates
        describing the 4 observer points to 5 by setting one of the points to be the origin,
        (0,0), and one to be along the y-axis, so at a position (0,y). In the example we are then
        left with 605 unknowns. However, for each generated target position we have 4
        measurements, one from each of the observer coordinates. This therefore gives us
        1200 measurements in total in our example. Since we have more measurements than
        unknowns we can again use the least squares method to give us a solution for all of
        the variables.

                           (0,y2)                                        (x3,y3)

Figure 4 – The set-
up for problem 2.
Here we have 3 of
the 8 variables set to
zero so we only have
to work with 5
unknowns



                                                                       (x4,y4)
                             (0,0)


        At this stage in the project I began to attack the problem on several fronts. Where
        before I had just used my own minimization function, I now started to use some of
        Matlab’s inbuilt functions and compare results. This led to many different functions
        being written, which are described below; while all are of some interest, the primary
        fucntions which I would use to find solutions to this problem are basicpat.m and its
        associated objective function somefun.m.

        Basicpat.m

        To use this function type ‘basicpat(N)’ in the Matlab command window, where N is
        the number of target points you wish to generate.

        This function generates a set of 5 true coordinates describing the positions of the
        observer points. We only generate 5 because we have a built-in assumption that one
        of the observer points is at (0,0) and the second is along the y-axis relative to this first
        point. The function then generates 2N coordinates to describe the N target points. By
        adding random errors to all of the target and observer coordinates the initial values are
        obtained to go into the vector X0.

        The function then calculates the actual distances from each observer point to each
        target point, and generates our measurements by adding an error to these distances.
        The error on the measurements (1E-5) is much less than the error used to give us our
        initial values in X0 (0.3).

        An anonymous function is then defined which calls a function called somefun, giving
        it the number of target points, the measurements, and the current values in the vector

                                                                                           Page 13
                                                                       Patrick Gloster


X. This anonymous function is necessary because all of Matlab’s inbuilt functions can
only have one input, so we define an anonymous function which only has one input
(in this case X) and calls another function which has several inputs (in this case X,
measurements and N). We then use Matlab’s function to minimize the anonymous
function.

The function then uses a built-in function called lsqnonlin to minimize the anonymous
function and give us a solution for the vector X. We use lsqnonlin since it is a non-
linear least squares minimizer.


Somefun.m

The function called in basicpat.m is somefun.m. This takes the vector X, which
contains our current answers for all variables, and works out what measurements we
would get if these were the values. It then works out the difference between these and
the actual measurements we have. By minimizing this we are therefore minimizing
the difference between the true measurements and these ‘predicted’ measurements.
Note – lsqnonlin sums and squares the residuals for us, which is why all somefun does
is calculate a vector of the differences.

Findobservertwod.m

 Here we come to a program using my own least squares minimization. It can be run
by typing ‘findobservertwod(N)’ in the Matlab command window. It solves the same
problem as basicpat, with 5 unknown observer coordinates. This was encouraging
because it gave the same answers as basicpat, indicating that they were both working.
However, it was much slower, which is why we decided to use Matlab’s least squares
minimizer instead. In writing this function I needed to write another function which
could calculate the symbolic form of the design matrix A. This was done in the
function makejacobian.m. Below is a comparison of the results using basicpat and
findobservertwod. Note the time taken for each function.

Makejacobian.m

This function is called in findobservertwod.m, and returns the symbolic form of the
matrix A. The values of the vector X are then substituted into the matrix. This
program proved difficult to write because I needed to be able to define an arbitrary
number of symbolic variables, which is a non-trivial problem in Matlab. Eventually
this was solved by asking for help on the internet forum Matlab central.

                    Y2          X3         Y3         X4         Y4         Time
                                                                            taken
                                                                            (seconds)
True value       0.219          0.047      0.678      0.679      0.934      -
Lsqnonlin        0..2204        0.0462     0.6827     0.6764     0.9406     0.186
calculated value
Findobservertwod 0.2204         0.0463     0.6827     0.6764     0.9405     112.5
calculated value



                                                                              Page 14
                                                                                         Patrick Gloster




               Basicpat2.m

               This function is identical to the previous function, basicpat, except it uses a different
               inbuilt minimizer called fminunc. This is not a least squares minimizer though;
               instead it is equipped with several types of search and therefore, while it is more
               powerful, it is not necessary to use it in our problem. I experimented with using this
               function only as a step to using fmincon, a minimization function which can handle
               constraints. To run this function type basicpat2(N) in the command window where N
               is the number of target points you want to generate.

               Somefun2.m

               Basicpat2 calls a different objective function to basicpat because, unlike lsqnonlin,
               fminunc does not accept an objective function which returns a vector. Somefun2
               therefore returns a scalar which is the sum of the squares of the difference between
               the ‘predicted’ measurements and the actual measurements. It is this that is
               minimized.

               Pat.m

               This function was an extension of the function basicpat.m. To use it type pat(N,M) in
               the Matlab command window, where N is the number of targets and M is the number
               of simulations you want to perform. The program performs the minimization using
               lsqnonlin, and repeats the entire thing M times. For each simulation it calculates the
               difference between the true position and the calculated position
               (e.g. true x3-calculated x3), and also the residuals of the objective function for each
               coordinate. Histograms are then plotted for all 10 of these, showing the data from all
               of the simulations.

               Basicpat3.m

               I wrote this function as an experiment in using constraints. Instead of including just 5
               variables for the coordinates of the observer points, I included 6, and applied the
               constraint that one of them was 0 as part of the minimization. This was done using a
               Matlab function called fmincon. The function works in a very similar way to
               basicpat.m and basicpat2.m, except where they have x2=0 built in as a condition, and
               therefore don’t work with it as a variable, basicpat3 includes it as a variable but then
               applies the constraint that x2=0 during the minimization.
                               (x2,y2)
                                                                           (x3,y3)



Figure 5 – The set-up for
the program basicpat3. We
include x2 as an unknown
variable and constrain it to
be zero during the
minimization


                                                                                                 Page 15
                                                                         (x4,y4)
                                 (0,0)
                                                                               Patrick Gloster


This program was a success; provided the true position of x2 was indeed close to zero
the answers were correct. Even if the true position of x2 was not zero, the constraint
was still applied; however, this led to huge inaccuracies in all of the other results.
Since basicpat3 works with 6 variables for the observer points instead of 5, I had to
use a different objective function again, so I wrote somefun3.m. I also included a
matrix in the program which defines the constraints. This has to be passed to fmincon
as well as the function to be minimized.

Here we see how the constraints work if we are applying a sensible one (i.e. saying
that y2=0 when it is close to 0) but if y2 is nowhere near 0 applying the constraint
destroys all accuracy in the calculated values:


                  Y2              X3             Y3              X4              Y4
True              0.0010          0.0470         0.6789          0.9347          0.3835
position
Calculated        0               0.0476         0.6717          0.9303          0.3923
position

                  Y2              X3             Y3              X4              Y4
True              0.2190          0.0470         0.6789          0.9347          0.3835
position
Calculated        0               0.2238         -0.5201         -0.1776         0.9441
position


Somefun3.m

This is the objective function for basicpat3. It is different to somefun1 and somefun2
because when working out the predicted measurements it includes the variable x2,
where before this was set to zero. So where before the distance from the target point
(x,y) to point 2 (x2,y2) was      ( x 2  ( y  y 2) 2 , somefun3 calculates the distance using
  ( x  x 2) 2  ( y  y 2) 2 .


Newpat.m

This function was written to test the use of fmincon without constraints. It is identical
to basicpat2.m, working with just 5 variables for the observer coordinates, only it uses
fmincon instead of fminunc.


Findobservertwod2.m & findobservertwod3.m

Both of these functions are solving the same problem as basicpat3 (so 6 unknown
observer coordinates) using my own minimization. Again, they call functions to give
them symbolic design matrices (named makejacobian2.m and makejacobian3.m
respectively). They differ from findobservertwod.m because they apply constraints.
These programs were both successful in calculating answers and applying constraints,
but were abandoned due to the length of time they took.

                                                                                      Page 16
                                                                                 Patrick Gloster



        6. Problem 3 – The bow-tie problem in 2D

        The bow-tie set up is a 2D version of the real experimental apparatus. We have two
        launch pairs on a fixed plate, and use them to measure the distances to two retro-
        reflectors on a second plate. It is called the bow-tie set up because the 4 measurements
        form a bow-tie shape. By moving around the right hand plate (plate 2) to different
        positions and obtaining a large number of measurements we can determine the
        coordinates specifying the locations of both plates.
                                                                            1m
                                        (a2,b2)               L2


    y              10 cm                                                                     (x2,y2)
                                                         L4




               x                                         L3


z                                     (a1,b1)                 L1                                   (x1,y1)


             Figure 6 – The bow-tie set up in 2D with point source launch pairs


        In the diagram above the left hand plate (which contains the launch heads for the
        lasers) is fixed; the right hand plate (which contains the retroreflectors) can move to
        the left and right, up and down, and rotate about the z axis. There is an error on the
        laser measurements of 1 nm. The error on our initial guesses for positions is 1mm. For
        the first program, I modelled the launch pairs as point sources, which is why we only
        have two points on the left hand plate. By defining the bottom point as the origin, and
        the top point as lying along the y-axis, we obtain x1=0, y1=0, x2=0. Distance scales
        are marked on the diagram. The distance between the two reflector points on the right
        hand plate is constrained to be the same for every position of this plate.

        An initial attempt to write a function to solve this problem using my own
        minimization method was quickly abandoned due to incorrect results. Similarly, when
        I wrote a function to solve this using fmincon, Matlab’s minimization function which
        allows the application of constraints, we found that the results were only accurate to
        approximately one tenth of a millimetre. We also found that increasing the number of
        target points did not improve the accuracy. We thus decided to reformulate the
        problem in a way which did not require constraints, so that we could use lsqnonlin.




                                                                                        Page 17
                                                                        Patrick Gloster


Reformulation of the problem

To reformulate the problem I changed the way I defined the position of the second
plate. Where before it was defined using the 4 coordinates x1,y1,x2,y2, it is now
defined using midx, midy and alpha, where midx and midy are the coordinates of the
centre of the plate, and alpha is the angle of orientation of the plate.

Old formulation                                                            New formulation




                        (x2,y2)
                                             d=Distance
                                             between reflectors            α


                                                                                         (midx,midy)



      (x1,y1)

                  Figure 7 – The old and new
                  formulations of the bow tie problem



Before I worked with 4 unknowns for every plate position, x1, y1, x2, y2, and
constrained the distance between these two points to be the same for every position of
the plate. With the new formulation, there are 3 unknowns for every plate position,
midx, midy and alpha. There is also d, which is one variable for all positions (i.e.
there is not a new d for every plate position). From midx, midy, alpha and d we can
work out the positions x1, y1, x2 and y2. We can therefore solve for midx, midy and
alpha and never have to apply any constraints during the minimization; this is already
part of the new formulation since in working out x1, y1, x2 and y2 (which we need to
do so that we can calculate ‘predicted’ measurements) we always use the same d.

As before there are several similar programs for this problem. Each one adds a little
complexity to the set-up, so the final programs are bowtie5.m and
new_bowtie5_histogram.m.

Bowtie2.m

This program can be run by typing ‘bowtie2(N,filenumber)’ in the command window,
where N is the number of different plate positions you want to generate and
filenumber is the number of the file you want to save the results as.

The program has b2 set to be 0.1, and the distance between reflector points to be 0.1.
The other coordinates of plate 1 are already set to be 0, so we only have one
coordinate to determine on this plate. The program then generates our initial values

                                                                               Page 18
                                                                          Patrick Gloster


for these variables by adding a normally distributed random error to them. The error
on these initial values is 1mm.

Using random numbers I then generate the midx and midy values of the second plate
N times. I also generate an angle of orientation of plate 2 for each position. Using
these I calculate the true positions of the reflectors for each plate position. Again, by
adding a random error in mm the program generates our initial values for the midpoint
and angle of the plate.

The program then uses the true values for the positions of the reflectors to calculate
the distances to the observer points, and adds a random error to give us the
measurements. The error on the measurements is 1nm.

It is again necessary to use an anonymous function in this program. The anonymous
function calls muchfun2.m, the objective function, which calculates predicted
measurements based on the current values of all of the variables. The anonymous
function is then minimized using the Matlab minimization function lsqnonlin. To
achieve high accuracy on the results I set one of the optimization parameters, TolFun
to 10-18. The lsqnonlin function then calculates values for all of the variables – b2 and
then midx, midy and alpha for N different plate positions.

The program then calculates the difference between the true values and the calculated
values for midx, midy and alpha and calculates the standard deviation of each. It also
plots a histogram of these differences for each of the three parameters.

Finally the program saves all information in a text file called myfile_filenumber.txt,
where filenumber is the second number you entered in the brackets in the command
window. To access this information you need to type ‘type myfile0filenumber.txt’ in
the command window. For example, if you enter bowtie2(50,4) in the command
window, to see all of the information saved you need to enter ‘type myfile04.txt’ into
the command window. The information saved includes a table of initial values,
calculated values and true values for midx, midy and alpha, the exit flag, the number
of iterations, the number of plate positions, the real value of b2, the calculated value
of b2, the real separation of reflectors and the calculated separation of reflectors. The
histograms are also saved under 'myfile0filenumber_figure1.bmp'. These pictures can
then be found in the bowtie2 folder within the Matlab folder in ‘CSM Summer 2008’.
One problem experienced with saving the pictures is that it seems to save them by
taking a print screen. This means that if you are looking at something else while the
program runs, and it gets to the point in the program when it should save the
histograms, it may end up saving whatever is on the screen at the time.

Muchfun2

This is the objective function minimized in bowtie2. It uses the current values of
midx, midy and alpha to calculate current values of x1, y1, x2 and y2, and then uses
these and the current value of b2 to calculate predicted measurements. Finally it finds
the difference between these predicted measurements and the true measurements. This
is what is minimized.




                                                                                 Page 19
                                                                            Patrick Gloster


    Switching to finite sized launch pairs on plate 1

    In bowtie2.m the lasers are launched from point sources. The next step was to
    separate them so that each of the four beams was emitted from a different source, and
    to allow them to protrude slightly from the surface of the plate. We use a separation of
    the launch heads of 2cm and a protrusion of 2mm.



    Notice that in this set-up the cross measurements (L3 and L4) are performed by the
    inside launch heads. This is something that was changed for the next function.




                                                                          1m
                                          (0,b2)             L2


    y            10 cm
                            2cm
                                                                                            (x2,y2)
                                                        L4
                              (a4,b4)




             x                (a3,b3)
                                                        L3
                                        2mm
z                                       (0,0)                L1                                 (x1,y1)


           Figure 8 – The bow-tie set up in 2D with separated launch heads




    Bowtie4.m

    This program runs in the same way as bowtie2.m, but it includes a3, b3, a4 and b4.
    Due to the inclusion of this it also needs a different objective function, the aptly
    named muchfun4.m. As before, to run the program type bowtie4(N,filenumber) in the
    command window. Information and results are saved as bowtie4_filenumber.txt, and
    can be accessed by typing ‘type bowtie4_filenumber.txt’ in the Matlab command
    window. Histograms are saved as 'bowtie4_filenumber_figure1.bmp' and can be
    found in the Matlab folder in ‘CSM Summer 2008’.

    new_bowtie4_histogram.m

    I wrote new_bowtie4_histogram.m so that we could examine the results for the
    positions of the launch heads more accurately. This is run by typing
    ‘new_bowtie4_histogram(N,n,filenumber)’. This works by generating the true values
    for N plate positions, as before, and then generating n sets of initial positions and

                                                                                    Page 20
                                                                             Patrick Gloster


    measurements. The program then solves for each of these n sets of measurements. It
    then calculates the differences between true and calculated values for b2, a3, b3, a4,
    b4, midx, midy and alpha. We have n values for the observer points (a’s and b’s), and
    n x N values for plate 2 parameters (midx’s, midy’s and alphas). Histograms are then
    plotted for all of these variables. The histograms are automatically saved as
    'new_bowtie4_run(filenumber)_histogram.bmp'.



    Measuring cross distances with the outside launch heads

    In the previous programs I measured the cross distances L3 and L4 with the inside
    launch heads. The next step was to switch to measuring these with the outside launch
    heads and compare results.



                                                                             1m
                                            (0,b2)
                                                              L2

    y              10 cm
                              2cm
                                                                                                (x2,y2)
                                                         L4
                                (a4,b4)




               x                (a3,b3)                    L3

                                          2mm
z                                                               L1
                                          (0,0)                                                     (x1,y1)


             Figure 9 – The bow-tie set up in 2D with crossed beams


    Bowtie5.m

    This program solves the system with the new set-up. Files are saved as
    bowtie5_filenumber.txt and histograms are saved as
    bowtie5_filenumber_figure1.bmp.

    To compare the results of this with the results of bowtie4, which is exactly the same
    but with the outer heads doing the straight measurements and the inner heads doing
    the cross measurements, I set all true values to be identical, used the same seed for the
    random number generator for both, used the same number of targets for each, and put
    exactly the same true error on each measurement. I obtained the following results:




                                                                                     Page 21
                                                                         Patrick Gloster


Program         File          Seed           Midx standard          Midy standard
                number                       deviation              deviation
Bowtie4         72            0              1.94E-6                2.49E-6
Bowtie5         18            0              8.03E-10               7.55E-10
Bowtie4         73            12             8.76E-7                1.15E-6
Bowtie5         19            12             3.06E-9                5.32E-9
Bowtie4         74            678            3.82E-6                5.20E-6
Bowtie5         20            678            2.90E-9                5.52E-9
Bowtie4         75            1008           9.01E-10               2.43E-9
Bowtie5         21            1008           3.78E-9                5.93E-9
Bowtie4         76            66.6           1.94E-10               2.87E-9
Bowtie5         22            66.6           3.28E-9                6.26E-9

We see that while the results for bowtie4 fluctuate wildly, the results for bowtie5 are
consistently extremely good. We also notice that the results are much more accurate
in this unrealistic case, when the error on each measurement is exactly the same (by
which I mean the actual difference between the measurement and the true value, not
just the size of the possible error). This was further examined and produced the
following results:


Errors            File number        Seed             Midx               Midy
                                                      standard           Standard
                                                      deviation          deviation
Identical         18                 0                8.03E-10           7.55E-10
Non-identical     23                 0                3.46E-9            4.18E-8
Identical         19                 12               3.06E-9            5.32E-9
Non-identical     24                 12               6.02E-8            1.23E-7
Identical         20                 678              2.90E-9            5.52E-9
Non-identical     25                 678              2.10E-8            4.76E-8
Identical         21                 1008             3.78E-9            5.93E-9
Non-identical     26                 1008             2.30E-8            1.05E-7
Identical         22                 66.6             3.28E-9            6.26E-9
Non-identical     27                 66.6             1.68E-8            2.12E-8




                                                                                Page 22
                                                                          Patrick Gloster




Figure 10 - An example of the histograms plotted in bowtie5


  The mean and the standard deviation

  Upon closer inspection it was found that the results given by bowtie5 were not what
  we would expect; the standard deviation was of the order of nanometres, but the
  average size of the difference between the true and calculated values was significantly
  larger than this standard deviation. We found gaussian distributions which were not
  centred on 0. The table below shows some results; seed is the number given to the
  random number generator to act as the seed. The error on initial measurements was
  1mm. 200 targets were used.

  File number            Seed                    Mean x                  Standard deviation
                                                                         on x
  73                     2                       1.87E-9                 1.94E-9
  74                     8                       1.52E-7                 2.75E-8
  75                     42                      6.10E-7                 8.38E-8
  76                     4096                    -5.37E-9                3.40E-9
  77                     11                      5.58E-9                 1.87E-9




                                                                                 Page 23
                                                                            Patrick Gloster


   It was hypothesised that this was due to poor calculated values for the launch points
   on plate 1, i.e. b2, a3, b3, a4 and b4. This was supported when we ran the program
   without any error on these parameters, and without including them as variables to be
   minimized. When this was done the gaussians were distributed about a mean much
   closer to zero than the standard deviation:




    Changing the error on initial values
Figure 11 – When we had no uncertainty on the positions of the launch
heads on plate 1, we got the distribution we expected

   For the majority of the project I worked with an error of 1mm on our ‘initial
   estimates’; that is, the error added to the true positions to give us our starting point
   (which we fed into the minimization program lsqnonlin) was of the order of 1E-3.
   When I experimented with this value, it produced some shocking results which I have
   not yet found an explanation for. In each case the same seed, 0, was used. Again there
   were 200 targets.

   File           Error on       Mean x         Standard       Mean y         Standard
   number         initial                       deviation                     deviation
                  values                        on x                          on y
   113            1E-3           1.28E-10       1.36E-10       7.95E-8        1.52E-8
   114            1E-5           -9.93E-9       5.89E-9        -7.30E-7       1.89E-7
   115            1E-7           -3.93E-10      1.20E-9        3.54E-8        1.33E-8


                                                                                   Page 24
                                                                          Patrick Gloster


This was done with the expectation of nothing happening, since no matter what
position you start from I would have expected the final convergence to be the same.
However, previous tests (page 80 of the logbook), which had some implicit errors and
are thus not included here, support this trend.

Further tests were carried out using 50 repeats (rather than just looking at one sample
with the same seed each time), and again showed a large and unexpected change in
accuracy. Graphs of these results are saved in the Matlab folders ‘bowtie5’.


Finding the optimum range of movement of plate 2

The next step was to use the program bowtie5.m and vary the positions of the launch
heads, the angle which the second plate could be orientated at, and the distance the
second plate could move in the x and y directions, in an effort to determine the
optimum arrangement.

Initially I proceeded by looking at 5 samples, each with a different seed, and every
time I changed a parameter (i.e. the angle or the distance that plate 2 could have) I
recorded various values for each of these 5 samples. The values recorded were

   a) the difference between the calculated value for b2 and the true value
   b) the difference between the calculated value for a3 and the true value
   c) the difference between the calculated value for b3 and the true value
   d) the difference between the calculated value for a4 and the true value
   e) the difference between the calculated value for b4 and the true value
   f) the mean of the differences between the calculated and true midx values
   g) the standard deviation of the differences between the calculated and true midx
      values
   h) the mean of the differences between the calculated and true midy values
   i) the standard deviation of the differences between the calculated and true midy
      values

This information is saved in EXCEL in a file called OPTIMIZATION TABLES.

This test was, however, limited in scope; looking at 5 samples wasn’t enough. I
therefore used the program new_bowtie5_histogram.m to look at 200 targets with 50
different random seeds. This program is similar to the other new_bowtie_histogram.m
files; for each of the 50 repeats the same ‘true’ 200 plate positions are used, but each
time a different random error is added to the positions (to give us our initial estimates)
and to the distances (to give us our measurements). Using this method I varied the
range of angles the plate could have, the range of movement of the plate and the size
of the plates (keeping the plates the same size as each other for each case).

The results of these tests are saved as variables and can be imported into the Matlab
workspace. The variables from the test of varying the range of motion that plate 2 can
have in the x and y direction are saved as varying_range.m. The variables from the
various angle test are saved as varying_angle.m and the variables of changing the size
of the plates are saved as varying_platesize.m.. I also used this method to investigate
changing the the error on the initial estimates; these variables are saved under
varying_initial_error.m.

                                                                                 Page 25
                                                                        Patrick Gloster



In each of these cases the results are saved in a matrix: range_matrix, angle_matrix,
platesize_matrix or error_matrix. Each column of the matrix represents a different
parameter (each column represents the same thing in each matrix). They are defined
below:

            Column number                                   Parameter
                  1                                       Range (+/- m)
                  2                                      Angle (radians)
                  3                            Mean difference between true b2 and
                                                       calculated value (m)
                     4                          Standard deviation of b2 difference
                     5                                Mean of a3 difference
                     6                          Standard deviation of a3 difference
                     7                                Mean of b3 difference
                     8                          Standard deviation of b3 difference
                     9                                Mean of a4 difference
                    10                          Standard deviation of a4 difference
                    11                                Mean of b4 difference
                    12                          Standard deviation of b4 difference
                    13                        Mean of mean difference between true
                                                     and calculated x values *
                    14                       Standard deviation of mean x differences
                    15                           Mean of standard deviation of x
                                                            differences
                    16                       Standard deviation of standard deviation
                                                         of x differences
                    17                            Mean of mean of y differences
                    18                           Standard deviation of mean of y
                                                            differences
                    19                           Mean of standard deviation of y
                                                            differences
                    20                       Standard deviation of standard deviation
                                                         of y differences
                    21                          Mean of mean of alpha differences
                    22                         Standard deviation of mean of alpha
                                                            differences
                    23                         Mean of standard deviation of alpha
                                                            differences
                    24                       Standard deviation of standard deviation
                                                       of alpha differences
                    25                               Mean of all x differences
                    26                        Standard deviation of all x differences
                    27                               Mean of all y differences
                    28                        Standard deviation of all y differences
                    29                             Mean of all alpha differences
                    30                            Standard deviation of all alpha
                                                            differences
                    31                                      Initial error


                                                                               Page 26
                                                                        Patrick Gloster


* We have these odd values because for each of the 50 repeats, we work out an
average difference in both x and y, and we also calculate the standard deviation of
these differences. In the matrices we therefore include the averages and standard
deviations of all of these averages and standard deviations, and also include the mean
and standard deviation of all of the values across all of the repeats.

Graphs plotting parameters 3-30 against range, angle and initial error are saved as
Matlab figures e.g. to see how the mean difference between the true and calculated
values for b2 changed with the amount that the second plate could move see the graph
entitled ‘b2_mean_vs_range.fig’.

The results show a general trend of increasing accuracy with increasing angle, range
and plate size. As we continue to increase these parameters still further the accuracy
decreases sharply, but I believe this is due to the optimization no longer being stopped
by TolFun; instead it is halted by TolX. It is difficult to make specific statements
about which values of these parameters would be the best to use because increasing
angle from 10 degrees to 20 degrees could for example make the accuracy in working
out a2 better, but worsen the accuraracy of b4.

Sparsity pattern

When implementing the new_bowtie5_histogram.m program it could be very slow
due to the 50 repetitions, and the fact that we were varying the parameters to values
which led to extremely long convergence times. This prompted me to add an
ammendement to the code which returns the jacobian for the first convergence. All of
the non-zero elements of this jacobian are then converted to ones, giving us a matrix
representing the sparsity pattern of the jacobian. In every subsequent use of lsqnonlin
this sparsity pattern is handed to the function. Using the sparsity pattern makes
convergence approximately 10 times faster, thus although the first convergence is
slow, the rest of the program runs very quickly.




                                                                                Page 27
                                                                           Patrick Gloster


  7. Problem 4 – The bow-tie problem in 3D

  Plate 1 is modelled as having 8 launch heads:
                                                                    Looking at plate 1 from
                                                                    plate 2
                          (a4,b4,c4)



                                       (a3,b3,c3)
             (a6,b6,c6)                                (a8,b8,c8)


     (a5,b5,c5)                                (a7,b7,c7)



                                       (a2,b2,c2)



                                       (a1,b1,c1)



  I define a1=0, b1=0, c1=0, a4=0, c4=0, c5=0, giving us a plane in z=0.

  Plate 2 lies in the positive z direction:




      10cm




Plate 1
                                                                                Plate 2
                                         z




                                                                                  Page 28
                                                                                   Patrick Gloster


Plate 2 is modelled as having 4 retroreflectors:



                                                                               Looking at plate 2 from
               (x2,y2,z2) = point 2                                            plate 1



                                                d1
                                                            (x4,y4,z4) = point 4
                         d2

      (x3,y3,z3) = point 3
                                                                                    z          y

                                                                                    x


                                          (x1,y1,z1) = point 1




Plate 2 can rotate around 3 axes:

1) Rotation about the x-axis where the angle alpha is the angle between the plate and
the y-axis
                                                              y
                x           y
                                                       α

                                      z                                        z




2) Rotation about the y-axis where beta is the angle between plate 2 and the x-axis

                                                                    x
       x             y
                                                           β

                              z                                            z


3) Rotation about the z-axis where gamma is the angle between plate 2 and the x-axis

                                                           y
           x             y
                                                                    γ

                                  z                                       x



                                                                                          Page 29
                                                                       Patrick Gloster


If we define the midpoint of plate 2 as (midx, midy, midz), then using angles alpha,
beta and gamma, and the distances d1 and d2, we can define the locations of the 4
retroreflectors:

x1 = midx + (d1/2).cosα.sinγ

y1 = midy – (d1/2).cosα.cosγ

z1 = midz + (d1/2).sinα


x2 = midx – (d1/2).cosα.sinγ

y2 = midy + (d1/2).cosα.cosγ

z2 = midz – (d1/2).sinα


x3 = midx + (d2/2).cosβ.cosγ

y3 = midy + (d2/2).cosβ.sinγ

z3 = midz – (d2/2).sinβ


x4 = midx – (d2/2).cosβ.cosγ

y4 = midy – (d2/2).cosβ.sinγ

z4 = midz + (d2/2).sinβ


We define our 8 measurements as follows:

L1  ( x1  a 2) 2  ( y1  b2) 2  ( z1  c 2) 2
L2  ( x1  a 4) 2  ( y1  b4) 2  ( z1  c 4) 2
L3  ( x 2  a3) 2  ( y 2  b3) 2  ( z 2  c3) 2
L4  ( x 2  a1) 2  ( y 2  b1) 2  ( z 2  c1) 2
L5  ( x3  a7) 2  ( y3  b7) 2  ( z3  c7) 2
L6  ( x3  a5) 2  ( y3  b5) 2  ( z3  c5) 2
L7  ( x4  a6) 2  ( y 4  b6) 2  ( z 4  c6) 2
L8  ( x 4  a8) 2  ( y 4  b8) 2  ( z 4  c8) 2

Remembering that a1=0, b1=0, c1=0, a4=0, c4=0 and c5=0.




                                                                               Page 30
                                                                          Patrick Gloster


threeDbowtie.m

The program threeDbowtie.m is the program to solve for the positions of the plates in
three dimensions. It is very much an extension of the 2D programs; where before I
defined 4 points on plate 1 and 2 on plate 2, I now define 8 points on plate 1 and 4 on
plate 2. Again the objective function calculates the difference between the predicted
measurements and the simulated measurements, and is minimized. It is run by
entering ‘threeDbowtie(N,n)’ in the Matlab command window, where N is the
number of plate 2 positions and n is the number of times you want to generate these
points and minimise them. In previous programs which included repeats the same N
true positions were used in all of the repeats; it was only the initial estimates and the
measurements which changed. In this program, however, an entirely new set of N
plate 2 positions is generated. To look at one run you just set n=1.

NOTE – An important point to make is that the variables rangeinmidX, rangeinmidY,
rangeinmidZ are the amount that can be added or taken away from the starting (x,y,z)
value 1m away from the centre of plate 1. So if rangeinmidZ is set to 0.1m the middle
of plate 2 can vary from 0.9m away from the centre of plate 1 to 1.1m away from the
centre of plate 1. This is the same in the following programs as well.

The manner in which the positions of plate 2 are generated is as follows:

   1) We take a point 1m away from the centre of plate 1 along the Z-axis
   2) We add a random number between -1 and 1 multiplied by the variables
      rangeinmidX/Y/Z to this point to give us the centre of plate 2


ThreeDbowtie.m minimizes all of the variables at once, and plots histograms of the
differences between the true and calculated values. The next step is to do run the
minimization in two steps; first we find values for the positions of the launch heads on
plate 1, then we use these calculated values as fixed positions and can calculate the
positions of plate 2 with a systematic error.

It should be noted that for this program and for calibrationbowtie.m I have an error on
initial estimates of angle of 1 degree. This can be changed in line 19 of the m-file.

calibrationbowtie3.m

This program takes a set of N plate 2 positions and calculates the positions of the
launch heads on plate 1. These are then used for the rest of the simulation. A further
M plate 2 positions are generated and the program calculates the values of midx,
midy, midz, alpha, beta and gamma for these positions, always using the same
calculated positions for the a’s, b’s and c’s. This program is run by typing
‘calibrationbowtie3(N,M) in the Matlab command window, where N is the number of
plate positions you want to use to determine the calibration constants and M is the
number of plate positions you want to calculate the positions of using these
calibration constants.

My initial set up of the launch heads on plate 1 was completely symmetric. This gave
results of:


                                                                                 Page 31
                                                                       Patrick Gloster


          midX           midY         midZ         Alpha        Beta         Gamma
          difference     difference   difference   difference   difference   difference
Mean      1.70E-6        -7.82E-6     1.63E-6      -2.75E-6     -8.34E-6     1.52E-4
Standard 3.59E-6         2.62E-6      5.66E-8      7.22E-5      8.34E-5      1.06E-3
deviation

We note that gamma is particularly bad.

The positions of the launch heads making up one of the launch pairs was then lowered
5 mm relative to the other pair, so where before b5=0.05, b6=0.05, b7=0.05 and
b8=0.05 we now had b5=0.045, b6=0.045, b7=0.05 and b8=0.05. This gave the
following results:


          midX           midY         midZ         Alpha        Beta         Gamma
          difference     difference   difference   difference   difference   difference
Mean      1.29E-7        -3.49E-7     6.09E-8      -3.67E-7     1.16E-8      8.76E-7
Standard 4.08E-8         4.88E-8      1.10E-9      2.56E-6      2.44E-6      3.22E-6
deviation

And when repeated:


          midX           midY         midZ         Alpha        Beta         Gamma
          difference     difference   difference   difference   difference   difference
Mean      -7.50E-8       -2.46E-7     4.60E08      8.51E-10     4.91E-8      1.53E-6
Standard 3.96E-8         3.79E-8      9.98E-10     1.84E-6      1.77E-6      3.19E-6
deviation

This indicated that removing the symmetry greatly improved the accuracy of the
calibrations.

Calibrations.m

The next target was to find out how many targets were required to give us a consistent
set of calibrations. This was done using the program calibrations.m, which calculates
the calibration constants 50 times and plots histograms of the difference between the
true and calculated values.To run calibrations.m type ‘calibrations(N)’ into the Matlab
command window where N is the number of plate 2 positions you wish to generate.

Histograms are saved in the folder ‘Calibration histograms’ in the Matlab folder.




                                                                               Page 32
                                                                             Patrick Gloster


             100 plate 2 positions:

              a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference
    Mean      -2.62E-8      -6.80E-9      -9.73E-8      -6.97E-8      4.26E-8       9.81E-8
    Standard 8.72E-8        9.08E-8       3.74E-7       2.21E-7       2.42E-7       3.77E-7
    deviation

          b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference
Mean      3.92E-8        1.53E-7      1.88E-7       9.14E-8       1.00E-7       8.13E-8        9.18E-8
Standard 1.70E-7         6.03E-7      7.51E-7       3.61E-7       3.79E-7       3.56E-7        3.53E-7
deviation

                       c2_difference c3_difference c6_difference c7_difference c8_difference
             Mean      -1.76E-8      -1.75E-8      -1.80E-8      -1.71E-8      9.31E-10
             Standard 6.72E-8        7.30E-8       6.97E-8       7.05E-8       6.01E-9
             deviation



             200 plate 2 positions:


              a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference
    Mean      7.20E-9       1.21E-8       4.67E-9       8.57E-9       7.77E-9       2.25E-9
    Standard 6.05E-8        6.06E-8       1.04E-7       7.37E-8       7.72E-8       1.06E-7
    deviation

          b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference
Mean      1.46E-8        8.59E-9      4.38E-9       -5.72E-11     1.30E-8       7.02E-9        2.66E-10
Standard 5.98E-8         1.70E-7      2.06E-7       1.01E-7       1.25E-7       1.09E-7        9.69E-8
deviation

                       c2_difference c3_difference c6_difference c7_difference c8_difference
             Mean      1.02E-9       1.43E-10      8.10E-10      2.53E-10      3.12E-10
             Standard 1.83E-8        2.01E-8       1.92E-8       1.94E-8       4.06E-9
             deviation




                                                                                    Page 33
                                                                             Patrick Gloster


             300 plate 2 positions:

              a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference
    Mean      -1.15E-9      4.22E-9       -7.71E-9      -3.34E-9      8.13E-9       7.24E-9
    Standard 4.72E-8        5.73E-8       5.48E-8       4.12E-8       4.11E-8       5.36E-8
    deviation

          b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference
Mean      8.67E-9        1.55E-8      2.24E-8       1.07E-8       1.41E-8       9.88E-9        1.01E-8
Standard 4.22E-8         9.45E-8      1.13E-7       5.50E-8       6.44E-8       8.62E-8        5.40E-8
deviation

                       c2_difference c3_difference c6_difference c7_difference c8_difference
             Mean      -1.12E-9      -1.93E-9      -8.32E-10     -2.27E-9      -9.09E-10
             Standard 1.04E-8        1.03E-8       9.94E-9       1.08E-8       2.87E-9
             deviation

             400 plate 2 positions:

              a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference
    Mean      -2.02E-9      -4.88E-9      5.54E-10      1.15E-9       -5.73E-9      1.66E-9
    Standard 3.49E-8        4.62E-8       4.19E-8       3.71E-8       4.29E-8       4.27E-8
    deviation

          b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference
Mean      6.70E-9        -3.24E-9     4.63E-10      -4.34E-10     1.38E-9       3.86E-9        3.42E-11
Standard 3.01E-8         6.78E-8      8.47E-8       3.95E-8       4.69E-8       5.75E-8        3.99E-8
deviation

                       c2_difference c3_difference c6_difference c7_difference c8_difference
             Mean      1.71E-10      4.57E-10      4.31E-10      1.77E-10      2.42E-10
             Standard  7.48E-9       8.68E-9       7.76E-9       8.64E-9       2.66E-9
             deviation

             500 plate 2 positions:

              a2_difference a3_difference a5_difference a6_difference a7_difference a8_difference
    Mean      -5.80E-9      -4.99E-10     -4.83E-10     1.88E-9       -4.10E-9      1.58E-9
    Standard 3.70E-8        3.76E-8       4.35E-8       3.73E-8       3.04E-8       4.20E-8
    deviation

          b2_difference b3_difference b4_difference b5_difference b6_difference b7_difference b8_difference
Mean      2.53E-9        -4.95E-9     -1.25E-9      -1.32E-9      2.05E-9       -3.99E-9       -1.36E-9
Standard 2.97E-8         6.46E-8      7.83E-8       3.90E-8       5.44E-8       4.97E-8        3.88E-8
deviation

                       c2_difference c3_difference c6_difference c7_difference c8_difference
             Mean      3.04E-10      1.13E-10      5.84E-11      2.99E-10      7.18E-11
             Standard 6.99E-9        7.55E-9       7.36E-9       7.09E-9       2.36E-9
             deviation

                                                                                    Page 34
                                                                       Patrick Gloster



When attempting to do 600 plate 2 positions Matlab gave an error ‘Out of memory’.
We can see that the accuracy of the locations of the launch heads is increasing with an
increasing number of plate 2 positions, so it would be advisable to use as many plate 2
positions as possible. However, time is an issue. When using calibrationbowtie3 to
calculate the calibration constants using 500 targets, and then using these calculated
values to calculate the position of 100 plate 2 positions, the time taken was 416
seconds.




Introducing further offsets on the launch heads for calibrationbowtie3.m

We saw earlier that lowering one of the launch pairs with respect to its opposite
results in a far better accuracy when calculating the positions of plate 2 due to the
removal of the symmetry. Subsequent testing should include adding further offsets,
with launch heads sticking out of the plate varying amounts and different distances to
the centre of the plate for each launch head.

The previous settings for calibrationbowtie3.m were:

a1                0           b1            0              c1              0
a2                0           b2            0.02           c2              0
a3                0           b3            0.08           c3              0
a4                0           b4            0.1            c4              0
a5                -0.05       b5            0.045          c5              0
a6                -0.03       b6            0.045          c6              0
a7                0.03        b7            0.05           c7              0
a8                0.05        b8            0.05           c8              0

giving results:


          midX            midY         midZ         Alpha        Beta          Gamma
          difference      difference   difference   difference   difference    difference
Mean      1.29E-7         -3.49E-7     6.09E-8      -3.67E-7     1.16E-8       8.76E-7
Standard 4.08E-8          4.88E-8      1.10E-9      2.56E-6      2.44E-6       3.22E-6
deviation

And when repeated:


          midX            midY         midZ         Alpha        Beta          Gamma
          difference      difference   difference   difference   difference    difference
Mean      -7.50E-8        -2.46E-7     4.60E08      8.51E-10     4.91E-8       1.53E-6
Standard 3.96E-8          3.79E-8      9.98E-10     1.84E-6      1.77E-6       3.19E-6
deviation




                                                                                Page 35
                                                                                      Patrick Gloster


The next test used the following settings for the launch head positions:

a1                 0                    b1                0                c1            0
a2                 0                    b2                0.02             c2            0
a3                 0                    b3                0.08             c3            0
a4                 0                    b4                0.1              c4            0
a5                 -0.05                b5                0.045            c5            0
a6                 -0.03                b6                0.045            c6            0
a7                 0.028                b7                0.05             c7            0
a8                 0.048                b8                0.05             c8            0

I have highlighted in red the alterations. We have shifted the two launch pair on the
right of the diagram 2mm towards the centre. This means the launch heads on the
hoizontal are no longer equidistant from the centre.

                                                                                   Looking at plate 1 from
                                                                                   plate 2
                           (a4,b4,c4)



                                             (a3,b3,c3)
            (a6,b6,c6)                                            (a8,b8,c8)


      (a5,b5,c5)                                     (a7,b7,c7)



                                             (a2,b2,c2)



                                             (a1,b1,c1)




The results for this set up were:


          midX                  midY           midZ               Alpha         Beta         Gamma
          difference            difference     difference         difference    difference   difference
Mean      -1.86E-9              -1.07E-7       1.47E-8            -2.43E-8      2.01E-8      4.22E-7
Standard 1.80E-8                1.48E-8        5.09E-10           6.44E-7       6.91E-7      1.67E-7
deviation

Figure saved as ‘figure1’ in the calibrationbowtie3 folder in the Matlab folder.

When repeated:


             midX       midY       midZ       Alpha      Beta       Gamma
             difference difference difference difference difference difference
                                                                                              Page 36
                                                                         Patrick Gloster


Mean      -7.91E-9       -6.22E-8      1.05E-8       -6.14E-9      2.21E-8       4.26E-7
Standard 1.50E-8         1.56E-8       4.19E-10      5.37E-7       4.27E-7       1.58E-6
deviation

Figure saved as ‘figure2’ in the calibrationbowtie3 folder.

We see that this appears to have increased accuracy further when compared to the
previous results shown on page 34.



The next test was to introduce an asymmetry in the vertical components on plate 1. I
chose to shift the upper launch pair closer to the centre of the plate. However, I did
not try intoducing a horizontal shift because the top launch head (the launch head at
(a4, b4, c4)) is defined to be at (0, b4, 0). While it is possible to change this, and
worth testing in the future, this must be done carefully because it would involve
rewriting the objective function; this is because currently the objective function has
the value of a4=0 hardwired into it. The new value would have to replace this in the
objective function if we varied a4.

The new set-up was therefore:


a1             0              b1             0                c1             0
a2             0              b2             0.02             c2             0
a3             0              b3             0.078            c3             0
a4             0              b4             0.098            c4             0
a5             -0.05          b5             0.045            c5             0
a6             -0.03          b6             0.045            c6             0
a7             0.028          b7             0.05             c7             0
a8             0.048          b8             0.05             c8             0


This produced the following results:

          midX           midY          midZ          Alpha         Beta          Gamma
          difference     difference    difference    difference    difference    difference
Mean      1.77E-8        -9.52E-8      5.28E-9       -1.25E-7      5.68E-8       4.99E-7
Standard 1.43E-8         1.46E-8       6.42E-10      1.43E-7       1.90E-7       1.28E-6
deviation

Figure saved as ‘figure3’ in the calibrationbowtie3 folder in the Matlab folder.

When repeated:


          midX           midY          midZ          Alpha         Beta          Gamma
          difference     difference    difference    difference    difference    difference
Mean      -1.86E-9       -1.07E-7      1.47E-8       -2.43E-8      2.01E-8       4.22E-7
Standard 1.80E-8         1.48E-8       5.09E-10      7.44E-7       6.91E-7       1.67E-6
deviation
                                                                                   Page 37
                                                                        Patrick Gloster



Figure saved as ‘figure4’ in the calibrationbowtie3 folder in the Matlab folder.


This has not had any obvious effects. It should be noted that these are preliminary
tests to give us an idea of any change; rigorous testing would involve the use of
calibrationbowtie4 with these difference settings for a’s, b’s and c’s.



The next test was introducing some protrusion from the plate:

a1             0              b1             0              c1             0
a2             0              b2             0.02           c2             0
a3             0              b3             0.078          c3             0
a4             0              b4             0.098          c4             0
a5             -0.05          b5             0.045          c5             0
a6             -0.03          b6             0.045          c6             0
a7             0.028          b7             0.05           c7             0.002
a8             0.048          b8             0.05           c8             0.002




          midX           midY         midZ           Alpha        Beta         Gamma
          difference     difference   difference     difference   difference   difference
Mean      3.69E-9        -4.66E-9     -3.66E-9       -2.74E-8     6.35E-8      -2.10E-7
Standard 1.25E-8         1.29E-8      4.80E-10       1.63E-7      1.59E-7      1.98E-6
deviation

Figure saved as ‘figure5’ in the calibrationbowtie3 folder in the Matlab folder.

When repeated:


          midX           midY         midZ           Alpha        Beta         Gamma
          difference     difference   difference     difference   difference   difference
Mean      2.19E-8        -6.18E-9     4.56E-9        3.96E-8      2.82E-8      5.65E-9
Standard 1.28E-8         1.20E-8      4.49E-10       2.01E-7      2.57E-7      1.02E-6
deviation

Figure saved as ‘figure6’ in the calibrationbowtie3 folder in the Matlab folder.


Again, there are no huge differences and this needs to be investigated using
calibrationbowtie4.




Time issues
                                                                                   Page 38
                                                                              Patrick Gloster



A problem I am currently experiencing with my programs is the amount of time it
takes to run each one. Since more plate 2 positions corresponds to a better calibration,
we want to use a large number of them. However, even using just 300 targets takes
131 seconds per calibration. This means looking at 50 calibrations using
calibrationbowtie4.m will be incredibly time consuming. If it is possible to speed up
the calibration stage this would be extremely beneficial.

Due to these problems with the amount of time taken for each calibration, I am
looking at calibrationbowtie4 with just 300 plate 2 positions for the calibration, and a
further 100 plate 2 positions calculated with that set of calibration constants, for 50
different calibrations. In the future we would like to have more plate 2 positions for
the calibration, but this is the largest number of positions which I can do in the time
available.

The true positions of the launch heads on plate 1 were set up as follows:

a1             0              b1             0                 c1               0
a2             0              b2             0.02              c2               0
a3             0              b3             0.08              c3               0
a4             0              b4             0.1               c4               0
a5             -0.05          b5             0.045             c5               0
a6             -0.03          b6             0.045             c6               0
a7             0.03           b7             0.05              c7               0
a8             0.05           b8             0.05              c8               0



We obtained the following results:

             midX         midY         midZ          alpha          beta            gamma
mean of      3.02E-9      -2.80E-8     3.99E-9       -3.10E-9       -1.00E-8        -4.04E-8
mean
standard     4.42E-8      6.12E-8      9.42E-9       6.71E-8        8.75E-8         4.37E-7
deviation
of mean
mean of      1.52E-8      1.56E-8      5.79E-10      3.39E-7        3.32E-7         1.72E-6
standard
deviation
standard     2.41E-9      2.77E-9      1.60E-10      3.02E-7        2.78E-7         3.30E-7
deviation
of
standard
deviation


Figures are saved in the folder calibrationbowtie4 as means1 and SD1.

If we look at the mean and standard deviation of the standard deviations we can see
that once we have calibrated plate 1 the difference between the calculated and true
values has a very consistent standard deviation of less than a nanometre on midZ and
                                                                                      Page 39
                                                                          Patrick Gloster


about 15 or 16nm on midX and midY. The average disparity between the true and
calculated values is approximately 30nm on midX and midY, and 4nm on midZ. I
find this encouraging since only 300 positions were used for the calibration and using
400 or 500 targets has been shown to give better calibrations.

I also did a quick check using only 200 targets for each calibration, with the same true
positions of the launch heads as above, and obtained the following results:




            midX         midY         midZ         alpha        beta          gamma
mean of     2.89E-7      -9.30E-7     1.16E-7      -5.01E-7     5.63E-7       1.23E-4
mean
standard    2.37E-6      3.90E-6      6.75E-7      3.40E-6      4.22E-6       2.14E-4
deviation
of mean
mean of     5.69E-7      6.29E-7      1.58E-8      2.00E-5      2.08E-5       1.88E-3
standard
deviation
standard    8.00E-7      1.20E-6      2.20E-8      2.71E-5      3.05E-5       1.88E-3
deviation
of
standard
deviation

Figures are saved as means2 and SD2 in calibrationbowtie4.

We can see here the huge difference made by using just 200 plate 2 positions for the
calibration instead of 300.




                                                                                 Page 40
                                                                         Patrick Gloster




8) Saving data

Any programs I have written which record data automatically save two types: pictures
of the histograms and text files containing different variables calculated/present in the
function.

Saving figures

      When the figure window is created a figure handle is assigned to it by
       h1=figure('Name',Whatever name you want’), where the figure handle is h1.
       This is defined BEFORE the figure is plotted
      The figure is then plotted by the command hist(vector, number of bins)
      The figure is saved by writing ‘saveas(figure handle, filename)’ e.g.
       ‘saveas(h1,myfigure)’
      To avoid overwriting saved files many of my functions have an input
       argument which is a number assigned to the filename you save the picture as.
       For example in bowtie5 the figures are saved by
       ‘saveas(h1,sprintf('bowtie5_%d_figure1.bmp',andy));’ where sprintf is used to
       read the number represented by andy and put this in the place of the ‘%d’. If
       andy was 10, the histogram would be saved as bowtie5_10_figure1.bmp
      Warning – this method of automatic saving appears to be a fairly poor method,
       since it seems to be some form of print screen and sometimes it takes a picture
       of some other application open at the time, rather than the picture you want
       saved

To access the pictures just open the Matlab folder. If the program you were running is
saved in the Code folder (as it probably will be), the figure will be automatically
saved in this folder. Double click on it to open the figure.


Saving text files

      First we open a file to save the information to with the line fid =
       fopen(filename, 'w'); where the w tells Matlab we want are opening the file
       with the intention of writing information to it
      We then write information using ‘fprintf(fid,’words = %d’, variable name)’.
       As an example, to write the word ‘midX’ next to the value for ‘midX’ we
       would have the line ‘fprintf(fid,’midX = %d’, midX)’
      Use \n for a new line
      To close the file we have the line fclose(fid);

To access the text files which are in a particular folder make sure you have selected
that folder in the current directory at the top centre of the screen then type ‘type
filename.txt’ in the Matlab command window.

If you are running programs saved in the Matlab/Code folder then you will already be
in this directory and all text files will be automatically saved here. If you haven’t
moved them then just typing ‘type filename.txt’ will open them


                                                                                 Page 41
                                                                           Patrick Gloster


If the files are in a different folder and you do not want to change directory type ‘type
foldername/filename.txt’ to open them
9) Current Matlab code

This section contains copies of the 3 most recent programs I have written, which will
be used in the future. To change the range of movement of plate 2 alter platerangeinX,
platerangeinY and platerangeinZ. All distance units are in metres, and all angle units
are in radians.

Just to reiterate, the variables rangeinmidX, rangeinmidY and rangeinmidZ define
how far in either direction from a central point the plates can move along that axis, so
if these are all set to 0.1 then the full range of the midpoints of the plate will be 0.2m3
about a central point.

The manner in which the positions of plate 2 are generated is as follows:

   3) We take a point 1m away from the centre of plate 1 along the Z-axis
   4) We add a random number between -1 and 1 multiplied by the variables
      rangeinmidX/Y/Z to this point to give us the centre of plate 2

CALIBRATIONBOWTIE3.M

function calibrationbowtie3(N,M)

format long

%This sets the seed for the random number by the clock
time=clock;
seed=time(6)+time(5)+time(4)+time(3)+time(2)+time(1);
rand('seed',seed);
randn('seed',seed);

%These are the distances between retroreflectors on plate 2
d1=0.1;
d2=0.1;


initialerror=1E-3; %This is the error on our initial estimates of the
positions, 1mm
lasererror=1E-9; %This is the error on the laser measurements, 1nm
distancebetweenplates=1; %This is the rough distance between the
plates along the Z axis, 1 metre
angleerror=(pi/180); %This is the error on how well we can determine
the angle, currently 1 degree

%These are the ranges of movement along each axis that plate 2 has
platerangeinX=0.01;
platerangeinY=0.01;
platerangeinZ=0.01;

%These are the sizes of the angles that plate 2 can rotate through
alphavariation=2.*pi/9;
betavariation=2.*pi/9;
gammavariation=2.*pi/9;

%Define the points on plate 1

                                                                                   Page 42
                                                           Patrick Gloster


a1=0;
b1=0;
c1=0;

a2=0;
b2=0.02;
c2=0;

a3=0;
b3=0.078;
c3=0;

a4=0;
b4=0.098;
c4=0;

a5=-0.05;
b5=0.045;
c5=0;

a6=-0.03;
b6=0.045;
c6=0;

a7=0.028;
b7=0.05;
c7=0.002;

a8=0.048;
b8=0.05;
c8=0.002;

%Generate the midpoint values of plate 2 for N positions

midx=platerangeinX.*2.*(rand(N,1)-0.5);
midy=(b4/2)+platerangeinY.*2.*(rand(N,1)-0.5);
midz=distancebetweenplates+platerangeinZ.*2.*(rand(N,1)-0.5);

%Generate the angles of orientation of plate 2 for each of N
positions

alpha=2.*(rand(N,1)-0.5).*alphavariation;
beta=2.*(rand(N,1)-0.5).*betavariation;
gamma=2.*(rand(N,1)-0.5).*gammavariation;

%Generate the true positions of points 1, 2, 3 and 4

x1=midx+(d1/2).*cos(alpha).*sin(gamma);
y1=midy-(d1/2).*cos(alpha).*cos(gamma);
z1=midz+(d1/2).*sin(alpha);

x2=midx-(d1/2).*cos(alpha).*sin(gamma);
y2=midy+(d1/2).*cos(alpha).*cos(gamma);
z2=midz-(d1/2).*sin(alpha);

x3=midx+(d2/2).*cos(beta).*cos(gamma);
y3=midy+(d2/2).*cos(beta).*sin(gamma);
z3=midz-(d2/2).*sin(beta);

                                                                  Page 43
                                                           Patrick Gloster


x4=midx-(d2/2).*cos(beta).*cos(gamma);
y4=midy-(d2/2).*cos(beta).*sin(gamma);
z4=midz+(d2/2).*sin(beta);

truevalues=[midx,midy,midz,alpha,beta,gamma];
truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4];


truevalues=truevalues';
truevalues=truevalues(:);

truepositions=truepositions';
truepositions=truepositions(:);




%Generate our starting approximations for a's, b's and c's. Remember,
a1=0,
%b1=0, c1=0, a4=0, c4=0 and c6=0 so aren't included as variables

approx_a2=a2+randn.*initialerror;
approx_b2=b2+randn.*initialerror;
approx_c2=c2+randn.*initialerror;

approx_a3=a3+randn.*initialerror;
approx_b3=b3+randn.*initialerror;
approx_c3=c3+randn.*initialerror;


approx_b4=b4+randn.*initialerror;


approx_a5=a5+randn.*initialerror;
approx_b5=b5+randn.*initialerror;


approx_a6=a6+randn.*initialerror;
approx_b6=b6+randn.*initialerror;
approx_c6=c6+randn.*initialerror;

approx_a7=a7+randn.*initialerror;
approx_b7=b7+randn.*initialerror;
approx_c7=c7+randn.*initialerror;

approx_a8=a8+randn.*initialerror;
approx_b8=b8+randn.*initialerror;
approx_c8=c8+randn.*initialerror;

%Generate our initial approximations to d1 and d2:
approx_d1=d1+randn.*initialerror;
approx_d2=d2+randn.*initialerror;


%Generate our approximations to the midpoints of plate 2

approx_midx=midx+randn(N,1).*initialerror;
approx_midy=midy+randn(N,1).*initialerror;

                                                                  Page 44
                                                           Patrick Gloster

approx_midz=midz+randn(N,1).*initialerror;

%Generate our approximations to the angles of orientation of plate 2
approx_alpha=alpha+(randn(N,1).*angleerror);
approx_beta=beta+(rand(N,1).*angleerror);
approx_gamma=gamma+(rand(N,1).*angleerror);


approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro
x_beta,approx_gamma];
approx_values=approx_values';
approx_values=approx_values(:);

%Calculate the true distances
D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2);
D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2);
D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2);
D4=sqrt((x2).^2+(y2).^2+(z2).^2);
D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2);
D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2);
D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2);
D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2);

%Generate the measurements by adding a random error to each distance

L1=D1+randn(N,1).*lasererror;
L2=D2+randn(N,1).*lasererror;
L3=D3+randn(N,1).*lasererror;
L4=D4+randn(N,1).*lasererror;
L5=D5+randn(N,1).*lasererror;
L6=D6+randn(N,1).*lasererror;
L7=D7+randn(N,1).*lasererror;
L8=D8+randn(N,1).*lasererror;

measurements=[L1,L2,L3,L4,L5,L6,L7,L8];


%Form a vector of our initial estimates of all variables

X0=[approx_a2;approx_b2;approx_c2;approx_a3;approx_b3;approx_c3;appro
x_b4;...

approx_a5;approx_b5;approx_a6;approx_b6;approx_c6;approx_a7;approx_b7
;...

approx_c7;approx_a8;approx_b8;approx_c8;approx_d1;approx_d2;approx_va
lues];

%We need to use an anonymous function as the objective function
because
%superhappyfun has 2 input arguments and lsqnonlin can only minimize
%functions with 1 input argument

anon = @(X) superhappyfun(measurements,X);

%Set lsqnonlin to do one iteration. This will quickly perform the
%minimzation with just one iteration, giving us a result with very
bad
%accuracy. However, what we are interested in is the jacobian which
is

                                                                  Page 45
                                                           Patrick Gloster

%returned. Using this jacobian the sparsity structure is obtained and
then
%handed to lsqnonlin for the next minimization, which has options
suitable
%for our minimization and will run much faster using the sparsity
%structure.

options=optimset('MaxIter',1);
[dummy_X,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,...
    dummy_lambda,dummy_jacobian]=lsqnonlin(anon,X0,[],[],options);

%Derive the sparsity structure

jac=dummy_jacobian;
    for a=1:8.*N
        for b=1:6.*N+20
            if jac(a,b)~=0
                 jac(a,b)=1;
            end
        end
    end


%Set the options for lsqnonlin suitable for a good minimization

options=optimset('TolFun',1E-18,'TolX',1E-
10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac
obPattern',jac);

%Use the minimization function lsqnonlin to calculate values for
%the positions of the launch heads on plate 1. These values are only
%calculated once

[X,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(anon,X
0,[],[],options);

calibrations=X(1:20);

%Clear our previous variables. We can now calculate the positions of
a set
%of plate 2 positions using the calibration constants we have just
%calculated

clear   midx midy midz alpha beta gamma x1 y1 z1 x2 y2 z2 x3 y3 z3 x4
y4 z4
clear   truevalues truepositions approx_midx approx_midy approx_midz
clear   approx_alpha approx_beta approx_gamma approx_values
clear   D1 D2 D3 D4 D5 D6 D7 D8 L1 L2 L3 L4 L5 L6 L7 L8 measurements

%Now generate M plate 2 positions and, using the calibrations we have
just
%worked out, calculate their positions

%Generate the midpoint values of plate 2 for M positions

midx=platerangeinX.*2.*(rand(M,1)-0.5);
midy=(b4/2)+platerangeinY.*2.*(rand(M,1)-0.5);
midz=distancebetweenplates+platerangeinZ.*2.*(rand(M,1)-0.5);



                                                                  Page 46
                                                           Patrick Gloster

%Generate the angles of orientation of plate 2 for each of M
positions

alpha=2.*(rand(M,1)-0.5).*alphavariation;
beta=2.*(rand(M,1)-0.5).*betavariation;
gamma=2.*(rand(M,1)-0.5).*gammavariation;

%Generate the true positions of points 1, 2, 3 and 4

x1=midx+(d1/2).*cos(alpha).*sin(gamma);
y1=midy-(d1/2).*cos(alpha).*cos(gamma);
z1=midz+(d1/2).*sin(alpha);

x2=midx-(d1/2).*cos(alpha).*sin(gamma);
y2=midy+(d1/2).*cos(alpha).*cos(gamma);
z2=midz-(d1/2).*sin(alpha);

x3=midx+(d2/2).*cos(beta).*cos(gamma);
y3=midy+(d2/2).*cos(beta).*sin(gamma);
z3=midz-(d2/2).*sin(beta);

x4=midx-(d2/2).*cos(beta).*cos(gamma);
y4=midy-(d2/2).*cos(beta).*sin(gamma);
z4=midz+(d2/2).*sin(beta);

truevalues=[midx,midy,midz,alpha,beta,gamma];
truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4];


truevalues=truevalues';
truevalues=truevalues(:);

truepositions=truepositions';
truepositions=truepositions(:);

%Generate our approximations to the midpoints of plate 2

approx_midx=midx+randn(M,1).*initialerror;
approx_midy=midy+randn(M,1).*initialerror;
approx_midz=midz+randn(M,1).*initialerror;

%Generate our approximations to the angles of orientation of plate 2
approx_alpha=alpha+(randn(M,1).*angleerror);
approx_beta=beta+(rand(M,1).*angleerror);
approx_gamma=gamma+(rand(M,1).*angleerror);


approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro
x_beta,approx_gamma];
approx_values=approx_values';
approx_values=approx_values(:);

%Calculate the true distances
D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2);
D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2);
D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2);
D4=sqrt((x2).^2+(y2).^2+(z2).^2);
D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2);
D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2);
D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2);
                                                                  Page 47
                                                            Patrick Gloster

D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2);

%Generate the measurements by adding a random error to each distance

L1=D1+randn(M,1).*lasererror;
L2=D2+randn(M,1).*lasererror;
L3=D3+randn(M,1).*lasererror;
L4=D4+randn(M,1).*lasererror;
L5=D5+randn(M,1).*lasererror;
L6=D6+randn(M,1).*lasererror;
L7=D7+randn(M,1).*lasererror;
L8=D8+randn(M,1).*lasererror;

measurements=[L1,L2,L3,L4,L5,L6,L7,L8];

X00=approx_values;

anon2 = @(XX) calibrationfun(measurements,XX,calibrations);

%Obtain the sparsity structure for this next minimization

options=optimset('MaxIter',1);

[dummy_XX,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,..
.
    dummy_lambda,dummy_jacobian2]=lsqnonlin(anon2,X00,[],[],options);

    jac2=dummy_jacobian2;

    for a=1:8.*M
        for b=1:6.*M
            if jac2(a,b)~=0
                 jac2(a,b)=1;
            end
        end
    end


%Calculate the values for the positions of plate 2 using the
following
%options with lsqnonlin

options=optimset('TolFun',1E-18,'TolX',1E-
10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac
obPattern',jac2);

[XX,resnorm,residual,exitflag,output,lambda,jacobian2]=lsqnonlin(anon
2,X00,[],[],options);


%Calculate the differences between the true and calculated values for
midx,
%midy, midz, alpha, beta and gamma
difference=XX-truevalues;
midx_difference=difference(1:6:end);
midy_difference=difference(2:6:end);
midz_difference=difference(3:6:end);
alpha_difference=difference(4:6:end);
beta_difference=difference(5:6:end);
gamma_difference=difference(6:6:end);
                                                                   Page 48
                                                         Patrick Gloster


%Calculate the standard deviation of these differences
midx_SD=std(midx_difference);
midy_SD=std(midy_difference);
midz_SD=std(midz_difference);
alpha_SD=std(alpha_difference);
beta_SD=std(beta_difference);
gamma_SD=std(gamma_difference);

%Calculate the mean of these differences
midx_mean=mean(midx_difference);
midy_mean=mean(midy_difference);
midz_mean=mean(midz_difference);
alpha_mean=mean(alpha_difference);
beta_mean=mean(beta_difference);
gamma_mean=mean(gamma_difference);


%Plot histograms for the results.. Histograms show the differences
between the true and
%calculated values of midx, midy, midz, alpha, beta and gamma

h1=figure('Name','calibration bowtie 3
results','NumberTitle','off','Units'...
     ,'normalized','Position',[0 0 1 0.8]);
subplot(2,3,1);
hist(midx_difference,20)
title('midx difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',midx_mean,midx_SD),'Location','SouthOutside')

subplot(2,3,2);
hist(midy_difference,20)
title('midy difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',midy_mean,midy_SD),'Location','SouthOutside')

subplot(2,3,3);
hist(midz_difference,20)
title('midz difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',midz_mean,midz_SD),'Location','SouthOutside')

subplot(2,3,4);
hist(alpha_difference,20)
title('alpha difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',alpha_mean,alpha_SD),'Location','SouthOutside')

subplot(2,3,5);
hist(beta_difference,20)
title('beta difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',beta_mean,beta_SD),'Location','SouthOutside')

subplot(2,3,6);
hist(gamma_difference,20)
title('gamma difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',gamma_mean,gamma_SD),'Location','SouthOutside')


                                                                Page 49
                                                         Patrick Gloster




CALIBRATIONBOWTIE4.M

function calibrationbowtie4(N,M,repeats)

format long

midx_SD=zeros(repeats,1);
midy_SD=zeros(repeats,1);
midz_SD=zeros(repeats,1);
alpha_SD=zeros(repeats,1);
beta_SD=zeros(repeats,1);
gamma_SD=zeros(repeats,1);

midx_mean=zeros(repeats,1);
midy_mean=zeros(repeats,1);
midz_mean=zeros(repeats,1);
alpha_mean=zeros(repeats,1);
beta_mean=zeros(repeats,1);
gamma_mean=zeros(repeats,1);

%This sets the seed for the random number by the clock
time=clock;
seed=time(6)+time(5)+time(4)+time(3)+time(2)+time(1);
rand('seed',seed);
randn('seed',seed);

%These are the distances between retroreflectors on plate 2
d1=0.1;
d2=0.1;


initialerror=1E-3; %This is the error on our initial estimates of the
positions, 1mm
lasererror=1E-9; %This is the error on the laser measurements, 1nm
distancebetweenplates=1; %This is the rough distance between the
plates along the Z axis, 1 metre
angleerror=(pi/180); %This is the error on how well we can determine
the angle, currently 1 degree

%These are the ranges of movement along each axis that plate 2 has
platerangeinX=0.01;
platerangeinY=0.01;

                                                                Page 50
                                                           Patrick Gloster

platerangeinZ=0.01;

%These are the sizes of the angles that plate 2 can rotate through
alphavariation=2.*pi/9;
betavariation=2.*pi/9;
gammavariation=2.*pi/9;

%Define the points on plate 1

a1=0;
b1=0;
c1=0;

a2=0;
b2=0.02;
c2=0;

a3=0;
b3=0.08;
c3=0;

a4=0;
b4=0.1;
c4=0;

a5=-0.05;
b5=0.045;
c5=0;

a6=-0.03;
b6=0.045;
c6=0;

a7=0.03;
b7=0.05;
c7=0;

a8=0.05;
b8=0.05;
c8=0;

for count=1:repeats

%Generate the midpoint values of plate 2 for N positions

midx=platerangeinX.*2.*(rand(N,1)-0.5);
midy=(b4/2)+platerangeinY.*2.*(rand(N,1)-0.5);
midz=distancebetweenplates+platerangeinZ.*2.*(rand(N,1)-0.5);

%Generate the angles of orientation of plate 2 for each of N
positions

alpha=2.*(rand(N,1)-0.5).*alphavariation;
beta=2.*(rand(N,1)-0.5).*betavariation;
gamma=2.*(rand(N,1)-0.5).*gammavariation;

%Generate the true positions of points 1, 2, 3 and 4

x1=midx+(d1/2).*cos(alpha).*sin(gamma);

                                                                  Page 51
                                                         Patrick Gloster

y1=midy-(d1/2).*cos(alpha).*cos(gamma);
z1=midz+(d1/2).*sin(alpha);

x2=midx-(d1/2).*cos(alpha).*sin(gamma);
y2=midy+(d1/2).*cos(alpha).*cos(gamma);
z2=midz-(d1/2).*sin(alpha);

x3=midx+(d2/2).*cos(beta).*cos(gamma);
y3=midy+(d2/2).*cos(beta).*sin(gamma);
z3=midz-(d2/2).*sin(beta);

x4=midx-(d2/2).*cos(beta).*cos(gamma);
y4=midy-(d2/2).*cos(beta).*sin(gamma);
z4=midz+(d2/2).*sin(beta);

truevalues=[midx,midy,midz,alpha,beta,gamma];
truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4];


truevalues=truevalues';
truevalues=truevalues(:);

truepositions=truepositions';
truepositions=truepositions(:);




%Generate our starting approximations for a's, b's and c's. Remember,
a1=0,
%b1=0, c1=0, a4=0, c4=0 and c6=0 so aren't included as variables

approx_a2=a2+randn.*initialerror;
approx_b2=b2+randn.*initialerror;
approx_c2=c2+randn.*initialerror;

approx_a3=a3+randn.*initialerror;
approx_b3=b3+randn.*initialerror;
approx_c3=c3+randn.*initialerror;


approx_b4=b4+randn.*initialerror;


approx_a5=a5+randn.*initialerror;
approx_b5=b5+randn.*initialerror;


approx_a6=a6+randn.*initialerror;
approx_b6=b6+randn.*initialerror;
approx_c6=c6+randn.*initialerror;

approx_a7=a7+randn.*initialerror;
approx_b7=b7+randn.*initialerror;
approx_c7=c7+randn.*initialerror;

approx_a8=a8+randn.*initialerror;
approx_b8=b8+randn.*initialerror;
approx_c8=c8+randn.*initialerror;


                                                                Page 52
                                                           Patrick Gloster

%Generate our initial approximations to d1 and d2:
approx_d1=d1+randn.*initialerror;
approx_d2=d2+randn.*initialerror;


%Generate our approximations to the midpoints of plate 2

approx_midx=midx+randn(N,1).*initialerror;
approx_midy=midy+randn(N,1).*initialerror;
approx_midz=midz+randn(N,1).*initialerror;

%Generate our approximations to the angles of orientation of plate 2
approx_alpha=alpha+(randn(N,1).*angleerror);
approx_beta=beta+(rand(N,1).*angleerror);
approx_gamma=gamma+(rand(N,1).*angleerror);


approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro
x_beta,approx_gamma];
approx_values=approx_values';
approx_values=approx_values(:);

%Calculate the true distances
D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2);
D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2);
D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2);
D4=sqrt((x2).^2+(y2).^2+(z2).^2);
D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2);
D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2);
D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2);
D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2);

%Generate the measurements by adding a random error to each distance

L1=D1+randn(N,1).*lasererror;
L2=D2+randn(N,1).*lasererror;
L3=D3+randn(N,1).*lasererror;
L4=D4+randn(N,1).*lasererror;
L5=D5+randn(N,1).*lasererror;
L6=D6+randn(N,1).*lasererror;
L7=D7+randn(N,1).*lasererror;
L8=D8+randn(N,1).*lasererror;

measurements=[L1,L2,L3,L4,L5,L6,L7,L8];


%Form a vector of our initial estimates of all variables

X0=[approx_a2;approx_b2;approx_c2;approx_a3;approx_b3;approx_c3;appro
x_b4;...

approx_a5;approx_b5;approx_a6;approx_b6;approx_c6;approx_a7;approx_b7
;...

approx_c7;approx_a8;approx_b8;approx_c8;approx_d1;approx_d2;approx_va
lues];

%We need to use an anonymous function as the objective function
because
%superhappyfun has 2 input arguments and lsqnonlin can only minimize

                                                                  Page 53
                                                           Patrick Gloster

%functions with 1 input argument

anon = @(X) superhappyfun(measurements,X);

%Set lsqnonlin to do one iteration. This will quickly perform the
%minimzation with just one iteration, giving us a result with very
bad
%accuracy. However, what we are interested in is the jacobian which
is
%returned. Using this jacobian the sparsity structure is obtained and
then
%handed to lsqnonlin for the next minimization, which has options
suitable
%for our minimization and will run much faster using the sparsity
%structure.

options=optimset('MaxIter',1);
[dummy_X,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,...
    dummy_lambda,dummy_jacobian]=lsqnonlin(anon,X0,[],[],options);

%Derive the sparsity structure

jac=dummy_jacobian;
    for a=1:8.*N
        for b=1:6.*N+20
            if jac(a,b)~=0
                 jac(a,b)=1;
            end
        end
    end


%Set the options for lsqnonlin suitable for a good minimization

options=optimset('TolFun',1E-18,'TolX',1E-
10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac
obPattern',jac);

%Use the minimization function lsqnonlin to calculate values for
%the positions of the launch heads on plate 1. These values are only
%calculated once

[X,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(anon,X
0,[],[],options);

calibrations=X(1:20);


clear   midx midy midz alpha beta gamma x1 y1 z1 x2 y2 z2 x3 y3 z3 x4
y4 z4
clear   truevalues truepositions approx_midx approx_midy approx_midz
clear   approx_alpha approx_beta approx_gamma approx_values
clear   D1 D2 D3 D4 D5 D6 D7 D8 L1 L2 L3 L4 L5 L6 L7 L8 measurements

%Now generate M plate 2 positions and, using the calibrations we have
just
%worked out, calculate their positions

%Generate the midpoint values of plate 2 for M positions


                                                                  Page 54
                                                           Patrick Gloster

midx=platerangeinX.*2.*(rand(M,1)-0.5);
midy=(b4/2)+platerangeinY.*2.*(rand(M,1)-0.5);
midz=distancebetweenplates+platerangeinZ.*2.*(rand(M,1)-0.5);

%Generate the angles of orientation of plate 2 for each of M
positions

alpha=2.*(rand(M,1)-0.5).*alphavariation;
beta=2.*(rand(M,1)-0.5).*betavariation;
gamma=2.*(rand(M,1)-0.5).*gammavariation;

%Generate the true positions of points 1, 2, 3 and 4

x1=midx+(d1/2).*cos(alpha).*sin(gamma);
y1=midy-(d1/2).*cos(alpha).*cos(gamma);
z1=midz+(d1/2).*sin(alpha);

x2=midx-(d1/2).*cos(alpha).*sin(gamma);
y2=midy+(d1/2).*cos(alpha).*cos(gamma);
z2=midz-(d1/2).*sin(alpha);

x3=midx+(d2/2).*cos(beta).*cos(gamma);
y3=midy+(d2/2).*cos(beta).*sin(gamma);
z3=midz-(d2/2).*sin(beta);

x4=midx-(d2/2).*cos(beta).*cos(gamma);
y4=midy-(d2/2).*cos(beta).*sin(gamma);
z4=midz+(d2/2).*sin(beta);

truevalues=[midx,midy,midz,alpha,beta,gamma];
truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4];


truevalues=truevalues';
truevalues=truevalues(:);

truepositions=truepositions';
truepositions=truepositions(:);

%Generate our approximations to the midpoints of plate 2

approx_midx=midx+randn(M,1).*initialerror;
approx_midy=midy+randn(M,1).*initialerror;
approx_midz=midz+randn(M,1).*initialerror;

%Generate our approximations to the angles of orientation of plate 2
approx_alpha=alpha+(randn(M,1).*angleerror);
approx_beta=beta+(rand(M,1).*angleerror);
approx_gamma=gamma+(rand(M,1).*angleerror);


approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro
x_beta,approx_gamma];
approx_values=approx_values';
approx_values=approx_values(:);

%Calculate the true distances
D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2);
D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2);

                                                                  Page 55
                                                            Patrick Gloster

D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2);
D4=sqrt((x2).^2+(y2).^2+(z2).^2);
D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2);
D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2);
D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2);
D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2);

%Generate the measurements by adding a random error to each distance

L1=D1+randn(M,1).*lasererror;
L2=D2+randn(M,1).*lasererror;
L3=D3+randn(M,1).*lasererror;
L4=D4+randn(M,1).*lasererror;
L5=D5+randn(M,1).*lasererror;
L6=D6+randn(M,1).*lasererror;
L7=D7+randn(M,1).*lasererror;
L8=D8+randn(M,1).*lasererror;

measurements=[L1,L2,L3,L4,L5,L6,L7,L8];

X00=approx_values;

anon2 = @(XX) calibrationfun(measurements,XX,calibrations);

%Obtain the sparsity structure for this next minimization

options=optimset('MaxIter',1);

[dummy_XX,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,..
.
    dummy_lambda,dummy_jacobian2]=lsqnonlin(anon2,X00,[],[],options);

    jac2=dummy_jacobian2;

    for a=1:8.*M
        for b=1:6.*M
            if jac2(a,b)~=0
                 jac2(a,b)=1;
            end
        end
    end


%Calculate the values for the positions of plate 2 using the
following
%options with lsqnonlin

options=optimset('TolFun',1E-18,'TolX',1E-
10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac
obPattern',jac2);

[XX,resnorm,residual,exitflag,output,lambda,jacobian2]=lsqnonlin(anon
2,X00,[],[],options);


%Calculate the differences between the true and calculated values for
midx,
%midy, midz, alpha, beta and gamma
difference=XX-truevalues;
midx_difference=difference(1:6:end);
                                                                   Page 56
                                                         Patrick Gloster

midy_difference=difference(2:6:end);
midz_difference=difference(3:6:end);
alpha_difference=difference(4:6:end);
beta_difference=difference(5:6:end);
gamma_difference=difference(6:6:end);

%Calculate the standard deviation of these differences
midx_SD(count)=std(midx_difference);
midy_SD(count)=std(midy_difference);
midz_SD(count)=std(midz_difference);
alpha_SD(count)=std(alpha_difference);
beta_SD(count)=std(beta_difference);
gamma_SD(count)=std(gamma_difference);

%Calculate the mean of these differences
midx_mean(count)=mean(midx_difference);
midy_mean(count)=mean(midy_difference);
midz_mean(count)=mean(midz_difference);
alpha_mean(count)=mean(alpha_difference);
beta_mean(count)=mean(beta_difference);
gamma_mean(count)=mean(gamma_difference);

calibration_numer=count

end

mean_midx_mean=mean(midx_mean);
mean_midy_mean=mean(midy_mean);
mean_midz_mean=mean(midz_mean);
mean_alpha_mean=mean(alpha_mean);
mean_beta_mean=mean(beta_mean);
mean_gamma_mean=mean(gamma_mean);

SD_midx_mean=std(midx_mean);
SD_midy_mean=std(midy_mean);
SD_midz_mean=std(midz_mean);
SD_alpha_mean=std(alpha_mean);
SD_beta_mean=std(beta_mean);
SD_gamma_mean=std(gamma_mean);

SD_midx_SD=std(midx_SD);
SD_midy_SD=std(midy_SD);
SD_midz_SD=std(midz_SD);
SD_alpha_SD=std(alpha_SD);
SD_beta_SD=std(beta_SD);
SD_gamma_SD=std(gamma_SD);

mean_midx_SD=mean(midx_SD);
mean_midy_SD=mean(midy_SD);
mean_midz_SD=mean(midz_SD);
mean_alpha_SD=mean(alpha_SD);
mean_beta_SD=mean(beta_SD);
mean_gamma_SD=mean(gamma_SD);




%Plot histograms for the results.. Histograms show the differences
between the true and
%calculated values of midx, midy, midz, alpha, beta and gamma

                                                                Page 57
                                                         Patrick Gloster


h1=figure('Name','means','NumberTitle','off','Units'...
     ,'normalized','Position',[0 0 1 0.8]);
subplot(2,3,1);
hist(midx_mean,20)
title('midx mean');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_midx_mean,SD_midx_mean),'Location','SouthOutside')

subplot(2,3,2);
hist(midy_mean,20)
title('midy mean');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_midy_mean,SD_midy_mean),'Location','SouthOutside')

subplot(2,3,3);
hist(midz_mean,20)
title('midz mean');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_midz_mean,SD_midz_mean),'Location','SouthOutside')

subplot(2,3,4);
hist(alpha_mean,20)
title('alpha mean');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_alpha_mean,SD_alpha_mean),'Location','SouthOutside
')

subplot(2,3,5);
hist(beta_mean,20)
title('beta mean');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_beta_mean,SD_beta_mean),'Location','SouthOutside')

subplot(2,3,6);
hist(gamma_mean,20)
title('gamma mean');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_gamma_mean,SD_gamma_mean),'Location','SouthOutside
')


h2=figure('Name','Standard deviations','NumberTitle','off','Units'...
     ,'normalized','Position',[0 0 1 0.8]);
subplot(2,3,1);
hist(midx_SD,20)
title('midx SD');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_midx_SD,SD_midx_SD),'Location','SouthOutside')

subplot(2,3,2);
hist(midy_SD,20)
title('midy SD');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_midy_SD,SD_midy_SD),'Location','SouthOutside')

subplot(2,3,3);
hist(midz_SD,20)
title('midz SD');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_midz_SD,SD_midz_SD),'Location','SouthOutside')

                                                                Page 58
                                                         Patrick Gloster


subplot(2,3,4);
hist(alpha_SD,20)
title('alpha SD');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_alpha_SD,SD_alpha_SD),'Location','SouthOutside')

subplot(2,3,5);
hist(beta_SD,20)
title('beta SD');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_beta_SD,SD_beta_SD),'Location','SouthOutside')

subplot(2,3,6);
hist(gamma_SD,20)
title('gamma SD');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_gamma_SD,SD_gamma_SD),'Location','SouthOutside')




                                                                Page 59
                                                               Patrick Gloster


SUPERHAPPYFUN.M – The objective function for calibrationbowtie3 and 4, and
calibrations

function out=superhappyfun(measurements,X)

%This is the objective function for threeDbowtie.

%First we define the current values of the a's, b's and c's on plate
1 from
%the vector X:

a2=X(1);
b2=X(2);
c2=X(3);
a3=X(4);
b3=X(5);
c3=X(6);
b4=X(7);
a5=X(8);
b5=X(9);
a6=X(10);
b6=X(11);
c6=X(12);
a7=X(13);
b7=X(14);
c7=X(15);
a8=X(16);
b8=X(17);
c8=X(18);

%Define the current values of d1 and d2 from the vector X:

d1=X(19);
d2=X(20);

%Define the current values of midx, midy, midz, alpha, beta and gamma
from
%the vector X:

midx=X(21:6:end);
midy=X(22:6:end);
midz=X(23:6:end);
alpha=X(24:6:end);
beta=X(25:6:end);
gamma=X(26:6:end);

%Calculate the positions of points 1, 2, 3 and 4 based on our current
%values for midx, midy, midz, alpha, beta and gamma:

x1=midx+(d1/2).*cos(alpha).*sin(gamma);
y1=midy-(d1/2).*cos(alpha).*cos(gamma);
z1=midz+(d1/2).*sin(alpha);

x2=midx-(d1/2).*cos(alpha).*sin(gamma);
y2=midy+(d1/2).*cos(alpha).*cos(gamma);
z2=midz-(d1/2).*sin(alpha);

x3=midx+(d2/2).*cos(beta).*cos(gamma);
y3=midy+(d2/2).*cos(beta).*sin(gamma);
z3=midz-(d2/2).*sin(beta);

                                                                      Page 60
                                                         Patrick Gloster


x4=midx-(d2/2).*cos(beta).*cos(gamma);
y4=midy-(d2/2).*cos(beta).*sin(gamma);
z4=midz+(d2/2).*sin(beta);

%Calculate what measurements we would observe if the 4 points were at
these
%calculated positions. These are the 'predicted measurements':

M1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2);
M2=sqrt((x1).^2+(y1-b4).^2+(z1).^2);
M3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2);
M4=sqrt((x2).^2+(y2).^2+(z2).^2);
M5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2);
M6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2);
M7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2);
M8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2);


predicted=[M1,M2,M3,M4,M5,M6,M7,M8];

%Calculate the difference between these 'predicted' measurements and
the
%true measurements:

out=predicted-measurements;




                                                                Page 61
                                                              Patrick Gloster


CALIBRATIONFUN.M – A second objective function used by calbration3 and 4
after they have calibrated plate 1

function out=calibrationfun(measurements,X,calibrations)

%This is the objective function for threeDbowtie.

%First we define the current values of the a's, b's and c's on plate
1 from
%the vector calibrations:

a2=calibrations(1);
b2=calibrations(2);
c2=calibrations(3);
a3=calibrations(4);
b3=calibrations(5);
c3=calibrations(6);
b4=calibrations(7);
a5=calibrations(8);
b5=calibrations(9);
a6=calibrations(10);
b6=calibrations(11);
c6=calibrations(12);
a7=calibrations(13);
b7=calibrations(14);
c7=calibrations(15);
a8=calibrations(16);
b8=calibrations(17);
c8=calibrations(18);

%Define the current values of d1 and d2 from the vector calibrations:

d1=calibrations(19);
d2=calibrations(20);

%Define the current values of midx, midy, midz, alpha, beta and gamma
from
%the vector X:

midx=X(1:6:end);
midy=X(2:6:end);
midz=X(3:6:end);
alpha=X(4:6:end);
beta=X(5:6:end);
gamma=X(6:6:end);

%Calculate the positions of points 1, 2, 3 and 4 based on our current
%values for midx, midy, midz, alpha, beta and gamma:

x1=midx+(d1/2).*cos(alpha).*sin(gamma);
y1=midy-(d1/2).*cos(alpha).*cos(gamma);
z1=midz+(d1/2).*sin(alpha);

x2=midx-(d1/2).*cos(alpha).*sin(gamma);
y2=midy+(d1/2).*cos(alpha).*cos(gamma);
z2=midz-(d1/2).*sin(alpha);

x3=midx+(d2/2).*cos(beta).*cos(gamma);
y3=midy+(d2/2).*cos(beta).*sin(gamma);
z3=midz-(d2/2).*sin(beta);

                                                                     Page 62
                                                         Patrick Gloster


x4=midx-(d2/2).*cos(beta).*cos(gamma);
y4=midy-(d2/2).*cos(beta).*sin(gamma);
z4=midz+(d2/2).*sin(beta);

%Calculate what measurements we would observe if the 4 points were at
these
%calculated positions. These are the 'predicted measurements':

M1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2);
M2=sqrt((x1).^2+(y1-b4).^2+(z1).^2);
M3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2);
M4=sqrt((x2).^2+(y2).^2+(z2).^2);
M5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2);
M6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2);
M7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2);
M8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2);


predicted=[M1,M2,M3,M4,M5,M6,M7,M8];

%Calculate the difference between these 'predicted' measurements and
the
%true measurements:

out=predicted-measurements;




                                                                Page 63
                                                         Patrick Gloster


CALIBRATIONS.M

function calibrations(N)

format long

%This sets the seed for the random number by the clock
time=clock;
seed=time(6)+time(5)+time(4)+time(3)+time(2)+time(1);
rand('seed',seed);
randn('seed',seed);

%These are the distances between retroreflectors on plate 2
d1=0.1;
d2=0.1;


initialerror=1E-3; %This is the error on our initial estimates of the
positions, 1mm
lasererror=1E-9; %This is the error on the laser measurements, 1nm
distancebetweenplates=1; %This is the rough distance between the
plates along the Z axis, 1 metre
angleerror=(pi/180); %This is the error on how well we can determine
the angle, currently 1 degree

%These are the ranges of movement along each axis that plate 2 has
platerangeinX=0.01;
platerangeinY=0.01;
platerangeinZ=0.01;

%These are the sizes of the angles that plate 2 can rotate through
alphavariation=2.*pi/9;
betavariation=2.*pi/9;
gammavariation=2.*pi/9;

%Define the points on plate 1

a1=0;
b1=0;
c1=0;

a2=0;
b2=0.02;
c2=0;

a3=0;
b3=0.08;
c3=0;

a4=0;
b4=0.1;
c4=0;

a5=-0.05;
b5=0.045;
c5=0;

a6=-0.03;
b6=0.045;

                                                                Page 64
                                                           Patrick Gloster

c6=0;

a7=0.03;
b7=0.05;
c7=0;

a8=0.05;
b8=0.05;
c8=0;

for count=1:50

%Generate the midpoint values of plate 2 for N positions

midx=platerangeinX.*2.*(rand(N,1)-0.5);
midy=(b4/2)+platerangeinY.*2.*(rand(N,1)-0.5);
midz=distancebetweenplates+platerangeinZ.*2.*(rand(N,1)-0.5);

%Generate the angles of orientation of plate 2 for each of N
positions

alpha=2.*(rand(N,1)-0.5).*alphavariation;
beta=2.*(rand(N,1)-0.5).*betavariation;
gamma=2.*(rand(N,1)-0.5).*gammavariation;

%Generate the true positions of points 1, 2, 3 and 4

x1=midx+(d1/2).*cos(alpha).*sin(gamma);
y1=midy-(d1/2).*cos(alpha).*cos(gamma);
z1=midz+(d1/2).*sin(alpha);

x2=midx-(d1/2).*cos(alpha).*sin(gamma);
y2=midy+(d1/2).*cos(alpha).*cos(gamma);
z2=midz-(d1/2).*sin(alpha);

x3=midx+(d2/2).*cos(beta).*cos(gamma);
y3=midy+(d2/2).*cos(beta).*sin(gamma);
z3=midz-(d2/2).*sin(beta);

x4=midx-(d2/2).*cos(beta).*cos(gamma);
y4=midy-(d2/2).*cos(beta).*sin(gamma);
z4=midz+(d2/2).*sin(beta);

truevalues=[midx,midy,midz,alpha,beta,gamma];
truepositions=[x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4];


truevalues=truevalues';
truevalues=truevalues(:);

truepositions=truepositions';
truepositions=truepositions(:);




%Generate our starting approximations for a's, b's and c's. Remember,
a1=0,
%b1=0, c1=0, a4=0, c4=0 and c6=0 so aren't included as variables


                                                                  Page 65
                                                           Patrick Gloster

approx_a2=a2+randn.*initialerror;
approx_b2=b2+randn.*initialerror;
approx_c2=c2+randn.*initialerror;

approx_a3=a3+randn.*initialerror;
approx_b3=b3+randn.*initialerror;
approx_c3=c3+randn.*initialerror;


approx_b4=b4+randn.*initialerror;


approx_a5=a5+randn.*initialerror;
approx_b5=b5+randn.*initialerror;


approx_a6=a6+randn.*initialerror;
approx_b6=b6+randn.*initialerror;
approx_c6=c6+randn.*initialerror;

approx_a7=a7+randn.*initialerror;
approx_b7=b7+randn.*initialerror;
approx_c7=c7+randn.*initialerror;

approx_a8=a8+randn.*initialerror;
approx_b8=b8+randn.*initialerror;
approx_c8=c8+randn.*initialerror;

%Generate our initial approximations to d1 and d2:
approx_d1=d1+randn.*initialerror;
approx_d2=d2+randn.*initialerror;


%Generate our approximations to the midpoints of plate 2

approx_midx=midx+randn(N,1).*initialerror;
approx_midy=midy+randn(N,1).*initialerror;
approx_midz=midz+randn(N,1).*initialerror;

%Generate our approximations to the angles of orientation of plate 2
approx_alpha=alpha+(randn(N,1).*angleerror);
approx_beta=beta+(rand(N,1).*angleerror);
approx_gamma=gamma+(rand(N,1).*angleerror);


approx_values=[approx_midx,approx_midy,approx_midz,approx_alpha,appro
x_beta,approx_gamma];
approx_values=approx_values';
approx_values=approx_values(:);

%Calculate the true distances
D1=sqrt((x1-a2).^2+(y1-b2).^2+(z1-c2).^2);
D2=sqrt((x1).^2+(y1-b4).^2+(z1).^2);
D3=sqrt((x2-a3).^2+(y2-b3).^2+(z2-c3).^2);
D4=sqrt((x2).^2+(y2).^2+(z2).^2);
D5=sqrt((x3-a7).^2+(y3-b7).^2+(z3-c7).^2);
D6=sqrt((x3-a5).^2+(y3-b5).^2+(z3).^2);
D7=sqrt((x4-a6).^2+(y4-b6).^2+(z4-c6).^2);
D8=sqrt((x4-a8).^2+(y4-b8).^2+(z4-c8).^2);

                                                                  Page 66
                                                           Patrick Gloster


%Generate the measurements by adding a random error to each distance

L1=D1+randn(N,1).*lasererror;
L2=D2+randn(N,1).*lasererror;
L3=D3+randn(N,1).*lasererror;
L4=D4+randn(N,1).*lasererror;
L5=D5+randn(N,1).*lasererror;
L6=D6+randn(N,1).*lasererror;
L7=D7+randn(N,1).*lasererror;
L8=D8+randn(N,1).*lasererror;

measurements=[L1,L2,L3,L4,L5,L6,L7,L8];


%Form a vector of our initial estimates of all variables

X0=[approx_a2;approx_b2;approx_c2;approx_a3;approx_b3;approx_c3;appro
x_b4;...

approx_a5;approx_b5;approx_a6;approx_b6;approx_c6;approx_a7;approx_b7
;...

approx_c7;approx_a8;approx_b8;approx_c8;approx_d1;approx_d2;approx_va
lues];

%We need to use an anonymous function as the objective function
because
%superhappyfun has 2 input arguments and lsqnonlin can only minimize
%functions with 1 input argument

anon = @(X) superhappyfun(measurements,X);

%Set lsqnonlin to do one iteration. This will quickly perform the
%minimzation with just one iteration, giving us a result with very
bad
%accuracy. However, what we are interested in is the jacobian which
is
%returned. Using this jacobian the sparsity structure is obtained and
then
%handed to lsqnonlin for the next minimization, which has options
suitable
%for our minimization and will run much faster using the sparsity
%structure.

options=optimset('MaxIter',1);
[dummy_X,dummy_resnorm,dummy_residual,dummy_exitflag,dummy_output,...
    dummy_lambda,dummy_jacobian]=lsqnonlin(anon,X0,[],[],options);

%Derive the sparsity structure

jac=dummy_jacobian;
    for a=1:8.*N
        for b=1:6.*N+20
            if jac(a,b)~=0
                 jac(a,b)=1;
            end
        end
    end


                                                                  Page 67
                                                         Patrick Gloster


%Set the options for lsqnonlin suitable for a good minimization

options=optimset('TolFun',1E-18,'TolX',1E-
10,'MaxFunEvals',1E6,'MaxIter',1E6,'PlotFcns',@optimplotresnorm2,'Jac
obPattern',jac);

%Use the minimization function lsqnonlin to calculate values for
%the positions of the launch heads on plate 1. These values are only
%calculated once

[X,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(anon,X
0,[],[],options);

calibrations=X(1:20);

a2_difference(count)=a2-calibrations(1);
b2_difference(count)=b2-calibrations(2);
c2_difference(count)=c2-calibrations(3);
a3_difference(count)=a3-calibrations(4);
b3_difference(count)=b3-calibrations(5);
c3_difference(count)=c3-calibrations(6);
b4_difference(count)=b4-calibrations(7);
a5_difference(count)=a5-calibrations(8);
b5_difference(count)=b5-calibrations(9);
a6_difference(count)=a6-calibrations(10);
b6_difference(count)=b6-calibrations(11);
c6_difference(count)=c6-calibrations(12);
a7_difference(count)=a7-calibrations(13);
b7_difference(count)=b7-calibrations(14);
c7_difference(count)=c7-calibrations(15);
a8_difference(count)=a8-calibrations(16);
b8_difference(count)=b8-calibrations(17);
c8_difference(count)=c8-calibrations(18);
d1_difference(count)=d1-calibrations(19);
d2_difference(count)=d2-calibrations(20);

repeat=count

end

mean_a2=mean(a2_difference);
SD_a2=std(a2_difference);
mean_b2=mean(b2_difference);
SD_b2=std(b2_difference);
mean_c2=mean(c2_difference);
SD_c2=std(c2_difference);
mean_a3=mean(a3_difference);
SD_a3=std(a3_difference);
mean_b3=mean(b3_difference);
SD_b3=std(b3_difference);
mean_c3=mean(c3_difference);
SD_c3=std(c3_difference);
mean_b4=mean(b4_difference);
SD_b4=std(b4_difference);
mean_a5=mean(a5_difference);
SD_a5=std(a5_difference);
mean_b5=mean(b5_difference);
SD_b5=std(b5_difference);
mean_a6=mean(a6_difference);
SD_a6=std(a6_difference);

                                                                Page 68
                                                          Patrick Gloster

mean_b6=mean(b6_difference);
SD_b6=std(b6_difference);
mean_c6=mean(c6_difference);
SD_c6=std(c6_difference);
mean_a7=mean(a7_difference);
SD_a7=std(a7_difference);
mean_b7=mean(b7_difference);
SD_b7=std(b7_difference);
mean_c7=mean(c7_difference);
SD_c7=std(c7_difference);
mean_a8=mean(a8_difference);
SD_a8=std(a8_difference);
mean_b8=mean(b8_difference);
SD_b8=std(b8_difference);
mean_c8=mean(c8_difference);
SD_c8=std(c8_difference);
mean_d1=mean(d1_difference);
SD_d1=std(d1_difference);
mean_d2=mean(d2_difference);
SD_d2=std(d2_difference);




h2=figure('Name','calibrations - a
differences','NumberTitle','off','Units'...
     ,'normalized','Position',[0 0 1 0.8]);
subplot(2,3,1);
hist(a2_difference,5)
title('a2 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_a2,SD_a2),'Location','SouthOutside')

subplot(2,3,2);
hist(a3_difference,5)
title('a3 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_a3,SD_a3),'Location','SouthOutside')

subplot(2,3,3);
hist(a5_difference,5)
title('a5 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_a5,SD_a5),'Location','SouthOutside')

subplot(2,3,4);
hist(a6_difference,5)
title('a6 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_a6,SD_a6),'Location','SouthOutside')

subplot(2,3,5);
hist(a7_difference,5)
title('a7 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_a7,SD_a7),'Location','SouthOutside')

subplot(2,3,6);
hist(a8_difference,5)
title('a8 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_a8,SD_a8),'Location','SouthOutside')

                                                                 Page 69
                                                          Patrick Gloster




%Plot histograms of the differences between the true and calculated
values
%of b2, b3, b4, b5, b6, b7 and b8

h3=figure('Name','calibrations - b
differences','NumberTitle','off','Units'...
     ,'normalized','Position',[0 0 1 0.8]);
subplot(2,4,1);
hist(b2_difference,5)
title('b2 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_b2,SD_b2),'Location','SouthOutside')

subplot(2,4,2);
hist(b3_difference,5)
title('b3 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_b3,SD_b3),'Location','SouthOutside')

subplot(2,4,3);
hist(b4_difference,5)
title('b4 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_b4,SD_b4),'Location','SouthOutside')

subplot(2,4,4);
hist(b5_difference,5)
title('b5 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_b5,SD_b5),'Location','SouthOutside')

subplot(2,4,5);
hist(b6_difference,5)
title('b6 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_b6,SD_b6),'Location','SouthOutside')

subplot(2,4,6);
hist(b7_difference,5)
title('b7 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_b7,SD_b7),'Location','SouthOutside')

subplot(2,4,7);
hist(b8_difference,5)
title('b8 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_b8,SD_b8),'Location','SouthOutside')


%Plot histograms of the differences between the true and calculated
values
%of c2, c3, c6, c7 and c8

h4=figure('Name','calibrations - c
differences','NumberTitle','off','Units'...
     ,'normalized','Position',[0 0 1 0.8]);
subplot(2,3,1);

                                                                 Page 70
                                                          Patrick Gloster

hist(c2_difference,5)
title('c2 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_c2,SD_c2),'Location','SouthOutside')

subplot(2,3,2);
hist(c3_difference,5)
title('c3 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_c3,SD_c3),'Location','SouthOutside')

subplot(2,3,3);
hist(c6_difference,5)
title('c6 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_c6,SD_c6),'Location','SouthOutside')

subplot(2,3,4);
hist(c7_difference,5)
title('c7 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_c7,SD_c7),'Location','SouthOutside')

subplot(2,3,5);
hist(c8_difference,5)
title('c8 difference');
legend(sprintf('Mean=%e\nStandard
deviation=%e',mean_c8,SD_c8),'Location','SouthOutside')




                                                                 Page 71

				
DOCUMENT INFO