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_circdiattenuator()
19## @deftypefnx {Function File} {@var{M} =} mueller_circdiattenuator(@var{d})
20## @deftypefnx {Function File} {@var{M} =} mueller_circdiattenuator(@var{pl},@var{pr})
21## @deftypefnx {Function File} {@var{M} =} mueller_circdiattenuator(..., @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{pl} is the transmittance in x direction.
30## @item @var{pr} 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{pl} or @var{pr} 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_lindiattenuator, mueller_absorber}
52## @end deftypefn
53
54function M = mueller_circdiattenuator(varargin)
55
56  if nargin<1
57
58    pr = 1;
59    pl = 1;
60    was_cell = false;
61
62  elseif nargin==1
63
64    pr = varargin{1};
65    [pr, was_cell] = __c2n__(pr, 0);
66    pl = (1-pr)./(1+pr);
67    pr(:) = 1;
68
69  else
70
71    pr = varargin{1};
72    [pr, was_cellr] = __c2n__(pr, 1);
73
74    pl = varargin{2};
75    [pl, was_celll] = __c2n__(pl, 1);
76
77    was_cell = was_cellr || was_celll;
78
79  end
80
81  % check mode
82  s_function = @s_circdiattenuator_int;
83  if nargin>=2 && ischar(varargin{end})
84    if strncmpi(varargin{end},'amp',3)
85      s_function = @s_circdiattenuator_amp;
86    end
87  end
88
89  % any matrix in parameters?
90  if (any([numel(pr),numel(pl)] > 1)) || was_cell
91
92    % adjust dimensions, i.e. fill missing dimensions with 1
93    spr = size(pr);
94    spl = size(pl);
95
96    maxdim = max([length(spr),length(spl)]);
97    if length(spr) < maxdim
98      spr = [spr, ones(1,maxdim-length(spr))];
99    end
100    if length(spl) < maxdim
101      spl = [spl, ones(1,maxdim-length(spl))];
102    end
103
104    % generate Mueller matrices
105    maxsize = max([spr;spl]);
106    M = cell(maxsize);
107    M_subs = cell(1,ndims(M));
108    numelM = numel(M);
109
110    % flatten parameter arrays
111    pr = pr(:);
112    pl = pl(:);
113    numelpr = numel(pr);
114    numelpl = numel(pl);
115
116    for mi=1:numelM
117      [M_subs{:}] = ind2sub(size(M),mi);
118      M{M_subs{:}} = s_function(pr(mod(mi-1,numelpr)+1), pl(mod(mi-1,numelpl)+1));
119    end
120
121  else
122
123    M = s_function(pr, pl);
124
125  end
126
127end
128
129% helper function
130function M = s_circdiattenuator_amp(pr_in_amplitude_domain,pl_in_amplitude_domain)
131
132  M = zeros(4,4);
133
134  M(1,1) = (pr_in_amplitude_domain^2+pl_in_amplitude_domain^2)/2;
135  M(1,4) = (pr_in_amplitude_domain^2-pl_in_amplitude_domain^2)/2;
136  M(2,2) = pr_in_amplitude_domain*pl_in_amplitude_domain;
137  M(3,3) = M(2,2);
138  M(4,1) = M(1,4);
139  M(4,4) = M(1,1);
140
141end
142
143% helper function 2
144function M = s_circdiattenuator_int(kr_in_intensity_domain,kl_in_intensity_domain)
145
146  M = zeros(4,4);
147
148  M(1,1) = (kr_in_intensity_domain+kl_in_intensity_domain)/2;
149  M(1,4) = (kr_in_intensity_domain-kl_in_intensity_domain)/2;
150  M(2,2) = sqrt(kr_in_intensity_domain*kl_in_intensity_domain);
151  M(3,3) = M(2,2);
152  M(4,1) = M(1,4);
153  M(4,4) = M(1,1);
154
155end
156
157%!test
158%! % test default return value
159%! A = mueller_circdiattenuator();
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%! kr = 1;
167%! kl = (1-d)./(1+d);
168%! A1 = mueller_circdiattenuator(d);
169%! A2 = mueller_circdiattenuator(kr, kl);
170%! R = A1-A2;
171%! assert(norm(R,inf), 0, 1e-9);
172%!
173%!test
174%! % test serial application of circular diattenuators
175%! kr = rand(1, 1);
176%! kl = rand(1, 1);
177%! A1 = mueller_circdiattenuator(kr,kl);
178%! A2 = mueller_circdiattenuator(kr*kr,kl*kl);
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%! kr1 = rand(1, 1);
185%! kr2 = rand(1, 1);
186%! kl1 = rand(1, 1);
187%! kl2 = rand(1, 1);
188%! A1 = mueller_circdiattenuator(kr1,kl1);
189%! A2 = mueller_circdiattenuator(kr2,kl2);
190%! A12 = mueller_circdiattenuator(kr1*kr2,kl1*kl2);
191%! R = A1*A2-A12;
192%! assert(norm(R,inf), 0, 1e-9);
193%!
194%!test
195%! % test mode
196%! kr = rand(1, 1);
197%! kl = rand(1, 1);
198%! A1 = mueller_circdiattenuator(kr,kl);
199%! A2 = mueller_circdiattenuator(kr,kl,'int');
200%! A3 = mueller_circdiattenuator(sqrt(kr),sqrt(kl),'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_circdiattenuator(R);
215%!   csize = size(C);
216%!   assert(rsize == csize);
217%! end
218%!
219%!test
220%! % another test correct size of return values
221%! kr = rand(3,4,5);
222%! kl = rand(5,4,3);
223%! C = mueller_circdiattenuator(kr,kl);
224%! csize = size(C);
225%! assert(csize == [5,4,5]);
226