Matlab example of CFL Calibration for PLab's SWB

swbcfl.m

% File:     swbcfl.m
% Org:      20150218 Dave Stoft
% Revid:    20150323
% Descr:    Example code to calculate pixel:wavelength mapping when pixel#
%           for 436nm and 546nm peaks are known
%
% Concept:
% With an incidence angle to the grating of 45-deg, and the camera image
% plane thus parallel to the grating, the 524nm diffraction angle will be
% 0-deg; relative to the camera normal (perpendicular to the image plane
% and through the center of the lens). The 524nm pixel will be roughly at
% the mid-point 'X' of the camera image.
%
% As the wavelength extends below and above 524nm, the change in nm per
% pixel change will, itself, change -- it is not constant. However, the
% rate of change is linear.
%
% Note: When limited to the PLab spectrometer's useful range of 400-650nm,
% the correction per pixel is relatively small - on the order of 2nm. Is
% this still worthwhile to do? I think so, just because 1) it can be
% corrected, 2) that's what computers are good at and 3) it would seem to
% be totally transparent to the user.
%
% The spectrum will be spread either side of 524nm and the 436nm
% and 546nm CFL peaks will appear with 'X' pixel values below and above
% this 'center'. However, since the CFL has no peak at 524nm, the 436nm
% and 546nm peaks are used to calculate the 524nm pixel via simple trig
% functions. If the incedence angle is not exactly 45-deg, an error can
% result, but that error is much smaller than the amount of this
% correction so can probably be ignored.
%
% Given this configuration, the relationship between wavelength,
% diffraction angle and pixel 'X' position in the image plane can be
% calculated. Given only the pixel# of the CFL 436nm and 546nm peaks, a 1:1
% array can thus be generated which maps pixels to wavelength(nm).
%
% Assume:
% - Incidence angle (IA) from slit to diffraction grating is 45-deg
% - Camera normal is perpendicular to the diffraction grating plane
% - Visible spectrum is roughly centered within the camera image plane
% - The spectrum 'average' data has been taken and stored in an array
% - Peaks at 436nm (blue) and 546nm (upper green of doublet) are known
% - The pixel numbers of the max value of 436 and 546 peaks are known
%
% Calc:
% - With an IA of 45-deg, 524nm has a rel diffraction angle (DA) of 0-deg
% - For any pixel, we can calc the pixel 'x' offset to get the angle
% - Then, back-calc from the angle to get the wavelength
% - Using known 436 and 546 pixels, use trig to find pixel# of 524nm
% - Then, calc the wavlength for any pixel
% - So we get a 1:1 mapping of pixels to nm !

% ------------------------------------------------------------------------
function swbcfl( )

% ------------------------------------------------------------------------
% Declarations:
IA = 45.00;             % From spectrometer construction
DGnm = 740;              % Diffraction grating line spacing in nm
Zero_nm = 523.78;       % Ideal center for 45.00 deg IA
Cal1_nm = 435.83;       % Blue Hg line
Cal2_nm = 546.08;       % Higher-nm peak of Grn Hg doublet
Pix436 = Pixel436;      % Image pixel# of 436nm peak
Pix546 = Pixel546;      % Image pixel# of 546nm peak
PixMax = 1600;          % Total 'X' pix of the camera image

% ------------------------------------------------------------------------
% DEBUG: Read the REF image data from the file
[cfl_image, w_pix, h_pix] = read_image( [ 'P3U-CFL-03_TLapse.jpg'] );

% DEBUG: Read/average three line of colors from the image
y_line = 630;
for x_pix = 1:1600
    rsum = 0;
    gsum = 0;
    bsum = 0;
    for y_pix = -1:+1
        rsum = rsum + int32(cfl_image(y_line+y_pix,x_pix,1));
        gsum = gsum + int32(cfl_image(y_line+y_pix,x_pix,2));
        bsum = bsum + int32(cfl_image(y_line+y_pix,x_pix,3));
    end
    pix(x_pix,1) = int16(double(rsum)/double(3));
    pix(x_pix,2) = int16(double(gsum)/double(3));
    pix(x_pix,3) = int16(double(bsum)/double(3));
    cfl_uimage(x_pix) = ( pix(x_pix,1)+pix(x_pix,2)+pix(x_pix,3) ) / 3.0;
end

% ------------------------------------------------------------------------
% First, knowing the 436 and 546 pixels, find the 524nm pixel
%
%          DVD
%  |       /|\
%  |      / | \
%  |     /  |  \
%  FD   / A | B \
%  |   /    |    \
%  |  /     |     \
%    --------------- Image
%    |DA436 | DA546|
%  X436   X524   X546
%
% Notice that angles A and B share a common adjacent side
%   Tan(A) = (X524-X436) / FD    % A=abs(DA436), FD=focal distance
%   Tan(B) = (X546-X524) / FD    % B=abs(DA546)
% Since FD is in common
%   (X524-X436)/Tan(A) = (X546-X524)/Tan(B)
% Then....
% Find the number of pixels between 436 and 546nm
%   P = (Pix546-Pix436)
%   So, X546 = P - X436
% And, by simple substitution...
%   X436/Tan(A) = (P-X436)/Tan(B)
% Then, derive the pixel for 524nm...
%   X436 = Tan(A)/Tan(B) * (P-X436)
%   X436 = P*(Tan(A)/Tan(B)) - P436*(Tan(A)/Tan(B))
%   X436 * (1+(Tan(A)/Tan(B)) = P*(Tan(A)/Tan(B))
%   X436 = [ P*(Tan(A)/Tan(B)) ] / [ 1+(Tan(A)/Tan(B)) ]
%   X546 = P - X436
%   X524 = X436 + X436

% Next, derive the formula for diffraction angle
%   Include the effect of 45-deg incidence angle
%
% Start with the basic diffraction equasion for non-zero IA
%   d(sin(a) - sin(i)) = m*L
% where....
%   d = grating line spacing
%   a = grating exit angle
%   i = incidence angle
%   m = diffraction mode (1 here)
%   L = lambda = wavelength of light
% Solving for angle (b)
%   L / d = sin(a) - sin(i)
%   sin(a) = (L/d) + sin(i)
%   a = asin( (L/d) + sin(i) )
% Because of the polarity of the angle ...
%   DA(nm) = asin( (Lnm/DGnm) - sin( IA*(2*pi/360) ) ) * (360/(2*pi));

% ========================================================================
% ========================================================================
% Now, do the calculations for real

% Diffraction angle for 524nm -- we've declared it to be zero
DA524 = 0.0;
% Diffraction angle for 436nm (trig functions work in radians)
DA436 = abs( asin( (Cal1_nm/DGnm) - sin(IA*(2*pi/360) ) ) * (360/(2*pi)) );
% Diffraction angle for 546nm
DA546 = abs( asin( (Cal2_nm/DGnm) - sin(IA*(2*pi/360) ) ) * (360/(2*pi)) );

% Find the image pixel span from Cal points: 436nm to 546nm
PixSpan = Pix546 - Pix436;
% Find the tangent ratio -- because the ZeroDA distance is common
TanAB = (tan(DA436*(2*pi/360))) / (tan(DA546*(2*pi/360)));
% Calc the #pix on the 436nm side -- the pix distance at image from 524nm
NumPix436 = int16( (PixSpan*TanAB) / (1+TanAB) );
% Now calc the ZeroDA distance which is in common
ZPixFD = double(NumPix436) / ( tan(DA436*(2*pi/360)) );
% Find the #pix from ZeroDA to 546nm image pix
NumPix546 = PixSpan - Pix436;
% Find the image pix for theoretical ZeroDA 524nm
Pix524 = Pix436 + NumPix436;

% Now, we have all we need to calc all pixle wavelengths
%
% Calc all Negative diffraction angles (per pixel)
f400 = 0;       % debug use
for pix = Pix524:-1:1
    % We now know focal distance 'ZPixFD' so we can find the diffraction
    % angle for any pixel distance relative to 524nm
    % Remember, for nm<ZeroDA, the diffraction angles are all negative
    %
    % First, find the DA angle relative to the ZeroDA (camera normal)
    %   L = atan( pixel_offset_from_524nm / FD )
    DA(pix) = atan( double(Pix524-pix) / ZPixFD ) * (360/(2*pi));
    % Then find the nm for this angle
    % Start with the basic diffraction equasion for non-zero IA
    %   d(sin(a) - sin(i)) = m*L
    %   L = d(sin(a) - sine(i))
    %   NM(pix) = DGnm * ( sin(DA(pix)*(2*pi/360)) - sin( IA*(2*pi/360) )
    PIXnm(pix) = DGnm * ...
                    ( sin(IA*(2*pi/360)) - sin(DA(pix)*(2*pi/360)) );
    % DEBUG: Identify the 400nm pixel
    if ( f400 == 0 )
        if ( PIXnm(pix) < 400 )
            Pix400 = pix;
            f400 = 1;
        end
    end
end

% Calc all Positive diffraction angles (per pixel)
f650 = 0;
for pix = (Pix524+1):PixMax
    DA(pix) = atan( double(pix-Pix524) / ZPixFD ) * (360/(2*pi));
    PIXnm(pix) = DGnm * ...
                    ( sin(IA*(2*pi/360)) + sin(DA(pix)*(2*pi/360)) );
    % DEBUG: Identify the 650nm pixel
    if ( f650 == 0 )
        if ( PIXnm(pix) > 650 )
            Pix650 = pix;
            f650 = 1;
        end
    end
end

% PIXnm(1:1600) now contains the (nm) values for every 'X' pixel in the
% spectrum. This is the end of the necessary calculations.
% ========================================================================
% ========================================================================

% ------------------------------------------------------------------------
% DEBUG: Find the near-linear nm/pix rate of change (look near 524nm)
Zslope = (PIXnm(Pix524+10)-PIXnm(Pix524-10)) / 20.0;

% DEBUG: Calc the offset from a linear estimate due to disperaion
for pix = Pix524:-1:1
    % Linear estimate (nm) = 524nm - (nm/pix * pix_offset_from_524nm)
    linear_est(pix) = double(PIXnm(Pix524)) - (Zslope * double(Pix524-pix));
    % Find nm difference from linear est
    nmdif(pix) = double(PIXnm(pix)) - linear_est(pix);
end
for pix = Pix524:1600
    linear_est(pix) = double(PIXnm(Pix524)) + (Zslope * double(pix-Pix524));
    nmdif(pix) = linear_est(pix) - double(PIXnm(pix));
end

pmin = Pix400;
pmax = Pix650;

% ------------------------------------------------------------------------
% DEBUG: Plot CFL and Dispersion on same plot

figure( 'position', [ 200 400 800 350] );
[AX,H1,H2] = ...
    plotyy( (pmin:pmax), cfl_uimage(pmin:pmax), (pmin:pmax), nmdif(pmin:pmax) );

xlabel( 'Pixel Position', 'FontSize',12 );

set(get(AX(1),'Ylabel'),'String','Intensity (unitless)');
set(get(AX(2),'Ylabel'),'String','Dispersion(nm)');

set(get(AX(1),'Ylabel'),'FontSize',14);
set(get(AX(2),'Ylabel'),'FontSize',14);

%ylabel( 'Intensity (unitless)', 'FontSize',12 );

title( 'PLab3U CFL Spectra + Dispersion', 'FontSize', 12 );

set(AX(1),'xlim',[pmin pmax]);
set(AX(1),'ylim',[0 150]);

set(AX(2),'xlim',[pmin pmax]);
set(AX(2),'ylim',[-1 +2]);

grid on;

return;

% ------------------------------------------------------------------------
% DEBUG: Read the image data
function [new_image, w_pix, h_pix] = read_image( image_name )

% Fetch the image from the file
new_image = imread( image_name );

% Find the image dimensions
[ h_pix, w_pix, colors ] = size(new_image);

return;