Created
May 18, 2015 17:23
-
-
Save Stoft1/2f3fd173f1ea3ccb2b45 to your computer and use it in GitHub Desktop.
Matlab example of CFL Calibration for PLab's SWB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
% 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; | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment