1## Copyright (C) 2013 Martin Vogel <octave@martin-vogel.info> 2## 3## This program is free software; you can redistribute it and/or modify it 4## under the terms of the GNU General Public License as published by the 5## Free Software Foundation; either version 3 of the License, or (at your 6## option) any later version. 7## 8## This program is distributed in the hope that it will be useful, but WITHOUT 9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 11## for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with this program; see the file COPYING. If not, see 15## <http://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {@var{M} =} mueller_lindiattenuator() 19## @deftypefnx {Function File} {@var{M} =} mueller_lindiattenuator(@var{d}) 20## @deftypefnx {Function File} {@var{M} =} mueller_lindiattenuator(@var{px},@var{py}) 21## @deftypefnx {Function File} {@var{M} =} mueller_lindiattenuator(..., @var{mode}) 22## Return the Mueller matrix for a linear diattenuator at zero 23## rotation. 24## 25## @itemize @minus 26## @item @var{d} is the diattenuation of the element, i.e. 27## @code{d=(px-py)/(px+py)}. Reversibly, transmission in y direction 28## is @code{(1-d)/(1+d)}, if transmission in x direction is 1. 29## @item @var{px} is the transmittance in x direction. 30## @item @var{py} is the transmittance in y direction. 31## @item @var{mode} is a string defining the interpretation of 32## transmittance values: 'intensity' (default) or 'amplitude'. 33## @end itemize 34## 35## Arguments @var{d}, @var{px} or @var{py} can be passed as a scalar 36## or as a matrix or as a cell array. In the two latter cases, a cell 37## array @var{M} of Mueller matrices is returned. The size of @var{M} 38## is set to the maximum of the parameters' size. 39## 40## References: 41## 42## @enumerate 43## @item E. Collett, Field Guide to Polarization, 44## SPIE Field Guides vol. FG05, SPIE (2005). ISBN 0-8194-5868-6. 45## @item R. A. Chipman, "Polarimetry," chapter 22 in Handbook of Optics II, 46## 2nd Ed, M. Bass, editor in chief (McGraw-Hill, New York, 1995) 47## @item @url{http://en.wikipedia.org/wiki/Mueller_calculus, "Mueller calculus"}, 48## last retrieved on Dec 17, 2013. 49## @end enumerate 50## 51## @seealso{mueller_circdiattenuator, mueller_absorber} 52## @end deftypefn 53 54function M = mueller_lindiattenuator(varargin) 55 56 if nargin<1 57 58 px = 1; 59 py = 1; 60 was_cell = false; 61 62 elseif nargin==1 63 64 px = varargin{1}; 65 [px, was_cell] = __c2n__(px, 0); 66 py = (1-px)./(1+px); 67 px(:) = 1; 68 69 else 70 71 px = varargin{1}; 72 [px, was_cellx] = __c2n__(px, 1); 73 74 py = varargin{2}; 75 [py, was_celly] = __c2n__(py, 1); 76 77 was_cell = was_cellx || was_celly; 78 79 end 80 81 % check mode 82 s_function = @s_lindiattenuator_int; 83 if nargin>=2 && ischar(varargin{end}) 84 if strncmpi(varargin{end},'amp',3) 85 s_function = @s_lindiattenuator_amp; 86 end 87 end 88 89 % any matrix in parameters? 90 if (any([numel(px),numel(py)] > 1)) || was_cell 91 92 % adjust dimensions, i.e. fill missing dimensions with 1 93 spx = size(px); 94 spy = size(py); 95 96 maxdim = max([length(spx),length(spy)]); 97 if length(spx) < maxdim 98 spx = [spx, ones(1,maxdim-length(spx))]; 99 end 100 if length(spy) < maxdim 101 spy = [spy, ones(1,maxdim-length(spy))]; 102 end 103 104 % generate Mueller matrices 105 maxsize = max([spx;spy]); 106 M = cell(maxsize); 107 M_subs = cell(1,ndims(M)); 108 numelM = numel(M); 109 110 % flatten parameter arrays 111 px = px(:); 112 py = py(:); 113 numelpx = numel(px); 114 numelpy = numel(py); 115 116 for mi=1:numelM 117 [M_subs{:}] = ind2sub(size(M),mi); 118 M{M_subs{:}} = s_function(px(mod(mi-1,numelpx)+1), py(mod(mi-1,numelpy)+1)); 119 end 120 121 else 122 123 M = s_function(px, py); 124 125 end 126 127end 128 129% helper function 130function M = s_lindiattenuator_amp(px_in_amplitude_domain,py_in_amplitude_domain) 131 132 M = zeros(4,4); 133 134 M(1,1) = (px_in_amplitude_domain^2+py_in_amplitude_domain^2)/2; 135 M(1,2) = (px_in_amplitude_domain^2-py_in_amplitude_domain^2)/2; 136 M(2,1) = M(1,2); 137 M(2,2) = M(1,1); 138 M(3,3) = px_in_amplitude_domain*py_in_amplitude_domain; 139 M(4,4) = M(3,3); 140 141end 142 143% helper function 2 144function M = s_lindiattenuator_int(kx_in_intensity_domain,ky_in_intensity_domain) 145 146 M = zeros(4,4); 147 148 M(1,1) = (kx_in_intensity_domain+ky_in_intensity_domain)/2; 149 M(1,2) = (kx_in_intensity_domain-ky_in_intensity_domain)/2; 150 M(2,1) = M(1,2); 151 M(2,2) = M(1,1); 152 M(3,3) = sqrt(kx_in_intensity_domain*ky_in_intensity_domain); 153 M(4,4) = M(3,3); 154 155end 156 157%!test 158%! % test default return value 159%! A = mueller_lindiattenuator(); 160%! R = A*A-A; 161%! assert(norm(R,inf), 0, 1e-9); 162%! 163%!test 164%! % test equality of providing diattenuation or kx and ky 165%! d = rand(1, 1); 166%! kx = 1; 167%! ky = (1-d)./(1+d); 168%! A1 = mueller_lindiattenuator(d); 169%! A2 = mueller_lindiattenuator(kx, ky); 170%! R = A1-A2; 171%! assert(norm(R,inf), 0, 1e-9); 172%! 173%!test 174%! % test serial application of linear diattenuators 175%! kx = rand(1, 1); 176%! ky = rand(1, 1); 177%! A1 = mueller_lindiattenuator(kx,ky); 178%! A2 = mueller_lindiattenuator(kx*kx,ky*ky); 179%! R = A1*A1-A2; 180%! assert(norm(R,inf), 0, 1e-9); 181%! 182%!test 183%! % another test of serial application of retarder elements 184%! kx1 = rand(1, 1); 185%! kx2 = rand(1, 1); 186%! ky1 = rand(1, 1); 187%! ky2 = rand(1, 1); 188%! A1 = mueller_lindiattenuator(kx1,ky1); 189%! A2 = mueller_lindiattenuator(kx2,ky2); 190%! A12 = mueller_lindiattenuator(kx1*kx2,ky1*ky2); 191%! R = A1*A2-A12; 192%! assert(norm(R,inf), 0, 1e-9); 193%! 194%!test 195%! % test mode 196%! kx = rand(1, 1); 197%! ky = rand(1, 1); 198%! A1 = mueller_lindiattenuator(kx,ky); 199%! A2 = mueller_lindiattenuator(kx,ky,'int'); 200%! A3 = mueller_lindiattenuator(sqrt(kx),sqrt(ky),'amp'); 201%! R1 = A1-A2; 202%! R2 = A1-A3; 203%! assert(norm(R1,inf)+norm(R2,inf), 0, 1e-9); 204%! 205%!test 206%! % test correct size of return values 207%! for dim = 1:5 208%! asize = randi([1 4], 1, dim); 209%! R = rand(asize); 210%! if numel(R) == 1 211%! R = {R}; 212%! end 213%! rsize = size(R); 214%! C = mueller_lindiattenuator(R); 215%! csize = size(C); 216%! assert(rsize == csize); 217%! end 218%! 219%!test 220%! % another test correct size of return values 221%! kx = rand(3,4,5); 222%! ky = rand(5,4,3); 223%! C = mueller_lindiattenuator(kx,ky); 224%! csize = size(C); 225%! assert(csize == [5,4,5]); 226 227 228 229