% vary analyzer angle; no quarter-wave plate

instatea          = [1; ];
instateb          = [1 0 ; 0 -1]*instatea; %180 degrees out of phase

maxV              = 0;
Vvec              = [];

anglevec          = linspace(0,90,100);

for analyzer_angle = anglevec
   outstatea      = pol_rotate(hlp,analyzer_angle)*instatea;
   totpowa        = sum(abs(outstatea).^2);
   outstateb      = pol_rotate(hlp,analyzer_angle)*instateb;
   totpowb        = sum(abs(outstateb).^2);
   visibility     = abs( (totpowa-totpowb)/(totpowa+totpowb) );
   Vvec           = [Vvec visibility];
   if maxV < visibility
      maxV        = visibility;
      maxVangle   = analyzer_angle;
   end
end

figure;plot(anglevec,Vvec);axis([0 90 0 1])

disp(['Maximum visibility = ' num2str(maxV) '  occurs at analyzer angle = ' num2str(maxVangle) ' degrees'])