(* File: compute-indices-pc by Hongyou Wu on 2006/08/27                 *)
(*                                                                      *)
(* Input: npieces, a[0], a[np], nterms, ev[i], mt[i]                    *)
(*                                                                      *)
(* Optional input: psiYY, inpYY, czYY                                   *)
(*                                                                      *)
(* Output: print-outs, beta0YY, beta0KYY, beta1YY, beta1KYY, indexYY[i] *)


Clear[sYY,mAYY,mBYY,mEYY,gammaYY,mKYY,iYY,jYY,aYY,bYY,cYY,dYY,zYY,beta0KYY,
      beta1KYY,\[Lambda],separatedYY,beta0YY,beta1YY,smallYY,betaYY,efZZ,defZZ,
      inpZZ,indexYY];


If[nterms>0,

(* Check self-adjointness of BC *)

Do[sYY[iYY,jYY]=mAB[[iYY,jYY]], {iYY,1,2}, {jYY,1,4}];

If[sYY[1,1]*sYY[2,3]-sYY[1,3]*sYY[2,1]==0 &&

   sYY[1,1]*sYY[2,4]-sYY[1,4]*sYY[2,1]==0 &&

   sYY[1,2]*sYY[2,3]-sYY[1,3]*sYY[2,2]==0 &&

   sYY[1,2]*sYY[2,4]-sYY[1,4]*sYY[2,2]==0,

   Print["The boundary condition is raelly NOT self-adjoint!"], {

mAYY={{sYY[1,1], sYY[1,2]},
      {sYY[2,1], sYY[2,2]}};

mBYY={{sYY[1,3], sYY[1,4]},
      {sYY[2,3], sYY[2,4]}};

mEYY={{0,-1},
      {1, 0}};

If[mAYY.mEYY.Conjugate[Transpose[mAYY]]!=mBYY.mEYY.Conjugate[Transpose[mBYY]],

   Print["The boundary condition is NOT self-adjoint!"], {


(* normalize BC: real if separated, [e^(I gamma) K | -I] if coupled; & *)
(* compute beta0K & beta1K if coupled                                  *)

If[sYY[1,1]*sYY[2,2]-sYY[1,2]*sYY[2,1]==0,

   {separatedYY=1;

    If[Abs[ sYY[1,1] ]^2+Abs[ sYY[1,2] ]^2<Abs[ sYY[2,1] ]^2+Abs[ sYY[2,2] ]^2,

       Do[{aYY=sYY[1,jYY]; sYY[1,jYY]=sYY[2,jYY]; sYY[2,jYY]=aYY;}, {jYY,1,4}] ];

    If[Abs[ sYY[1,1] ]>=Abs[ sYY[1,2] ],

       {aYY=Arg[ sYY[1,1] ]; sYY[1,1]=Abs[ sYY[1,1] ];

        sYY[1,2]=Simplify[ Exp[-I*aYY]*sYY[1,2] ]; cYY=sYY[2,1]/sYY[1,1]; },

       {aYY=Arg[ sYY[1,2] ]; sYY[1,1]=Simplify[ Exp[-I*aYY]*sYY[1,1] ];

        sYY[1,2]=Abs[ sYY[1,2] ]; cYY=sYY[2,2]/sYY[1,2]; } ];

    sYY[2,3]=sYY[2,3]-cYY*sYY[1,3]; sYY[2,4]=sYY[2,4]-cYY*sYY[1,4];

    If[Abs[ sYY[2,3] ]>=Abs[ sYY[2,4] ],

       {aYY=Arg[ sYY[2,3] ]; sYY[2,3]=Abs[ sYY[2,3] ];

        sYY[2,4]=Simplify[ Exp[-I*aYY]*sYY[2,4] ]; },

       {aYY=Arg[ sYY[2,4] ]; sYY[2,3]=Simplify[ Exp[-I*aYY]*sYY[2,3] ];

        sYY[2,4]=Abs[ sYY[2,4] ]; } ];

    mAB={{sYY[1,1], sYY[1,2],        0,        0},
         {       0,        0, sYY[2,3], sYY[2,4]}}; },

   {separatedYY=-1; mAB=-Inverse[mBYY].mAB;

    Do[sYY[iYY,jYY]=mAB[[iYY,jYY]], {iYY,1,2}, {jYY,1,4}];

    aYY=0;

    Do[{bYY=sYY[iYY,jYY]; bYY=Abs[bYY]; If[bYY>aYY, aYY=bYY; cYY=iYY; dYY=jYY;];},

       {iYY,1,2}, {jYY,1,2} ];

    gammaYY=sYY[cYY,dYY]; gammaYY=Arg[gammaYY];

    Do[sYY[iYY,jYY]=Simplify[ Exp[-I*gammaYY]*sYY[iYY,jYY] ],

       {iYY,1,2}, {jYY,1,2} ];

    mKYY={{sYY[1,1], sYY[1,2]},
          {sYY[2,1], sYY[2,2]}};

    mAB={{Exp[I*gammaYY]*sYY[1,1], Exp[I*gammaYY]*sYY[1,2], -1,  0},
         {Exp[I*gammaYY]*sYY[2,1], Exp[I*gammaYY]*sYY[2,2],  0, -1}};

    Do[{aYY=mKYY[[1,jYY]]; bYY=mKYY[[2,jYY]];

        If[Abs[aYY]<Abs[bYY], cYY=Re[ArcTan[aYY/bYY]], cYY=Re[ArcCot[bYY/aYY]]];

        While[cYY<=0, cYY=cYY+Pi]; While[cYY>Pi, cYY=cYY-Pi]; sYY[jYY,1]=cYY; },

       {jYY,1,2} ];

    beta0KYY=sYY[2,1]; beta1KYY=sYY[1,1]; } ];


(* Choose an eigenvalue with smallest absolute value, & evaluate Phi *)
(*                                                                   *)
(* psi: partial sequence index                                       *)

If[ValueQ[psiYY]==False, psiYY=1; aYY=ev[1]; aYY=Abs[aYY];

   Do[{bYY=ev[iYY]; bYY=Abs[bYY]; If[bYY<aYY, psiYY=iYY; aYY=bYY;];},

      {iYY,2,nterms} ]; ];

bYY=a[npieces]; zYY=\[Lambda]=ev[psiYY]; <<Phi-matrix-pc;

Do[sYY[iYY,jYY]=\[Phi][iYY,jYY,bYY,zYY], {iYY,1,2}, {jYY,1,2}];


(* If BC is separated, compute Cos[alpha] & Sin[alpha];                  *)
(* if not, calculate beta0 & beta1, find a corresponding separated BC, & *)
(*         compute Cos[alpha] & Sin[alpha]                               *)

If[separatedYY>0, {aYY=mAB[[1,1]]; bYY=-mAB[[1,2]];},

   {If[mt[psiYY]>1, {beta0YY=beta0KYY; beta1YY=beta1KYY;},

       {Do[{aYY=sYY[1,jYY]; bYY=sYY[2,jYY];

            If[Abs[aYY]<Abs[bYY], cYY=Re[ArcTan[aYY/bYY]],

               cYY=Re[ArcCot[bYY/aYY]] ];

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

            If[jYY>3/2, beta0YY=cYY, beta1YY=cYY]; },

           {jYY,1,2} ];

        If[beta0YY<0.0001 || beta0YY>0.9999*Pi, aYY=beta0YY//N;

           Print["Warning: \[Beta]0=",aYY," is close to 0 or ",Pi,","];

           Print["use a different eigenvalue to verify the indices."]; ]; } ];

    If[beta0YY<=Pi/2, {smallYY=1; betaYY=(beta0YY+Pi)/2;},

       {smallYY=-1; betaYY=beta0YY/2;} ];

    aYY=sYY[1,1]*Cos[betaYY]-sYY[2,1]*Sin[betaYY];

    bYY=-sYY[1,2]*Cos[betaYY]+sYY[2,2]*Sin[betaYY]; } ];

cYY=aYY*aYY+bYY*bYY; aYY=aYY/cYY; bYY=bYY/cYY;


(* Form eigenfunction for separated BC, & count its zeros                     *)
(*                                                                            *)
(* ef: eigenfunction              def: derivative of eigenfunction            *)
(* inp: initial number of points  cz: 1 if zeros are to be counted, -1 if not *)

efZZ[tYY_]=Re[ bYY*\[Phi][1,1,tYY,zYY]+aYY*\[Phi][1,2,tYY,zYY] ];

defZZ[tYY_]=Re[ bYY*\[Phi][2,1,tYY,zYY]+aYY*\[Phi][2,2,tYY,zYY] ];

aaZZ=a[0]; bbZZ=a[npieces]; If[ValueQ[inpYY]==True, inpZZ=inpYY];

<<compare-angles; cYY=npointsZZ;

If[ValueQ[czYY]==False, czYY=-1];

If[czYY>0, efZZ[tYY_]=Re[ bYY*\[Phi][1,1,tYY,zYY]+aYY*\[Phi][1,2,tYY,zYY] ];

   If[ValueQ[inpYY]==True, inpZZ=inpYY]; <<count-zeros;

   If[nzerosZZ!=changeZZ, cYY=Max[cYY,npointsZZ];

      Print["Failed to determine the number of zeros inside open interval!"];

      Print["Suggested initial number of points: inpYY=",cYY]; ]; ];


(* If double, determine smaller index;   *)
(* if coupled & simple, determine index; *)
(* if separated, determine index.        *)

If[mt[psiYY]>1, If[smallYY<0, changeZZ=changeZZ-1],

   If[separatedYY<0,

      aYY=Min[beta0KYY,beta0YY]; bYY=Max[beta0KYY,beta0YY]; cYY=bYY-aYY;

      If[cYY<10^(-3) && (aYY<=beta1KYY<=bYY || aYY<=beta1YY<=bYY),

         {aYY=beta0YY//N; bYY=beta0KYY//N;

          Print["Warning: \[Beta]0=",aYY," is close to \[Beta]0K=",bYY,","];

          Print["use a different eigenvalue to verify the indices."];

          If[smallYY<0 && beta0YY<beta0KYY, changeZZ=changeZZ-1];

          If[smallYY>0 && beta0KYY<beta0YY, changeZZ=changeZZ+1]; },

         {If[smallYY<0 &&

             ((cYY>=10^(-3) && beta0YY<beta0KYY) ||

              (cYY<10^(-3) &&

               (beta1KYY>beta1YY>beta0KYY || beta1YY>beta0YY>beta1KYY ||

                beta0YY>beta1KYY>beta1YY ) ) ), changeZZ=changeZZ-1 ];

          If[smallYY>0 &&

             ((cYY>=10^(-3) && beta0KYY<beta0YY) ||

              (cYY<10^(-3) &&

               (beta1YY>beta1KYY>beta0YY || beta1KYY>beta0YY>beta1YY ||

                beta0YY>beta1YY>beta1KYY ) ) ), changeZZ=changeZZ+1 ]; } ]; ] ];

changeZZ=changeZZ+1; aYY=changeZZ-Sum[mt[iYY],{iYY,1,psiYY-1}];

Do[{If[mt[iYY]<2, Print["index of ",ev[iYY]//N,": ",aYY],

       Print["indices of ",ev[iYY]//N,": from ",aYY," to ",aYY+mt[iYY]-1] ];

    indexYY[iYY]=aYY; aYY=aYY+mt[iYY]; },

   {iYY,1,nterms} ];

} ]; } ]; ];


Clear[psiYY,inpYY,czYY,sYY,mAYY,mBYY,mEYY,gammaYY,mKYY,iYY,jYY,aYY,bYY,cYY,dYY,zYY,
      \[Lambda],separatedYY,smallYY,betaYY,changeZZ,npointsZZ,nzerosZZ];