function [RadiusVCMGIP, DistanceInputPupil, DistanceExitPupil, DistanceExitPupilFixedVCM] = ComputeVCMParametersATBC(TelescopeStation, ...
                                                                                   DL, ...
                                                                                   LabInputChannel, ...
                                                                                   DLPos, ...
                                                                                   longPupilOffset, ...
                                                                                   RadiusFixedVCM)
                                                                                   
% A function to compute the DL-VCM radius of curvature for a given telescope (UT or AT) on a given
% station and a given DL at a certain position.
%
% Usage:
% ~~~~~~
%
% [RadiusVCM, DistanceInputPupil, DistanceExitPupil] = ComputeVCMParameters3(STATION, DL, LAB_CHANNEL, DL_POS, longPupilOffset, RadiusFixedVCM);
%
% with the following input parameters:
%    STATION     == station of the telescope, e.g. 'UT3', 'A0', 'J6', 'M0'
%    DL          == index of the DL used, between 1 ... 8
%    LAB_CHANNEL == index of the laboratory input channel, between 1 ... 8
%    DL_POS      == DL position in units of [m], between 0 ... 60 m
%    longPupilOffset == offset _after_ the BC (a multiplicative factor 4.44^2 applies)
%    RADIUS_FIXED_VCM == radius of curvature of the fixed VCM [m]
%
% and the output parameters:
%    RadiusVCM          == radius of curvature of the VCM required to re-image the pupil
%    DistanceInputPupil == distance between the relayed telescope exit pupil (which
%                          acts as input pupil for the DL) and the DL entrance aperture
%    DistanceExitPupil  == distance between the DL entrance aperture and the re-imaged
%                          exit pupil in the laboratory (behind M16). This does not
%                          take the BCs into account.
%    DistanceExitPupilFixedVCM == distance between the DL entrance aperture and the re-imaged
%                          exit pupil in the laboratory (behind M16) in the case of a VCM with the
%			   fixed radius of curvature. This does not take the BCs into account.
%
% Example:
% ~~~~~~~~
% [RadiusVCM, DistanceInputPupil, DistanceExitPupil, DistanceExitPupilFixedVCM] = ComputeVCMParametersATBC('J6', 3, 1, 17,0, -1);
%
% will compute the radius of the VCM for an AT located at station 'J6', DL 3 @ 17 m, lab input channel 1.
%
% Created from ComputeVCMParameters.m
% GIP
% Date: February 11, 2004
%
% The following information is automatically generated by CVS:
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% $Source: /usr/local/cvsrep/VLTI_Tools/ComputeVCMParameters.m,v $
% $Revision: 1.2 $
% $Date: 2003/07/01 15:06:10 $
% $Author: rwilhelm $
% 
% $Log: ComputeVCMParameters.m,v $
% Revision 1.2  2003/07/01 15:06:10  rwilhelm
% Changed lower limit for VCM radius from 84 mm to 150 mm
%
% Revision 1.1.1.1  2003/07/01 14:52:12  rwilhelm
% Adding new module to the repository
% 
% MODIFIED on November 14, 2005 to correct a bug (incorrect pupil distance
% calculation when BC were by-passed)
%

% Lateral compression factor for the BCs
compFactor = 4.444;
                                                                      
% Get (u, v) of station center
[StationPosition, TelescopeType] = GetStationPosition(TelescopeStation);

% Get u-coordinate (at zero position) and v-coordinate of DL
[DLPosition] = GetDLPosition(DL);

% Compute u-position of DL
if((DL==1)|(DL==2)|(DL==5)|(DL==6))
    DLPosition.u = DLPosition.u0 + DLPos;
else
    DLPosition.u = DLPosition.u0 - DLPos;
end
 
% Compute distance of relayed telescope exit pupil to DL input aperture
% This is now the same for ATs or UTs
    
% Get u-coordinate of DL
UTExitPupiluCoordinate = GetUTExitPupiluCoordinate(DL);
    
% Compute distance between u-position of DL and u-position of telescope relayed exit pupil:
DistanceInputPupil = abs(DLPosition.u - UTExitPupiluCoordinate);
 
% Compute distance between DL input aperture and re-imaged pupil ("behind M16")
% ---------------------------------------------------


% Get re-imaged pupil position (depending on the lab input channel selected and use of BC) 
ReimagedPupilPositionWithBC = GetReimagedPupilPositionWithBC(LabInputChannel);
%ReimagedPupilPositionNoBC = GetReimagedPupilPositionNoBC(LabInputChannel);
inputChannelUcoord = ReimagedPupilPositionWithBC.u;

% ATs (Astro beam) for PRIMA case = going through BC
DistanceExitPupil = abs(DLPosition.u - ReimagedPupilPositionWithBC.u) + ...
                      abs(DLPosition.vOutRef - ReimagedPupilPositionWithBC.v)+longPupilOffset*compFactor*compFactor;
                  

RadiusVCM = ComputeRadiusOfCuravtureVCM(DistanceInputPupil, DistanceExitPupil);
RadiusVCMGIP = ComputeRadiusOfCuravtureVCMGIP(DistanceInputPupil, DistanceExitPupil, DL); %uses optical parameters specific to each DL

DistanceExitPupilFixedVCM = ComputeOutputPupilDistance(DistanceInputPupil, RadiusFixedVCM);


if(strcmp(TelescopeType, 'UT')==1) 

    % UTs (Ref beam)
    % --------------
    
    OffsetPositionExitPupilBeforeBC = (DistanceExitPupilFixedVCM - DistanceExitPupil);
    OffsetPositionExitPupilAfterBC = (DistanceExitPupilFixedVCM - DistanceExitPupil)/(compFactor * compFactor);
else
    % ATs (Astro beam)
    % ----------------
    OffsetPositionExitPupil = (DistanceExitPupilFixedVCM - DistanceExitPupil);
end

fprintf('\n');
fprintf('Telescope = %s, Station = %s, DL %d @ Position %f m, Lab Input Channel #%d, Pupil offset after BCs %.1f m,  VCM radius adjusted \n', ...
        TelescopeType, TelescopeStation, DL, DLPos, LabInputChannel, longPupilOffset);
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
fprintf(' > Input pupil distance (*)                            = %e m\n', DistanceInputPupil);
fprintf(' > Exit pupil distance required                        = %e m\n', DistanceExitPupil);
fprintf(' > VCM ROC required                                 = %e mm\n', RadiusVCM*1000);
fprintf(' > VCM curv required    = %e mm-1\n', 1/(RadiusVCM*1000));
fprintf(' > VCM ROC required (with specific DL optical par.) = %e mm\n', RadiusVCMGIP*1000);
fprintf(' > VCM curv required (with specific DL optical par.) = %e mm-1\n', 1/(RadiusVCMGIP*1000));

if((RadiusVCM < 150.e-03)|(RadiusVCM > 2.8))
    fprintf('   WARNING: VCM radius out of range [150 mm ... 2800 mm]\n');
    fprintf('            150 mm = software limit; Remark: hardware lower limit = 84 mm\n');
    if(abs(RadiusVCM) > 2.8)
        RadiusVCM = 2.8;
    else
        RadiusVCM = 150.e-03;
    end
    fprintf('            It is suggested to use a radius R = %e m\n', RadiusVCM);
end
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
fprintf('  (*) from relayed telescope pupil to DL entrance aperture\n');

fprintf('\n');
fprintf('Telescope = %s, Station = %s, DL %d @ Position %f m, Lab Input Channel #%d, VCM fixed radius: %f m\n', ...
        TelescopeType, TelescopeStation, DL, DLPos, LabInputChannel, RadiusFixedVCM);
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
fprintf(' > Case fixed VCM: Distance from the DL front plate to the exit pupil  = %e m\n', DistanceExitPupilFixedVCM);

if(strcmp(TelescopeType, 'UT')==1) 

    % UTs (Ref beam)
    % --------------

    fprintf(' > Case fixed VCM: Distance of the exit pupil before BC wrt its nominal position  = %e m\n', OffsetPositionExitPupilBeforeBC);
    fprintf(' > Case fixed VCM: Absolute position of the exit pupil before BC: V = %e m\n', OffsetPositionExitPupilBeforeBC+ReimagedPupilPositionWithBC.v);
    fprintf(' > Case fixed VCM: Distance of the exit pupil after BC wrt its nominal position  = %e m\n', OffsetPositionExitPupilAfterBC);
    fprintf(' >                 (Positive offset = pupil farther (optically speaking) than nominally)\n');
else
    % ATs (Astro beam)
    % ----------------
%    fprintf(' > Case fixed VCM: Distance of the exit pupil wrt its nominal position  = %e m\n', OffsetPositionExitPupil);
%    fprintf(' > Case fixed VCM: Absolute position of the exit pupil: V = %e m\n', OffsetPositionExitPupil+ReimagedPupilPositionNoBC.v);
%    fprintf(' >                 (Positive offset = pupil farther (optically speaking) than nominally)\n');
end





function [StationPosition, TelescopeType] = GetStationPosition(TelescopeStation)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Station positions in u v coordinates:                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch TelescopeStation
    case 'UT1'
        StationPosition = struct('u',-16.,'v',-16.);
        TelescopeType = 'UT';
    case 'UT2'
        StationPosition = struct('u',24.,'v',24.);
        TelescopeType = 'UT';
    case 'UT3'
        StationPosition = struct('u',64.,'v',48.);
        TelescopeType = 'UT';
    case 'UT4'
        StationPosition = struct('u',112.,'v', 8.);
        TelescopeType = 'UT';
    case 'A0'
        %A Series (@ u = -32 m)
        StationPosition = struct('u',-32.,'v',-48.);
        TelescopeType = 'AT';
    case 'A1'
        StationPosition = struct('u',-32.,'v',-64.);
        TelescopeType = 'AT';
        %B Series (@ u = -24 m)
    case 'B0'
        StationPosition = struct('u',-24.,'v',-48.);
        TelescopeType = 'AT';
    case 'B1'
        StationPosition = struct('u',-24.,'v',-64.);
        TelescopeType = 'AT';
    case 'B2'
        StationPosition = struct('u',-24.,'v',-72.);
        TelescopeType = 'AT';
    case 'B3'
        StationPosition = struct('u',-24.,'v',-80.);
        TelescopeType = 'AT';
    case 'B4'
        StationPosition = struct('u',-24.,'v',-88.);
        TelescopeType = 'AT';
    case 'B5'
        StationPosition = struct('u',-24.,'v',-96.);
        TelescopeType = 'AT';
    %C Series (@ u = -16 m)
    case 'C0'
        StationPosition = struct('u',-16.,'v',-48.);
        TelescopeType = 'AT';
    case 'C1'
        StationPosition = struct('u',-16.,'v',-64.);
        TelescopeType = 'AT';
    case 'C2'
        StationPosition = struct('u',-16.,'v',-72.);
        TelescopeType = 'AT';
    case 'C3'
        StationPosition = struct('u',-16.,'v',-80.);
        TelescopeType = 'AT';
    %D Series (@ u = 0 m)
    case 'D0'
        StationPosition = struct('u',0.,'v',-48.);
        TelescopeType = 'AT';
    case 'D1'
        StationPosition = struct('u',0.,'v',-80.);
        TelescopeType = 'AT';
    case 'D2'
        StationPosition = struct('u',0.,'v',-96.);
        TelescopeType = 'AT';
    %E Series (@ u = 16 m)
    case 'E0'
        StationPosition = struct('u',16.,'v',-48.);
        TelescopeType = 'AT';
    %G Series (@ u = 32 m)
    case 'G0'
        StationPosition = struct('u',32.,'v',-48.);
        TelescopeType = 'AT';
    case 'G1'
        StationPosition = struct('u',32.,'v',-112.);
        TelescopeType = 'AT';
    case 'G2'
        StationPosition = struct('u',32.,'v',-24.);
        TelescopeType = 'AT';
    %H Series (@ u = 64 m)
    case 'H0'
        StationPosition = struct('u',64.,'v',-48.);
        TelescopeType = 'AT';
    %I Series (@ u = 72 m)
    case 'I1'
        StationPosition = struct('u',72.,'v',-88.);
        TelescopeType = 'AT';
    %J Series (@ u = 88 m)
    case 'J1'
        StationPosition = struct('u',88.,'v',-72.);
        TelescopeType = 'AT';
    case 'J2'
        StationPosition = struct('u',88.,'v',-96.);
        TelescopeType = 'AT';
    case 'J3'
        StationPosition = struct('u',88.,'v',8.);
        TelescopeType = 'AT';
    case 'J4'
        StationPosition = struct('u',88.,'v',24.);
        TelescopeType = 'AT';
    case 'J5'
        StationPosition = struct('u',88.,'v',48.);
        TelescopeType = 'AT';
    case 'J6'
        StationPosition = struct('u',88.,'v',72.);
        TelescopeType = 'AT';
    %K Series (@ u = 96 m)
    case 'K0'
        StationPosition = struct('u',96.,'v',-48.);
        TelescopeType = 'AT';
    %L Series (@ u = 104 m)
    case 'L0'
        StationPosition = struct('u',104.,'v',-48.);
        TelescopeType = 'AT';
    %M Series (@ u = 112 m)
    case 'M0'
        StationPosition = struct('u',112.,'v',-48.);
        TelescopeType = 'AT';
end


function [DLPosition] = GetDLPosition(DL)

switch DL
    case 1
        DLPosition = struct('v', -37.125, 'u0', 59.717);  % Mechanical pos of DL given in table A3 of ISS is U=58.920  (delta=0.797 cm ~d(M1-M2))  
        DLPosition.vInRef    = DLPosition.v - 0.12;       % DLPosition seems to include DL internal path
        DLPosition.vOutRef   = DLPosition.v + 0.12;       %
    case 2
        DLPosition = struct('v', -37.875, 'u0', 59.717);
        DLPosition.vInRef    = DLPosition.v - 0.12;
        DLPosition.vOutRef   = DLPosition.v + 0.12;
    case 3
        DLPosition = struct('v', -38.625, 'u0', 44.283);
        DLPosition.vInRef    = DLPosition.v + 0.12;
        DLPosition.vOutRef   = DLPosition.v - 0.12;
    case 4
        DLPosition = struct('v', -39.375, 'u0', 44.283);% Mechanical pos of DL given in table A3 of ISS is 58.920  (delta=0.797 cm ~d(M1-M2))
        DLPosition.vInRef    = DLPosition.v + 0.12;     % DLPosition seems to include DL internal path
        DLPosition.vOutRef   = DLPosition.v - 0.12;
    case 5
        DLPosition = struct('v', -40.625, 'u0', 59.717);
        DLPosition.vInRef    = DLPosition.v - 0.12;
        DLPosition.vOutRef   = DLPosition.v + 0.12;
    case 6
        DLPosition = struct('v', -41.375, 'u0', 59.717);
        DLPosition.vInRef    = DLPosition.v - 0.12;
        DLPosition.vOutRef   = DLPosition.v + 0.12;
    case 7
        DLPosition = struct('v', -42.125, 'u0', 44.283);
        DLPosition.vInRef    = DLPosition.v + 0.12;
        DLPosition.vOutRef   = DLPosition.v - 0.12;
    case 8
        DLPosition = struct('v', -42.875, 'u0', 44.283);
        DLPosition.vInRef    = DLPosition.v + 0.12;
        DLPosition.vOutRef   = DLPosition.v - 0.12;
end

DLPosition.vInAstro  = DLPosition.vOutRef;
DLPosition.vOutAstro = DLPosition.vInRef;


function [UTExitPupiluCoordinate] = GetUTExitPupiluCoordinate(DL)

switch DL
    case 1
        UTExitPupiluCoordinate = 50.625;  % Nominal value (called UTPupilBeforeDL in ISS)
        %UTExitPupiluCoordinate = -6.4763;     %Test Samuel AT3-STS-A0-DL1 with VCMROC=-730mm
        % UTExitPupiluCoordinate =-22.476      %Test Samuel AT3-STS-A1-DL1 with VCMROC=-730mm
    case 2
        UTExitPupiluCoordinate = 49.875;  % Nominal value with AT-VCM G2=-731mm
        %UTExitPupiluCoordinate = 49.875-32.1487+8.6664; %test Samuel G2/STS with flat VCM
        %UTExitPupiluCoordinate = 49.875-32.1487+8.6664-6.235; %test Samuel G2/STS with flat STS-VCM to get DL2-VCM=492 mm for OPL=14m (7m mechanical)
                   
    case 3
        UTExitPupiluCoordinate = 53.375;
    case 4
        UTExitPupiluCoordinate = 54.125;%nominal value (e.g. VCMROC=-210mm for AT4-J2)\
        %UTExitPupiluCoordinate = 119.125 ;%Test Samuel AT4-STS-J2-DL4 with VCMROC=-1000mm and theoretical distATPupFixVCM=25.5m
        %UTExitPupiluCoordinate = 76.88;%Test Samuel AT4-STS-J2-DL4 with VCMROC=-1000mm +PRIMET adjustment of COM12, distATPupFixVCM=67.77 instead of 25.5m
        %UTExitPupiluCoordinate = 68.08;%Test Samuel AT4-STS-J2-DL4 with VCMROC=-1000mm +PRIMET adjustment of COM15=ACOM1, distATPupFixVCM=76.54 instead of 25.5m

    % Added SLE based on ISS doc VLT-SPE-ESO-15410-2114 Issue 4
    case 5 
        UTExitPupiluCoordinate = 47.125;
    case 6
        UTExitPupiluCoordinate = 46.375 % Nominal value (called UTPupilBeforeDL in ISS p.159)
end

function [du, dv] = Get_dudv_AT(StationPosition)

if StationPosition.v > -40
    % Northern station
    du = -0.0889;
    dv =  0.3594;
else
    % Southern station
    du =  0.0889;
    dv = -0.3594;
end


function [ReimagedPupilPosition] = GetReimagedPupilPositionWithBC(LabInputChannel)

switch LabInputChannel
    case 1
        ReimagedPupilPosition.u =  52.320;
        ReimagedPupilPosition.v = -36.903;
    case 2
        ReimagedPupilPosition.u =  52.560;
        ReimagedPupilPosition.v = -36.903;
    case 3
        ReimagedPupilPosition.u =  52.800;
        ReimagedPupilPosition.v = -32.170;
    case 4
        ReimagedPupilPosition.u =  53.040;
        ReimagedPupilPosition.v = -32.170;
    case 5
        ReimagedPupilPosition.u =  53.280;
        ReimagedPupilPosition.v = -27.438;
    case 6
        ReimagedPupilPosition.u =  53.520;
        ReimagedPupilPosition.v = -27.438;
    case 7
        ReimagedPupilPosition.u =  53.760;
        ReimagedPupilPosition.v = -22.705;
    case 8
        ReimagedPupilPosition.u =  54.000;
        ReimagedPupilPosition.v = -22.705;
end


function [ReimagedPupilPosition] = GetReimagedPupilPositionNoBC(LabInputChannel)

ReimagedPupilPosition.u =  46.560;

switch LabInputChannel
    case 1
        ReimagedPupilPosition.v = -34.215;
    case 2
        ReimagedPupilPosition.v = -35.175;
    case 3
        ReimagedPupilPosition.v = -33.975;
    case 4
        ReimagedPupilPosition.v = -34.935;
    case 5
        ReimagedPupilPosition.v = -33.735;
    case 6
        ReimagedPupilPosition.v = -34.695;
    case 7
        ReimagedPupilPosition.v = -33.495;
    case 8
        ReimagedPupilPosition.v = -34.455;
end



function [RadiusVCM] = ComputeRadiusOfCuravtureVCM(DistanceInputPupil, DistanceExitPupil)

% Optical parameters of the DL
% ----------------------------
f1  = 0.9;
f2  = -266.66e-03;
d12 = 0.7;
d23 = 0.8;

% f1 = optical focal length of M13 (DL primary)
% f2 = optical focal length of M14 (DL secondary)
% d12 = distance between vertices of M13 and M14 (measured along u-axis)
% d23 = distance between vertices of M14 and M15 (the VCM) (measured along u-axis)

% Do the computation
% ------------------

% Pin = distance between relayed telescope pupil and M13 vertex (measured along u-axis)
% Pout = distance between exit pupil and M13 vertex (measured along u-and v-axis (folded around M16))
% Din, Dout = image distances (measured w. r. t. M14 vertex) of the pupil image
% DinDash, DoutDash = same as Din and Dout, but measured w. r. t. the vertex of M15 (VCM)

Pin  = DistanceInputPupil + d12;
Pout = DistanceExitPupil + d12;

Din = f2 * (d12 - f1*Pin / (Pin - f1)) / ...
      ((d12 - f1*Pin / (Pin - f1)) - f2);

Dout = f2 * (d12 - f1*Pout / (Pout - f1)) / ...
       ((d12 - f1*Pout / (Pout - f1) - f2)); 

DinDash  = d23 - Din;

DoutDash = d23 - Dout;

RadiusVCM = -2*(DinDash * DoutDash) / ...
            (DinDash + DoutDash);



function [DistanceExitPupilFixedVCM] = ComputeOutputPupilDistance(DistanceInputPupil, RadiusFixedVCM)

% Optical parameters of the DL
% ----------------------------
f1  = 0.9;
f2  = -266.66e-03;
d12 = 0.7;
d23 = 0.8;

% f1 = optical focal length of M13 (DL primary)
% f2 = optical focal length of M14 (DL secondary)
% d12 = distance between vertices of M13 and M14 (measured along u-axis)
% d23 = distance between vertices of M14 and M15 (the VCM) (measured along u-axis)

% Do the computation
% ------------------

% Pin = distance between relayed telescope pupil and M13 vertex (measured along u-axis)
% Pout = distance between exit pupil and M13 vertex (measured along u-and v-axis (folded around M16))
% Din, Dout = image distances (measured w. r. t. M14 vertex) of the pupil image
% DinDash, DoutDash = same as Din and Dout, but measured w. r. t. the vertex of M15 (VCM)

Pin  = DistanceInputPupil + d12;

Din = f2 * (d12 - f1*Pin / (Pin - f1)) / ...
      ((d12 - f1*Pin / (Pin - f1)) - f2);

DinDash  = d23 - Din;
 
DoutDash = -1 * RadiusFixedVCM * DinDash / ...
            (RadiusFixedVCM + 2 * DinDash);

Dout = d23 - DoutDash;
 
Pout = f1 * (Dout * (d12 - f2) - d12 * f2) / ...
       (Dout * (d12 - f1 - f2) + f1 * f2 - f2 * d12);

DistanceExitPupilFixedVCM = Pout - d12;



function [RadiusVCM, PinMeter, PoutMeter] = ComputeRadiusOfCuravtureVCMGIP(DistanceInputPupil, DistanceExitPupil, DL)

% Optical parameters
% ------------------

switch DL
    case 1   	 %    (701.54, 801.01, 901.38, 266.20,0),
        m1m2 = -1.0 * 701.54;
        r1 = -2.0 * 901.38;
        r2 = -2.0 * 266.20;
        ofs = 0.0;  
    case 2       %    (701.03, 800.64, 901.24, 266.94,0),
        m1m2 = -1.0 * 701.03;
        r1 = -2.0 * 901.24;
        r2 = -2.0 * 266.94;
        ofs = 0.0;  
    case 3       %    (701.27, 800.93, 901.29, 266.48,0),
        m1m2 = -1.0 * 701.27;
        r1 = -2.0 * 901.29;
        r2 = -2.0 * 266.48;
        ofs = 0.0;  
    case 4       %    (700.51, 800.85, 900.28, 266.24,0)
        m1m2 = -1.0 * 700.51;
        r1 = -2.0 * 900.28;
        r2 = -2.0 * 266.24;
        ofs = 0.0;  
    case 5       %    (701.06, 801.42, 900.98, 266.41,0)
        m1m2 = -1.0 * 701.06;
        r1 = -2.0 * 900.98;
        r2 = -2.0 * 266.41;
        ofs = 0.0;  
    case 6       %    (700.45, 800.80, 900.30, 266.34,0)
        m1m2 = -1.0 * 700.45;
        r1 = -2.0 * 900.30;
        r2 = -2.0 * 266.34;
        ofs = 0.0;  
    case 7
        m1m2 = -1.0 * 700;
        r1 = -2.0 * 900;
        r2 = -2.0 * 266.66;
        ofs = 0.0;  
    case 8
        m1m2 = -1.0 * 700;
        r1 = -2.0 * 900;
        r2 = -2.0 * 266.66;
        ofs = 0.0;  
end


% Do the computation
% ------------------

% Pin = distance between relayed telescope pupil and M13 vertex (measured along u-axis)
% Pout = distance between exit pupil and M13 vertex (measured along u-and v-axis (folded around M16))
% Din, Dout = image distances (measured w. r. t. M14 vertex) of the pupil image
% DinDash, DoutDash = same as Din and Dout, but measured w. r. t. the vertex of M15 (VCM)

Pin  = DistanceInputPupil*1000 - m1m2 + ofs;
Pout = DistanceExitPupil*1000  - m1m2 + ofs;

F   = 0.5 * (r1 * r2) / (r1 - r2 - 2.0 * m1m2);
D   = r1 * (0.5 - (F / r2));

cur = (D + 0.5 * (Pin + Pout)) / (F * F);
roc = 1.0 / cur; % in mm

RadiusVCM = roc/1000; % in m

PinMeter = Pin/1000;
PoutMeter = Pout/1000;