% pulse with resonance error, w != w0
%
%  c0(t) =   c0(0)*(cos(d*t) - (i*a/d)*sin(d*t)) + m*c1(0)*(i*b/d)*exp(i*p)*sin(d*t)
%
%  c1(t) = m*c1(0)*(cos(d*t) + (i*a/d)*sin(d*t)) +   c0(0)*(i*b/d)*exp(-i*p)*sin(d*t)
%
% for derivation of formulas without resonance error, see Eq. (13.9) page 79 of
%  Introduction to Quantum Computers,
%  https://www.worldscientific.com/worldscibooks/10.1142/3808
%  Gennady P Berman, Gary D Doolen, Ronnie Mainieri, Vladimir I Tsifrinovich,
%  World Scientific, 1998.
%
% here p=phi (phase shift), a=(w-w0)/2, b=Omega/2, d=sqrt(a*a+b*b), rb=b/d, ra=a/d
%
% m=1, p=0 for pi-pulse, applied until time t1 = pi/Omega = pi/(2*b)
%
% m=-1, p=pi/2 for Hadamard, with pulse applied until time t2 = pi/(2*Omega) = pi/(4*b)
%
% ra=0, rb=1 correspond to no resonance error
%---
if( !exist( "ex") )
  disp( "set ex = 1 or 2 or 3 or 4")
  return
end

ex

if( ex == 1) % example 1, pi-pulse with resonance error
  %
  %  at t=0: c0(0)=1, c1(0)=0
  %  m=1, p=0, and pulse is applied until time t1 = pi/Omega = pi/(2*b)
  %  let A = pi/(2*rb), then:
  %   c0(t1) = cos(A) - i*ra*sin(A)
  %   c1(t1) = i*rb*sin(A)
  %
  ra = [-50:50]*(sqrt(3)/2/50); % min rb = 0.5
  rb = sqrt(1-ra.*ra);
  %
  A = pi./(2*rb);
  S = i*sin(A);
  C = cos(A);
  c0 = C-ra.*S;
  c1 = rb.*S;
  %
  hold off
  plot(ra,imag(c1),"-.b")
  % set (gca (), "xdir", "reverse") 
  % set (gca(),'interpreter','tex')
  hold on
  grid on
  plot(ra,imag(c0),"--r")
  plot(ra,real(c1),":g")
  plot(ra,real(c0),"-k")
  legend( "imag(c1)", "imag(c0)", "real(c1)", "real(c0)" ) 
  title '(a) pi-pulse result with resonance error' % must use single quotes to get \pi properly
  xlabel "a/d"
  ylabel 'real, imag'
  print -dpng pulse1.png

elseif( ex == 2) % example 2, Hadamard with resonance error

  %  at t=0: c0(0)=c1(0)=1/sqrt(2)
  %  m=-1, p=pi/2 and pulse is applied until time t2 = pi/(2*Omega) = pi/(4*b)
  %  let A = pi/(4*rb), then:
  %   c0(t2) = ( cos(A) - i*ra*sin(A) + rb*sin(A)) / sqrt(2)
  %   c1(t2) = (-cos(A) - i*ra*sin(A) + rb*sin(A)) / sqrt(2)
  %
  ra = [-50:50]*(sqrt(3)/2/50); % min rb = 0.5
  rb = sqrt(1-ra.*ra);
  %
  r2 = sqrt(2)/2;
  A = pi./(4*rb);
  C = r2*cos(A);
  S = r2*(rb-i*ra).*sin(A);
  c0 = S+C;
  c1 = S-C;
  % note: imag(c0) == imag(c1)
  hold off
  plot(ra,real(c0),"-k")
  hold on
  grid on
  plot(ra,imag(c0),"--r")
  plot(ra,imag(c1),"-.b")
  plot(ra,real(c1),":g")
  legend( "real(c0)", "imag(c0)", "imag(c1)", "real(c1)", "location", "southwest" )
  title "(b) Hadamard result with resonance error"
  xlabel "a/d"
  ylabel 'real, imag'
  print -dpng pulse2.png

elseif( ex == 3) % example 3, Hadamard with pulse timing error

  %  at t=0: c0(0)=c1(0)=1/sqrt(2)
  %  m=-1, p=pi/2, a=0, b=d
  %  pulse is applied until time t, A0=d*t=pi/4 with no timing error
  %   c0(t) = ( cos(A) + sin(A)) / sqrt(2)
  %   c1(t) = (-cos(A) + sin(A)) / sqrt(2)
  %
  x = [0:50]/25;
  A = x*(pi/4); % 0 ... pi/2
  r2 = sqrt(2)/2;
  C = cos(A);
  S = sin(A);
  c0 = r2*(S+C);
  c1 = r2*(S-C);
  hold off
  plot(x,c0,"-k")
  hold on
  grid on
  plot(x,c1,":g")
  legend( "c0", "c1", "location", "southeast" )
  title "(c) Hadamard result with pulse timing error"
  xlabel("A/A0"); %  A_0 -> "interpreter", "tex");
  ylabel 'c0, c1'
  print -dpng pulse3.png

elseif( ex == 4) % example 4, Hadamard with phase shift error

  %  at t=0: c0(0)=c1(0)=1/sqrt(2)
  %  m=-1, a=0, b=d, and pulse is applied until time t2 = pi/(2*Omega) = pi/(4*b), A = pi/4
  %  p0 = pi/2 with no phase shift error
  %   c0(t2) = ( 1 - i*exp( i*p)) / 2 = (sin(p) + 1 - i*cos(p)) / 2
  %   c1(t2) = (-1 + i*exp(-i*p)) / 2 = (sin(p) - 1 + i*cos(p)) / 2
  %
  x = [0:50]/25;
  p = x*(pi/2); % 0 ... pi
  C = 1 - i*cos(p);
  S = sin(p);
  c0 = (S + C)/2;
  c1 = (S - C)/2;
  % note: imag(c0) == -imag(c1)
  hold off
  plot(x,real(c0),"-k")
  hold on
  grid on
  plot(x,imag(c0),"--r")
  plot(x,imag(c1),"-.b")
  plot(x,real(c1),":g")
  legend( "real(c0)", "imag(c0)", "imag(c1)", "real(c1)", "location", "southeast" )
  title "(d) Hadamard result with phase shift error"
  xlabel("p/p0"); % "interpreter", "tex");
  ylabel 'real, imag'
  print -dpng pulse4.png

else

  disp( "set ex = 1 or 2 or 3 or 4")

endif