%fix analyzer, vary qwp angle

analyzer_angle    = 45;

instatea          = [1; 1];

maxV              = 0;
Vvec              = [];

anglevec          = linspace(0,90,20);
delphasevec       = 5*linspace(-180,180,100);

for qwp_angle = anglevec
   powvec         = [];        
   for cosphase   = delphasevec
      instate     = [1 0; 0 exp(j*cosphase*pi/180)]*instatea;
      outstate    = pol_rotate(hlp,analyzer_angle)*pol_rotate(qwp,qwp_angle)*instate;
      totpow      = sum(abs(outstate).^2);
      powvec      = [powvec totpow];
   end
   dispmat        = ones(50,1)*powvec;
   figure(1);pcolor(dispmat);shading interp;colormap(gray);caxis([0 max(powvec)])
   pause(0.5)
   visibility     = abs( (max(powvec) - min(powvec) )/(max(powvec)+min(powvec)) );
   Vvec           = [Vvec visibility];
   if maxV < visibility
      maxV        = visibility;
      maxVangle   = qwp_angle;
   end
   figure(2);plot(anglevec(1:length(Vvec)),Vvec);
   title('Visibility')
end


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