(* File: compare-angles by Hongyou Wu on 2006/08/27 *)
(*                                                  *)
(* Input: aaZZ, bbZZ, efZZ[t_], defZZ[t_]           *)
(*                                                  *)
(* Optional input: inpZZ                            *)
(*                                                  *)
(* Output: npointsZZ, changeZZ                      *)


Clear[iZZ,jZZ,npointsZZ,stepZZ,valueZZ,dvalueZZ,aZZ,alphaZZ,betaZZ,oaZZ,isZZ,tsZZ,
      thetaZZ,oocZZ,ocZZ,changeZZ,ovZZ,odvZZ];


shiftZZ[iZZ_,jZZ_]:=

If[Abs[jZZ-2Pi-iZZ]<Min[ Abs[jZZ-iZZ], Abs[jZZ+2Pi-iZZ] ], -2Pi,

If[Abs[jZZ-iZZ]<=Min[ Abs[jZZ-2Pi-iZZ], Abs[jZZ+2Pi-iZZ] ], 0, 2Pi] ];


(* inp: initial number of points     oa: old angle   *)
(* is: initial shift                 ts: total shift *)
(* ooc: old old change               oc: old change  *)
(* change: chnage of theta on [a,b]                  *)

If[ValueQ[inpZZ]==False || inpZZ<9, inpZZ=9];

npointsZZ=inpZZ; stepZZ=(bbZZ-aaZZ)/(npointsZZ-1);

Do[{valueZZ[jZZ]=efZZ[aaZZ+(jZZ-1)*stepZZ]//N;

    dvalueZZ[jZZ]=defZZ[aaZZ+(jZZ-1)*stepZZ]//N; },

   {jZZ,1,npointsZZ} ];

aZZ=alphaZZ=Arg[ defZZ[aaZZ]+efZZ[aaZZ]*I ];

While[alphaZZ<0, alphaZZ=alphaZZ+Pi]; While[alphaZZ>=Pi, alphaZZ=alphaZZ-Pi];

betaZZ=Arg[ defZZ[bbZZ]+efZZ[bbZZ]*I ];

While[betaZZ<=0, betaZZ=betaZZ+Pi]; While[betaZZ>Pi, betaZZ=betaZZ-Pi];

oaZZ=alphaZZ; isZZ=tsZZ=Round[ (alphaZZ-aZZ)/Pi ]*Pi;

Do[{thetaZZ=Arg[ dvalueZZ[jZZ]+valueZZ[jZZ]*I ]+tsZZ;

    aZZ=shiftZZ[oaZZ,thetaZZ]; oaZZ=thetaZZ=thetaZZ+aZZ; tsZZ=tsZZ+aZZ; },

   {jZZ,2,npointsZZ} ];

oocZZ=ocZZ=-Infinity; changeZZ=(thetaZZ-betaZZ)/Pi; changeZZ=Round[changeZZ];


(* ov: old value  odv: old derivative value *)

iZZ=2;

While[iZZ<100 && (oocZZ!=changeZZ || ocZZ!=changeZZ),

Do[{ovZZ[jZZ]=valueZZ[jZZ]; odvZZ[jZZ]=dvalueZZ[jZZ];}, {jZZ,1,npointsZZ}];

stepZZ=stepZZ/2;

Do[{valueZZ[2*jZZ-1]=ovZZ[jZZ]; dvalueZZ[2*jZZ-1]=odvZZ[jZZ];

    valueZZ[2*jZZ]=efZZ[aaZZ+(2*jZZ-1)*stepZZ]//N;

    dvalueZZ[2*jZZ]=defZZ[aaZZ+(2*jZZ-1)*stepZZ]//N; },

   {jZZ,1,npointsZZ-1} ];

aZZ=2*npointsZZ-1; valueZZ[aZZ]=ovZZ[npointsZZ]; dvalueZZ[aZZ]=odvZZ[npointsZZ];

npointsZZ=aZZ; oaZZ=alphaZZ; tsZZ=isZZ;

Do[{thetaZZ=Arg[ dvalueZZ[jZZ]+valueZZ[jZZ]*I ]+tsZZ;

    aZZ=shiftZZ[oaZZ,thetaZZ]; oaZZ=thetaZZ=thetaZZ+aZZ; tsZZ=tsZZ+aZZ; },

   {jZZ,2,npointsZZ} ];

oocZZ=ocZZ; ocZZ=changeZZ; changeZZ=(thetaZZ-betaZZ)/Pi; changeZZ=Round[changeZZ];

iZZ=iZZ+1;

];


If[oocZZ!=changeZZ || ocZZ!=changeZZ, Print["Failed to compare angles!"]];


Clear[aaZZ,bbZZ,efZZ,defZZ,inpZZ,iZZ,jZZ,stepZZ,valueZZ,dvalueZZ,aZZ,alphaZZ,
      betaZZ,oaZZ,isZZ,tsZZ,thetaZZ,oocZZ,ocZZ,ovZZ,odvZZ];