opt_absorption.m

Matlab实现光线跟踪算法

function damping = opt_absorption(glass,wavelength)
% OPT_ABSORPTION - Calculate the absorption coefficient
% The data is my reading of Pr鋤isions Glas & Optik GmbH
% Table of properties of common glasses. Available on:
% http://www.pgo-online.com/intlframes/produkteset.html
% 
% Calling:
% damping = opt_absorption(glass,wavelength)

persistent glass_names absorption wl_named

if isempty(glass_names)
  
  qwe = '';
  fp = fopen('Glass.trans','r');
  while ~feof(fp)
    qwe = str2mat(qwe,fgetl(fp));
  end
  glass_names = qwe(3:end,1:7);
  transmission = str2num(qwe(3:end,8:end));
  wl_named = str2num(qwe(2,8:end));
  
end

I = strmatch(glass,glass_names);
if ~isempty(I)
  lambda = wl_named;
  alpha = log(1-transmission(I,:))/10;% units!!! ~ /mm?
else
  switch lower(glass)
   case 'b270'
    transmission = [.12424 .83731 .88507 .9 .9 .9 .9 .9 .9 .9 .9 .9 .9 .9 .9];
    alpha = [996.58 43.172 15.609 7.5415 7.646 8.3044 8.6947 8.941 ...
             9.1052 9.2198 9.3027 9.3646 9.4119 9.4489 9.4784];
    lambda =  [300 343.53 358.82 388.24 400:100:1400]*1e-9;
   case 'bk7'
    transmission = [0.09802 0.2981 0.50024 0.69869 0.85628 0.89346 ...
                    0.91346 0.91918 0.92775 0.93061 0.92204 0.92489 ...
                    0.93061 0.89918 0.87345];
    alpha = [223.02 111.83 60.079 26.696 6.409 2.1927 0.044224 0 0 0 ...
             0 0 0 2.0032 4.9095];
    lambda = [300 310.53 317.76 328.29 350 365.33 400 450 1000 1375 ...
              1453.1 1500 1546.2 1853.8 2000]*1e-9;
   case 'fusedsilica'
    transmission = [0.051752 0.19597 0.40228 0.60006 0.80112 0.81296 ...
                    0.82184 0.86032 0.88993 0.91362 0.93139 0.94028 ...
                    0.94028 0.94324 0.86328 0.81296 0.93436 0.93139 ...
                    0.92547 0.90177 0.61779 0.60597 0.61484 0.80112 ...
                    0.83368 0.84256 0.83664 0.60006 0.037045 0.022341 ...
                    0.042928 0.60006 0.68281 0.74491 0.77154 0.65917 ...
                    0.51146 0.43178 0.38163 0.20775 0.1518 0.054693 ...
                    0.0076379];
    alpha = [287.06 153.86 81.879 41.822 12.856 11.377 10.268 5.5872 ...
             2.2029 0 0 0 0 0 6.054 12.063 0 0 0 1.7196 39.541 41.473 ...
             40.022 13.557 9.5739  8.515 9.2205 42.457 320.95 371.52 ...
             306.21 42.461 29.542 20.838 17.327 33.069 58.439 75.376 ...
             87.722 148.54 179.91 281.99 478.86]; 
    lambda = [162.31 164.64 167.53 171.55 176.68 177.81 180.07 200.47 ...
              220.02 243.52 386.01 550 710 1210.9 1281.2 1334 1474.9 ...
              1828.7 2006.4 2113.2 2131 2166.7 2220.2 2291.6 2309.5 ...
              2363.1 2398.9 2470.6 2650 2739.9 2848 2974.3 3046.5 3137 ...
              3282 3427.3 3500.1 3591.2 3627.6 3664.1 3920.1 4140.3 ...
              4342.8]*1e-9;
   case 'pyrex'
    transmission = [0.099787 0.79775 0.90059 0.90994 0.91306 0.90371 ...
                    0.88188 0.88188 0.89124 0.89747 0.90059 0.90059 ...
                    0.88188 0.85383 0.8289 0.80087 0.75728 0.69815 ...
                    0.60176 0.52721 0.39691 0.33183 0.30085 0.2637 ...
                    0.25132 0.2606 0.27918 0.30085 0.35042 0.39691 ...
                    0.42481 0.43102 0.42481 0.40001 0.35351 0.30085 ...
                    0.19871 0.099787 0.050366 0.03184];
    alpha = [221.09 13.282 1.2425 0.32018 0.066598 1.3191 3.8822 3.9096 ...
             2.877 2.1863 1.8445 1.8482 3.9513 7.186 10.15 13.591 ...
             19.189 27.319 42.176 55.402 83.792 101.7 111.5 124.68 ...
             129.49 125.86 118.98 111.5 96.253 83.795 77 75.551 77 ...
             83.017 95.373 111.5 152.98 221.86 290.24 336.1];
    lambda = [253.2 274.46 302.81 345.32 387.83 586.06 925.43 1130.2 ...
              1454.6 1616.5 1785.4 1947.1 2200 2368.4 2466.6 2536.6 ...
              2592.7 2634.7 2669.7 2697.8 2725.8 2753.8 2781.8 2823.7 ...
              2914.7 2998.6 3033.6 3068.5 3117.4 3145.4 3166.3 3194.3 ...
              3229.2 3264.1 3306 3340.9 3410.6 3522.2 3605.9 3682.5]*1e-9;
   case 'air'
    lambda = [0 10000]*1e-9;
    alpha = [0 0];
   otherwise
    lambda = [0 10000]*1e-9;
    alpha = [0 0];
    warning(['No values for glass absorption for: ',glass])
  end

end

damping = interp1(lambda,alpha,wavelength);
ii = find(~isfinite(damping));
damping(ii) = inf;