%fix analyzer, vary qwp angle

analyzer_angle    = 45;

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

maxV              = 0;
Vvec              = [];
powavec           = [];
powbvec           = [];

anglevec          = linspace(0,45,10);

for qwp_angle = anglevec
   outstatea      = pol_rotate(hlp,analyzer_angle)*pol_rotate(qwp,qwp_angle)*instatea
   state          = outstatea;
   figure(1);pol_locus_plot
   title('A')
   pause
   totpowa        = sum(abs(outstatea).^2);
   powavec        = [powavec totpowa];
   outstateb      = pol_rotate(hlp,analyzer_angle)*pol_rotate(qwp,qwp_angle)*instateb
   state          = outstateb;
   figure(2);pol_locus_plot
   title('B')
   pause
   totpowb        = sum(abs(outstateb).^2);
   powbvec        = [powbvec totpowb];
   visibility     = abs( (totpowa-totpowb)/(totpowa+totpowb) );
   Vvec           = [Vvec visibility];
   if maxV < visibility
      maxV        = visibility;
      maxVangle   = qwp_angle;
   end
   figure(3);plot(anglevec(1:length(powavec)),[powavec;powbvec;Vvec]');
end

maxpow= max([powavec powbvec])

figure;plot(anglevec,[powavec/maxpow;powbvec/maxpow;Vvec]');axis([0 180 0 1])

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