* This file contains the followings
*
***  WQCfunc, WmCfunc,  VmCfunc
***  WGCfunc, WmuCfunc, VmuCfunc
***  WGdCfunc, VGdCfunc
***  WGNKCfunc, VLog, LogNK
***  WQintCfunc, WQintVCfunc, VQintCfunc
***  WGintCfunc, WGintVCfunc, VGintCfunc 
***  EXPINT
***  StrmLnPlt, traject, TOKOSN 
***  Use  gr7WB   gr7WBd   gr7WBBW   gr7WBdBW   to compile Main and Those PG's
***                                 Fixed on 2015/11/ 2016/07/
************************************************************************
      complex function WQCfunc(z,zeta,K,ThetaCut,WQS)
*
*     Complex Potential of Wave Source WQ
*
*      ThetaCut ( Angle of Cut of Log(z-zeta) ) must be POSITIVE
*
*      ThetaCut - 2*pi  <    Argument of Log(z-zeta)      =<  ThetaCut
*                      .lt.                              .le.
*             - 3/2 PI  < Argument of Log(z-zetaBar) & E1 =<  PI/2
*
      real    K
      complex i/(0,1)/
      complex z,zeta,zetaBar,miKZ,E1
      complex LogZ,LogZBar,EmiKZ,SkappaC,SkappaS,WQS
      pi=acos(0.0)*2
      pi2=2*pi
      pih=pi/2
*
      LogZ=log(z-zeta)
      Alog=real(LogZ)
      ArgZ=imag(LogZ)
      if(ArgZ.le.ThetaCut-pi2) then
         ArgZ=ArgZ+pi2
      elseif(ArgZ.gt.ThetaCut) then
         ArgZ=ArgZ-pi2
      endif
      LogZ=cmplx(ALog,ArgZ)
*
      zetaBar=conjg(zeta)
      LogZBar=log(z-zetaBar)
      Alog=real(LogZBar)
      ArgZ=imag(LogZBar)
      if(ArgZ.gt.pih) ArgZ=ArgZ-pi2
      LogZBar=cmplx(ALog,ArgZ)
*
      miKZ=-i*K*(z-zetaBar)
      EmiKZ=exp(miKZ)
      call expint(miKZ,E1)
*
      xsign=real(z-zeta)
      if(xsign.gt.0) then
         SkappaC=EmiKZ*(E1-pi*i)
      else
         SkappaC=EmiKZ*(E1+pi*i)
      endif
      SkappaS=-pi*EmiKZ
      WQCfunc=LogZ-LogZBar-2*SkappaC
      WQS=-2*SkappaS
      end
************************************************************************
      complex function WmCfunc(z,zeta,K,WmS)
*
*     Complex Potentail of Wave x-Doublet Wm
*               = -d/dz(WQ)
*
      real    K
      complex i/(0,1)/
      complex z,zeta,zetaBar,miKZ,E1
      complex InvZ,InvZBar,EmiKZ,SkappaC,SkappaS,WmS
      pi=acos(0.0)*2
*
      InvZ=1/(z-zeta)
      zetaBar=conjg(zeta)
      InvZBar=1/(z-zetaBar)
      miKZ=-i*K*(z-zetaBar)
      EmiKZ=exp(miKZ)
      call expint(miKZ,E1)
      xsign=real(z-zeta)
      if(xsign.gt.0) then
         SkappaC=EmiKZ*(E1-pi*i)
      else
         SkappaC=EmiKZ*(E1+pi*i)
      endif
      SkappaS=-pi*EmiKZ
      WmCfunc=-InvZ-InvZBar-2*i*K*SkappaC
      WmS=-2*i*K*SkappaS
      end
************************************************************************
      complex function VmCfunc(z,zeta,K,VmS)
*
*     Complex Velocity Induced by Wave x-Doublet
*               = d/dz(Wm)
*
      real    K
      complex i/(0,1)/
      complex z,zeta,zetaBar,miKZ,E1
      complex InvZ2,InvZBar,InvZBar2,EmiKZ,SkappaC,SkappaS,VmS
      pi=acos(0.0)*2
*
      InvZ2=1/(z-zeta)**2
      zetaBar=conjg(zeta)
      InvZBar=1/(z-zetaBar)
      InvZBar2=1/(z-zetaBar)**2
      miKZ=-i*K*(z-zetaBar)
      EmiKZ=exp(miKZ)
      call expint(miKZ,E1)
      xsign=real(z-zeta)
      if(xsign.gt.0) then
         SkappaC=EmiKZ*(E1-pi*i)
      else
         SkappaC=EmiKZ*(E1+pi*i)
      endif
      SkappaS=-pi*EmiKZ
      VmCfunc=InvZ2+InvZBar2+2*i*K*InvZBar-2*K**2*SkappaC
      VmS=-2*K**2*SkappaS
      end


***********************************************************************
      complex function WGCfunc(z,zeta,K,ThetaCut,WGS)
*
*     Complex Potential of Wave Vortex WG
*
*       ThetaCut ( Angle of Cut of Log(z-zeta) ) must be POSITIVE
*
*      ThetaCut - 2*pi  <    Argument of Log(z-zeta)      =<  ThetaCut
*                      .lt.                              .le.
*             - 3/2 PI  < Argument of Log(z-zetaBar) & E1 =<  PI/2
*
      real    K
      complex i/(0,1)/
      complex z,zeta,zetaBar,miKZ,E1
      complex LogZ,LogZBar,EmiKZ,SkappaC,SkappaS,WGS
*
      pi=acos(0.0)*2
      pi2=2*pi
      pih=pi/2
*
      LogZ=log(z-zeta)
      Alog=real(LogZ)
      ArgZ=imag(LogZ)
      if(ArgZ.le.ThetaCut-pi2) then
         ArgZ=ArgZ+pi2
      elseif(ArgZ.gt.ThetaCut) then
         ArgZ=ArgZ-pi2
      endif
      LogZ=cmplx(ALog,ArgZ)
*
      zetaBar=conjg(zeta)
      LogZBar=log(z-zetaBar)
      Alog=real(LogZBar)
      ArgZ=imag(LogZBar)
      if(ArgZ.gt.pih) ArgZ=ArgZ-pi2
      LogZBar=cmplx(ALog,ArgZ)
*
      miKZ=-i*K*(z-zetaBar)
      EmiKZ=exp(miKZ)
      call expint(miKZ,E1)
*
      xsign=real(z-zeta)
      if(xsign.gt.0) then
         SkappaC=EmiKZ*(E1-pi*i)
      else
         SkappaC=EmiKZ*(E1+pi*i)
      endif
      SkappaS=-pi*EmiKZ
      WGCfunc=-i*(LogZ+LogZBar)-2*i*SkappaC
      WGS=-2*i*SkappaS
      end
************************************************************************
      complex function WmuCfunc(z,zeta,K,WmuS)
*
*     Complex Potential of Wave -y-Doublet Wmu
*              = -d/dz(WG)
*
      real    K
      complex i/(0,1)/
      complex z,zeta,zetaBar,miKZ,E1
      complex InvZ,InvZBar,EmiKZ,SkappaC,SkappaS,WmuS
      pi=acos(0.0)*2
*
      InvZ=1/(z-zeta)
      zetaBar=conjg(zeta)
      InvZBar=1/(z-zetaBar)
      miKZ=-i*K*(z-zetaBar)
      EmiKZ=exp(miKZ)
      call expint(miKZ,E1)
      xsign=real(z-zeta)
      if(xsign.gt.0) then
         SkappaC=EmiKZ*(E1-pi*i)
      else
         SkappaC=EmiKZ*(E1+pi*i)
      endif
      SkappaS=-pi*EmiKZ
      WmuCfunc=i*InvZ-i*InvZBar+2*K*SkappaC
      WmuS=2*K*SkappaS
      end
************************************************************************
      complex function VmuCfunc(z,zeta,K,VmuS)
*
*     Complex Velocity Induced by Wave -y-Doublet
*             = d/dz(Wmu)
*
      real    K
      complex i/(0,1)/
      complex z,zeta,zetaBar,miKZ,E1
      complex InvZ2,InvZBar,InvZBar2,EmiKZ,SkappaC,SkappaS,VmuS
      pi=acos(0.0)*2
*
      InvZ2=1/(z-zeta)**2
      zetaBar=conjg(zeta)
      InvZBar=1/(z-zetaBar)
      InvZBar2=1/(z-zetaBar)**2
      miKZ=-i*K*(z-zetaBar)
      EmiKZ=exp(miKZ)
      call expint(miKZ,E1)
      xsign=real(z-zeta)
      if(xsign.gt.0) then
         SkappaC=EmiKZ*(E1-pi*i)
      else
         SkappaC=EmiKZ*(E1+pi*i)
      endif
      SkappaS=-pi*EmiKZ
      VmuCfunc=-i*InvZ2+i*InvZBar2-2*K*InvZBar-2*i*K**2*SkappaC
      VmuS=-2*i*K**2*SkappaS
      end


***********************************************************************
      complex function WGdCfunc(z,zeta2,zeta1,K,WGdS)
*
*     a couple of wave VORTEX WG(a,zeta2) - WG(z,zeta1)
*
*     Cut :  zeta1 -- zeta2
*
      real    K
      complex WG1C,WG1S,WG2C,WG2S
      complex WGCfunc,WGNKCfunc
      complex z,zeta2,zeta1,WGdS
*
      pi=acos(0.0)*2
      thetacut=pi/2
*  for zeta1
      WG1C=WGNKCfunc(z,zeta1,K,  zeta2, WG1S)
*  for zeta2
      WG2C=  WGCfunc(z,zeta2,K,thetacut,WG2S)
*  Difference
      WGdCfunc=WG2C-WG1C
      WGdS    =WG2S-WG1S
      end
***********************************************************************
      complex function VGdCfunc(z,zeta2,zeta1,K,VGdS)
*
*     Complex velocity induced by a couple of wave VORTEX
*             = d/dz(WG2-WG1)
*
      real    K
      complex VG1C,VG1S,VG2C,VG2S
      complex WmuCfunc,WmuS
      complex z,zeta2,zeta1,VGdS
*
*  for zeta1
      VG1C=-WmuCfunc(z,zeta1,K,WmuS)
      VG1S=-WmuS
*  for zeta2
      VG2C=-WmuCfunc(z,zeta2,K,WmuS)
      VG2S=-WmuS
*  Difference
      VGdCfunc=VG2C-VG1C
      VGdS    =VG2S-VG1S
      end


***********************************************************************
      complex function WGNKCfunc(z,zeta,K,zetaM,WGS)
*                        WB
*
*     Wave Vortex with cut of Log like in N-K problem
*
*     Cut of log(z-zeta) is set as
*        zeta -- zetaM -- (xsiM,inf)
*                                  
*     Cut of log(z-zetaBar) as
*        zetaBar -- (xsiBar,inf)
*
*     Cut of E1 as
*        zeatBar -- (xsiBar,inf)
*
*
      real    K
      complex i/(0,1)/
      complex z,zeta,zetaM,zetaBar,miKZ,E1,LogNK
      complex LogZ,LogZBar,EmiKZ,SkappaC,SkappaS,WGS
*
      pi=acos(0.0)*2
      pi2=2*pi
      pih=pi/2
*          _____
      LogZ=LogNK(z,zeta,zetaM)
*          ~~~~~
      zetaBar=conjg(zeta)
      LogZBar=log(z-zetaBar)
      Alog=real(LogZBar)
      ArgZ=imag(LogZBar)
      if(ArgZ.gt.pih) ArgZ=ArgZ-pi2
      LogZBar=cmplx(ALog,ArgZ)
*
      miKZ=-i*K*(z-zetaBar)
      EmiKZ=exp(miKZ)
      call expint(miKZ,E1)
*
      xsign=real(z-zeta)
      if(xsign.gt.0) then
         SkappaC=EmiKZ*(E1-pi*i)
      else
         SkappaC=EmiKZ*(E1+pi*i)
      endif
      SkappaS=-pi*EmiKZ
      WGNKCfunc=-i*(LogZ+LogZBar)-2*i*SkappaC
      WGS=-2*i*SkappaS
      end
***************************************************************************
      complex function VLog(z)
*
*     Complex Logarithmic Function for Img(z) > 0
*     Cut : z -- (x,inf)
*
      complex z,pi2i
*
      pi=acos(0.0)*2
      pi2i=cmplx(0,pi*2)
      pih=pi/2
*
      VLog=log(z)
      if(imag(VLog).gt.pih) VLog=VLog-pi2i
      end
***************************************************************************
      complex function LogNK(z,zeta,zetaM)
*
*     Complex Logarithmic Function for Neumann-Kelvin Problem
*     Cut : zeta -- zetaM -- (xsiM,inf)
*     ( See Section 3.3.4 )
*
      complex z,zeta,zetaM,LogZ
*
      pi=acos(0.0)*2
      pi2=2*pi
*
      LogZ=log(z-zeta)
      Alog=real(LogZ)
      ArgZ=imag(LogZ)
*
      x=real(z)
      y=imag(z)
      xsi=real(zeta)
      eta=imag(zeta)
      xsiM=real(zetaM)
      etaM=imag(zetaM)
* Case I
      if(xsi.ge.xsiM .and. eta.ge.etaM) then
         if(x.lt.xsiM .and. y.ge.eta) then
            ArgZ=ArgZ-pi2
         elseif(x.ge.xsiM .and. y.lt.eta .and. 
     *          (y-etaM)*(xsi-xsiM).ge.(eta-etaM)*(x-xsiM)) then
            ArgZ=ArgZ+pi2
         endif
* Case II
      elseif(xsi.lt.xsiM .and. eta.ge.etaM) then
         if(x.lt.xsiM .and.
     *  (y.ge.eta .or.(y-eta)*(xsiM-xsi).ge.(etaM-eta)*(x-xsi))) then
            ArgZ=ArgZ-pi2
         endif            
* Case III
      elseif(xsi.lt.xsiM .and. eta.lt.etaM) then
         if(x.lt.xsiM .and. y.ge.eta .and. 
     *  (y-eta)*(xsiM-xsi).ge.(etaM-eta)*(x-xsi)) then
            ArgZ=ArgZ-pi2
         endif            
* Case IV
      else
         if(y.ge.eta .and. 
     *  (x.lt.xsiM .or.(y-etaM)*(xsi-xsiM).lt.(eta-etaM)*(x-xsiM))) then
            ArgZ=ArgZ-pi2
         endif            
      endif
      LogNK=cmplx(ALog,ArgZ)
      end


************************************************************************
      complex function WQintCfunc(z,zeta2,zeta1,K,zetaM,WQintS)
*                      WQintWB_C,S
*                      WQintNK_C,S
*     
*     zeta1, zeta2 are exchangable
*
*     Cut of log(z-zeta) is set as
*        zeta1 -- zetaM -- zeta2
*
*     Cut of log(z-zetaBar) and E1 as
*        zetaBar -- (xsi,inf)
*
      real    K
      complex i/(0,1)/
      complex WQintC1,WQintC2
      complex WQintS1,WQintS2,WQintS
      complex z,zeta1,zeta2,zetaM
      complex zdif,zdifBar,dzeta
      complex LogNK,VLog
      complex LogZ,LogZBar
      complex EiG,miKz,EmiKz,E1,SKappaC
      pi=acos(0.0)*2
      pi2=pi*2
*
      dzeta=zeta2-zeta1
      EiG=dzeta/abs(dzeta)
*
*     For zeta2
*
      LogZ=LogNK(z,zeta2,zetaM)
**    LogZ=VLog(z-zeta2)
      zdif   =z-zeta2
      zdifBar=z-conjg(zeta2)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta2).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WQintC2=-conjg(EiG)*zdif   *(LogZ-1)
     *        +      EiG *zdifBar*(LogZBar-1)
     *        +i*2/K*EiG         *(LogZBar+SKappaC)
      WQintS2=-i*pi2/K*EiG*EmiKz
*
*     For zeta1
*
      LogZ=LogNK(z,zeta1,zetaM)
**    LogZ=LogNK(z,zeta1,zeta2)
      zdif=z-zeta1
*
      zdifBar=z-conjg(zeta1)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta1).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WQintC1=-conjg(EiG)*zdif   *(LogZ-1)
     *        +      EiG *zdifBar*(LogZBar-1)
     *        +i*2/K*EiG         *(LogZBar+SKappaC)
      WQintS1=-i*pi2/K*EiG*EmiKz
*
      WQintCfunc=WQintC2-WQintC1
      WQintS    =WQintS2-WQintS1
      end



************************************************************************
      complex function WQintVCfunc(z,zeta2,zeta1,K,WQintS)
*                      WQintV_C,S
*     
*     zeta1, zeta2 are exchangable
*
*     Cut of log(z-zeta) is set as
*        zeta1 -- zeta2 -- (xsi2,inf)
*
*     Cut of log(z-zetaBar) and E1 as
*        zetaBar -- (xsiBar,inf)
*
      real    K
      complex i/(0,1)/
      complex WQintC1,WQintC2
      complex WQintS1,WQintS2,WQintS
      complex z,zeta1,zeta2
      complex zdif,zdifBar,dzeta
      complex LogNK,VLog
      complex LogZ,LogZBar
      complex EiG,miKz,EmiKz,E1,SKappaC
      pi=acos(0.0)*2
      pi2=pi*2
*
      dzeta=zeta2-zeta1
      EiG=dzeta/abs(dzeta)
*
*     For zeta2
*
**    LogZ=LogNK(z,zeta2,zetaM)
**    LogZ=VLog(z-zeta2)
      zdif   =z-zeta2
      LogZ=VLog(zdif)
      zdifBar=z-conjg(zeta2)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta2).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WQintC2=-conjg(EiG)*zdif   *(LogZ-1)
     *        +      EiG *zdifBar*(LogZBar-1)
     *        +i*2/K*EiG         *(LogZBar+SKappaC)
      WQintS2=-i*pi2/K*EiG*EmiKz
*
*     For zeta1
*
**    LogZ=LogNK(z,zeta1,zetaM)
      LogZ=LogNK(z,zeta1,zeta2)
      zdif=z-zeta1
*
      zdifBar=z-conjg(zeta1)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta1).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WQintC1=-conjg(EiG)*zdif   *(LogZ-1)
     *        +      EiG *zdifBar*(LogZBar-1)
     *        +i*2/K*EiG         *(LogZBar+SKappaC)
      WQintS1=-i*pi2/K*EiG*EmiKz
*
      WQintVCfunc=WQintC2-WQintC1
      WQintS     =WQintS2-WQintS1
      end



************************************************************************
      complex function VQintCfunc(z,zeta2,zeta1,K,VQintS)
*                      VQintWB_C,S
*
*     zeta1, zeta2 are exchangable
*     
*     Cut of log(z-zeta) is set as
*        zeta1 -- zeta2 -- (xsi2,inf)
*
*     Cut of log(z-zetaBar) and E1 as
*        zetaBar -- (xsiBar,inf)
*
      real    K
      complex i/(0,1)/
      complex VQintC1,VQintC2
      complex VQintS1,VQintS2,VQintS
      complex z,zeta1,zeta2
      complex zdifBar,dzeta
      complex LogNK,VLog
      complex LogZ,LogZBar
      complex EiG,miKz,EmiKz,E1,SKappaC
      pi=acos(0.0)*2
      pi2=pi*2
*
      dzeta=zeta2-zeta1
      EiG=dzeta/abs(dzeta)
*
*     For zeta2
*
      LogZ=VLog(z-zeta2)
*
      zdifBar=z-conjg(zeta2)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta2).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      VQintC2=-conjg(EiG)*LogZ + EiG*(LogZBar+2*SKappaC)
      VQintS2=-pi2*EiG*EmiKz
*
*     For zeta1
*
      LogZ=LogNK(z,zeta1,zeta2)
*
      zdifBar=z-conjg(zeta1)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta1).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      VQintC1=-conjg(EiG)*LogZ + EiG*(LogZBar+2*SKappaC)
      VQintS1=-pi2*EiG*EmiKz
*
      VQintCfunc=VQintC2-VQintC1
      VQintS    =VQintS2-VQintS1
      end



************************************************************************
      complex function WGintCfunc(z,zeta2,zeta1,K,zetaM,WGintS)
*                      WGintWB_C,S
*                      WGintNK_C,S
*
*     zeta1, zeta2 are exchangable 
*
*     Cut of log(z-zeta) is set as
*        zeta1 -- zetaM -- zeta2
*
*     Cut of log(z-zetaBar) and E1 as
*        zetaBar -- (xsi,inf)
*
      real    K
      complex i/(0,1)/
      complex WGintC1,WGintC2
      complex WGintS1,WGintS2,WGintS
      complex z,zeta1,zeta2,zetaM
      complex zdif,zdifBar,dzeta
      complex LogNK,VLog
      complex LogZ,LogZBar
      complex EiG,miKz,EmiKz,E1,SKappaC
      pi=acos(0.0)*2
      pi2=pi*2
*
      dzeta=zeta2-zeta1
      EiG=dzeta/abs(dzeta)
*
*     For zeta2
*
      LogZ=LogNK(z,zeta2,zetaM)
**    LogZ=VLog(z-zeta2)
      zdif   =z-zeta2
      zdifBar=z-conjg(zeta2)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta2).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WGintC2= i*conjg(EiG)*zdif   *(LogZ-1)
     *        +i*      EiG *zdifBar*(LogZBar-1)
     *        -2/K*EiG             *(LogZBar+SKappaC)
      WGintS2=pi2/K*EiG*EmiKz
*
*     For zeta1
*
      LogZ=LogNK(z,zeta1,zetaM)
      zdif=z-zeta1
*
      zdifBar=z-conjg(zeta1)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta1).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WGintC1= i*conjg(EiG)*zdif   *(LogZ-1)
     *        +i*      EiG *zdifBar*(LogZBar-1)
     *        -2/K*EiG             *(LogZBar+SKappaC)
      WGintS1=pi2/K*EiG*EmiKz
*
      WGintCfunc=WGintC2-WGintC1
      WGintS    =WGintS2-WGintS1
      end




************************************************************************
      complex function WGintVCfunc(z,zeta2,zeta1,K,WGintS)
*                      WGintV_C,S
*
*     zeta1, zeta2 are exchangable 
*
*     Cut of log(z-zeta) is set as
*        zeta1 -- zeta2 -- (xsi2,inf)
*
*     Cut of log(z-zetaBar) and E1 as
*        zetaBar -- (xsi,inf)
*
      real    K
      complex i/(0,1)/
      complex WGintC1,WGintC2
      complex WGintS1,WGintS2,WGintS
      complex z,zeta1,zeta2
      complex zdif,zdifBar,dzeta
      complex LogNK,VLog
      complex LogZ,LogZBar
      complex EiG,miKz,EmiKz,E1,SKappaC
      pi=acos(0.0)*2
      pi2=pi*2
*
      dzeta=zeta2-zeta1
      EiG=dzeta/abs(dzeta)
*
*     For zeta2
*
**    LogZ=LogNK(z,zeta2,zetaM)
**    LogZ=VLog(z-zeta2)
      zdif   =z-zeta2
      LogZ=VLog(zdif)
      zdifBar=z-conjg(zeta2)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta2).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WGintC2= i*conjg(EiG)*zdif   *(LogZ-1)
     *        +i*      EiG *zdifBar*(LogZBar-1)
     *        -2/K*EiG             *(LogZBar+SKappaC)
      WGintS2=pi2/K*EiG*EmiKz
*
*     For zeta1
*
**    LogZ=LogNK(z,zeta1,zetaM)
      LogZ=LogNK(z,zeta1,zeta2)
      zdif=z-zeta1
*
      zdifBar=z-conjg(zeta1)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta1).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      WGintC1= i*conjg(EiG)*zdif   *(LogZ-1)
     *        +i*      EiG *zdifBar*(LogZBar-1)
     *        -2/K*EiG             *(LogZBar+SKappaC)
      WGintS1=pi2/K*EiG*EmiKz
*
      WGintVCfunc=WGintC2-WGintC1
      WGintS     =WGintS2-WGintS1
      end




************************************************************************
      complex function VGintCfunc(z,zeta2,zeta1,K,VGintS)
*                      VGintWB_C,S
*
*     zeta1, zeta2 are exchangable
*     
*     Cut of log(z-zeta) is set as
*        zeta1 -- zeta2 -- (xsi2,inf)
*
*     Cut of log(z-zetaBar) and E1 as
*        zetaBar -- (xsiBar,inf)
*
      real    K
      complex i/(0,1)/
      complex VGintC1,VGintC2
      complex VGintS1,VGintS2,VGintS
      complex z,zeta1,zeta2
      complex zdifBar,dzeta
      complex LogNK,VLog
      complex LogZ,LogZBar
      complex EiG,miKz,EmiKz,E1,SKappaC
      pi=acos(0.0)*2
      pi2=pi*2
*
      dzeta=zeta2-zeta1
      EiG=dzeta/abs(dzeta)
*
*     For zeta2
*
      LogZ=VLog(z-zeta2)
*
      zdifBar=z-conjg(zeta2)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta2).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      VGintC2=i*conjg(EiG)*LogZ + i*EiG*(LogZBar+2*SKappaC)
      VGintS2=-pi2*i*EiG*EmiKz
*
*     For zeta1
*
      LogZ=LogNK(z,zeta1,zeta2)
*
      zdifBar=z-conjg(zeta1)
      LogZBar=VLog(zdifBar)
      miKz=-i*K*zdifBar
      EmiKz=exp(miKz)
      call expint(miKZ,E1)
      if(real(z)-real(zeta1).gt.0) then
         SKappaC=EmiKz*(E1-pi*i)
      else
         SKappaC=EmiKz*(E1+pi*i)
      endif
      VGintC1=i*conjg(EiG)*LogZ + i*EiG*(LogZBar+2*SKappaC)
      VGintS1=-pi2*i*EiG*EmiKz
*
      VGintCfunc=VGintC2-VGintC1
      VGintS    =VGintS2-VGintS1
      end




************************************************************************
      SUBROUTINE    EXPINT(z,e1)
*
*     EXPONENTIAL INTEGRAL OF COMPLEX ARGUMENT
*     CODED BY KATSUO SUZUKI (BODAI)
*
      IMPLICIT  COMPLEX  (Z)
      REAL INF/9.99999999999999D09/
      COMPLEX  E1,E1GL,EXPE1,LOGZ
      REAL GAMMA/0.5772156649D00/
      REAL PAI/3.14159265359D0/
**    RETURN
*
**    ENTRY  EXPI1 (Z,E1)
**    ENTRY  EXPI2 (Z,E1GL)
**    ENTRY  EXPI3 (Z,EXPE1)
**    ENTRY  EXPI4 (Z,ZEXPE1)
**    ENTRY  EXPIAL (Z,E1,E1GL,EXPE1,ZEXPE1)
*
*          E1 =COMPLEX EXPONENTIAL INTEGRAL OF COMPLEX ARGUMENT Z
*        E1GL =E1 + GAMMA + LOG( Z )
*       EXPE1 = EXP( Z ) * E1
*      ZEXPE1 = Z * EXPE1
*
      R=ABS(Z)
      X=real(Z)
      Y=IMAG(Z)
      IF(R.EQ.0.0)   GO TO 10
      LOGZ=LOG(Z)
      IF(R.GT.25.0)  GO TO 40
      IF(R.GT.4.0.AND.X.GE.0.0)  GO TO 30
      IF(ABS(Y).GT.4.0)  GO TO 30
      GO  TO 20
*
*         R = 0
*
   10 CONTINUE
      E1=DCMPLX(INF,0.0)
      E1GL=(1.0,0.0)
      EXPE1=E1
      ZEXPE1=DCMPLX(GAMMA,0.0)
      WRITE(6,100)
  100 FORMAT('   *** (EXPINT)  Z=(0.0,0.0)  ***')
      RETURN
*
*         SERIES  EXPANSION
*
*              E1GL= - S (-Z)**N/(N*FACT(N))
*                     N=1
*
   20 CONTINUE
      I=1
      ZI=Z
      ZSUM=ZI
   21 I=I+1
      A=-DFLOAT(I-1)/DFLOAT(I*I)
      ZI=A*ZI*Z
      ZSUM=ZSUM+ZI
      RE=ABS(ZSUM)
      AD=ABS(ZI)*1.0E8
      IF(RE.LE.AD)  GO TO 21
      E1GL=ZSUM
      E1=E1GL-GAMMA-LOGZ
      EXPE1=EXP(Z)*E1
      ZEXPE1=Z*EXPE1
      RETURN
*
*         CONTINUED  FRACTION
*
*              EXPE1 = 1/Z+1/1+1/Z+2/1+2/Z+3/1+3/Z+
*
   30 CONTINUE
      I=11
      Z0=(0.0,0.0)
      ZI=Z0
   31 IF(I.EQ.1) GO TO 35
      I=I-1
      AI=DFLOAT(I)
      ZD=Z+ZI
      IF(ZD.NE.0.0) GO TO 32
      ZI=Z0
      GO TO 33
   32 ZI=AI/ZD
   33 CONTINUE
      ZD=1.0+ZI
      IF(ZD.NE.0.0) GO TO 34
      ZI=Z0
      GO TO 31
   34 CONTINUE
      ZI=AI/ZD
      GO TO 31
   35 CONTINUE
      EXPE1=1.0/(Z+ZI)
      E1=EXP(-Z)*EXPE1
      E1GL=E1+GAMMA+LOGZ
      ZEXPE1=Z*EXPE1
      RETURN
*
*         ASYMPTOTIC EXPANSION
*
*              ZEXPE1 = 1-1/Z+FACT(2)/Z**2-FACT(3)/Z**3+
*
   40 CONTINUE
      ZINV=1.0/Z
      I=0
      ZI=(1.0,0.0)
      ZSUM=ZI
   41 I=I+1
      AI=DFLOAT(I)
      ZI=-AI*ZINV*ZI
      ZSUM=ZSUM+ZI
      RE=ABS(ZSUM)
      AD=ABS(ZI)*1.0E8
      IF(RE.LE.AD)  GO TO 41
      IF(Y.EQ.0.0.AND.X.LT.0.0)  ZSUM=ZSUM-Z*EXP(Z)*DCMPLX(0.0,PAI)
      ZEXPE1=ZSUM
      EXPE1=ZEXPE1*ZINV
      IF(-X.GE.200.0) GO TO 42
      E1=EXP(-Z)*EXPE1
      E1GL=E1+GAMMA+LOGZ
      RETURN
   42 E1=INF*EXP(DCMPLX(0.0,-Y))*EXPE1
      E1GL=E1+GAMMA+LOGZ
      RETURN
      END

************************************************************************
      subroutine StrmLnPlt(NxD, Nx,Ny,xpos,ypos, Phi,Psi, L0,z00,ze0,dp)
*
      real    xpos(NxD,*),ypos(NxD,*)
      real    Phi(NxD,*),Psi(NxD,*)
*
      Nxp1=Nx+1
      Nyp1=Ny+1
      dp=(ze0-z00)/L0
*
      call lcolor(3,4)
      call minmaxV(NxD,Nxp1,Nyp1,Phi,PhiMax,PhiMin,0)
      write(6,'(a,2f10.3)') '*** Min,Max of Phi(i,j) =',PhiMax,PhiMin
      call minmaxV(NxD,Nxp1,Nyp1,Psi,PsiMax,PsiMin,0)
      write(6,'(a,2f10.3)') '*** Min,Max of Psi(i,j) =',PsiMax,PsiMin
*
* Phi > 0
      z0=0
      ze=ze0
      L=L0*(ze-z0)/(ze0-z00)
      call color(2)
      call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Phi, z0,ze,L)
      call display
* Phi =< 0
      z0=z00
      ze=0
      L=L0*(ze-z0)/(ze0-z00)
      call color(6)
      call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Phi, z0,ze,L)
      call display
* Psi > 0
      z0=0
      ze=ze0
      L=l0*(ze-z0)/(ze0-z00)
      call color(5)
      call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Psi, z0,ze,L)
      call display
* Psi =< 0
      z0=z00
      ze=0
      L=l0*(ze-z0)/(ze0-z00)
      call color(4)
      call TOKOSN(NxD,Nxp1,Nyp1,xpos,ypos,Psi, z0,ze,L)
      call display
      end
************************************************************************
      subroutine traject(z0,VIC,VIS)
*
      complex z0,VIC,VIS
      pi=acos(0.0)*2
      pi2=pi*2
      Ni=10
      NT=Ni*4
*
      x0=real(z0)
      y0=imag(z0)
      uC= real(VIC)
      vC=-imag(VIC)
      uS= real(VIS)
      vS=-imag(VIS)
      dwt=pi2/NT
      iwt=0
      call color(3)
      do it=1,4
         call newpen(it)
         do idt=1,Ni
            wt=iwt*dwt
            sinwt=sin(wt)
            coswt=cos(wt)
            xi=x0+uC*sinwt+uS*coswt
            yi=y0+vC*sinwt+vS*coswt
**            if(iwt.eq.0) call symplt(xi,yi,0.35,-1) 
            call draw(xi,yi)
            iwt=iwt+1
         enddo
      enddo
      wt=iwt*dwt
      sinwt=sin(wt)
      coswt=cos(wt)
      xi=x0+uC*sinwt+uS*coswt
      yi=y0+vC*sinwt+vS*coswt
      call draw(xi,yi)
      call penup
**    hgt=0.35
      hgt=0.25
      call symplt(xi,yi,hgt,-1) 
      call display
      end
************************************************************************
      SUBROUTINE TOKOSN(nd,n,m,xpos,ypos,zval, ZI,ZE,L)
*     
      real xpos(nd,*),ypos(nd,*),zval(nd,*)
      parameter ( lmax=400, lmaxp1=lmax+1 )
      DIMENSION X(3),Y(3),Z(3)                                          
      COMMON /T1/ZEQ(lmaxp1),zmax,zmin,LL
*                                                                       
      IF(L.GT.lmax.OR.L.LT.0) then
         WRITE(6,
     *  '(''   *** (TOKOSN)  9TH ARG. .GT.'',i4,'' .OR. .LT.0 ***'')')
     *        lmax
         STOP
      else
         IF(N.LE.0.OR.M.LE.0) then
            WRITE(16,
     *           '(''   *** (TOKOSN) ARGUMENT N OR M .LE. 0 ***'')')  
            STOP
         else                                                         
            LL=L
            llp1=ll+1
            ZMin=MIN(ZI,ZE)
            zmax=max(zi,ze)
            if(ll.eq.0) then
               dz=0
            else
               dz=(zmax-zmin)/ll
            endif
            DO I=1,Llp1                                              
               ZEQ(I)=ZMin+DZ*(I-1)
            enddo
            DO J=1,M-1         
               x(2)=xpos(1,j)
               y(2)=ypos(1,j)
               z(2)=zval(1,j)
               x(3)=xpos(1,j+1)
               y(3)=ypos(1,j+1)
               z(3)=zval(1,j+1)
               DO I=1,N-1 
**                write(6,*) ' i,j =',i,j
                  X(1)=X(3)                                     
                  Y(1)=Y(3)                                       
                  Z(1)=Z(3)                                             
                  X(3)=X(2)                                             
                  Y(3)=Y(2)                                            
                  Z(3)=Z(2)                                         
                  X(2)=xpos(i+1,j)
                  Y(2)=ypos(i+1,j)
                  Z(2)=zval(i+1,j)
                  CALL TPLOT(X,Y,Z)
**                call display
**                pause
                  X(3)=xpos(i+1,j+1)
                  Y(3)=ypos(i+1,j+1)
                  Z(3)=zval(i+1,j+1)
                  CALL TPLOT(X,Y,Z)                                  
**                call display
**                pause
               enddo
            enddo
         endif                                                      
      endif
      call penup
      RETURN                                                
      END                                                               
************************************************************************
      SUBROUTINE TPLOT(X,Y,T)                                           
*                                                                       
*     ***SUBPROGRAM-1 OF TOKOSN ***                                     
*                                                                       
      parameter ( lmax=400, lmaxp1=lmax+1 )
      COMMON /T1/ Z(lmaxp1),zmax,zmin,L
      DIMENSION X(3),Y(3),T(3)
      dimension XP(lmaxp1),YP(lmaxp1),XQ(lmaxp1),YQ(lmaxp1)          
      DATA XE,YE/1.0E30,1.0E30/                                         
*             
      ZMA=MAX(T(1),T(2),T(3))                                         
      ZMI=MIN(T(1),T(2),T(3))                                         
      IF(ZMA.ge.Zmin.and.ZMI.le.Zmax) then
         if(zma.eq.t(1)) then
            nma=1
            if(zmi.eq.t(3)) then
               nmi=3
               nmd=2
            else
               nmi=2
               nmd=3
            endif
         elseif(zma.eq.t(2)) then
            nma=2
            if(zmi.eq.t(1)) then
               nmi=1
               nmd=3
            else
               nmi=3
               nmd=1
            endif
         else
            nma=3
            if(zmi.eq.t(1)) then
               nmi=1
               nmd=2
            else
               nmi=2
               nmd=1
            endif
         endif
         ZMD=T(NMD)                                                    
         CALL TNODE(ZMA,ZMI,ZMD,IMA,IMI,IMAM,IMIM,IND)                 
**       write(6,*) 'ind,ima-imi =',ind,ima,imim,imam,imi
         IF(IND.ne.0) then
            XA=X(NMA)                                                 
            YA=Y(NMA)                                                  
            ZA=T(NMA)                                                
            XB=X(NMI)                                                
            YB=Y(NMI)                                                
            ZB=T(NMI)                                                
            XC=X(NMD)                                                
            YC=Y(NMD)                                                 
            ZC=T(NMD)                                               
            CX=(XA-XB)/(ZA-ZB)                                       
            CY=(YA-YB)/(ZA-ZB)                                       
            DO I=IMI,IMA                                          
               XP(I)=XB+CX*(Z(I)-ZB)
               YP(I)=YB+CY*(Z(I)-ZB)
**             write(6,*) 'P,i,x,y =',i,xp(i),yp(i)
            enddo 
            if(imi.le.imam.and.zc.ne.zb) then
               CX=(XC-XB)/(ZC-ZB)
               CY=(YC-YB)/(ZC-ZB)
               DO I=IMI,IMAM                                       
                  XQ(I)=XB+CX*(Z(I)-ZB)
                  YQ(I)=YB+CY*(Z(I)-ZB)
**                write(6,*) 'Q,i,x,y =',i,xq(i),yq(i)
               enddo
            endif
            if(imim.le.ima.and.za.ne.zc) then
               CX=(XA-XC)/(ZA-ZC)
               CY=(YA-YC)/(ZA-ZC)
               DO I=IMIM,IMA                                      
                  XQ(I)=XC+CX*(Z(I)-ZC)
                  YQ(I)=YC+CY*(Z(I)-ZC)
**                write(6,*) 'Q,i,x,y =',i,xq(i),yq(i)
               enddo             
            endif                                
            DO I=IMI,IMA
               if(mod(i-1,5).eq.0) then
                  call newpen(3)
               else
                  call newpen(1)
               endif
               IF(I.EQ.IMI.AND.XP(I).EQ.XE.AND.YP(I).EQ.YE) then
                  CALL DRAW(XQ(I),YQ(I))                                
                  XE=XQ(I)                                      
                  YE=YQ(I)                                         
               else
                  IF(I.ne.IMI.or.XQ(I).ne.XE.or.YQ(I).ne.YE) then
                     CALL PENUP                        
                     IF(XP(I).GT.XQ(I)) then
                        CALL DRAW(XQ(I),YQ(I))                     
                        CALL DRAW(XP(I),YP(I))                  
                        XE=XP(I)                                        
                        YE=YP(I)                
                     else
                        CALL DRAW(XP(I),YP(I))                    
                        CALL DRAW(XQ(I),YQ(I))                     
                        XE=XQ(I)                                 
                        YE=YQ(I)
                     endif                                              
                  else                                                 
                     CALL DRAW(XP(I),YP(I))                            
                     XE=XP(I)                                        
                     YE=YP(I)                
                  endif                                   
               endif                                
            enddo                    
         endif
      endif
      END                                                               
************************************************************************
      SUBROUTINE TNODE(ZMA,ZMI,ZMD,IMA,IMI,IMAM,IMIM,IND)               
*                                                                       
*     ***SUBPROGRAM-2 OF TOKOSN ***                                     
*                          
      parameter ( lmaxp1=400+1 )
      COMMON /T1/Z(lmaxp1),zmax,zmin,L
*
**    write(6,'(a,3f10.3)') 'ZMA,ZMI,ZMD =',ZMA,ZMI,ZMD
**    write(6,'(a/3(5(i3,f10.3)/))') 'i,z(i) =',(i,z(i),i=1,L+1)
      if(zmi.gt.zmax) then
         ind=0
      else
         Lp1=L+1
         do i=1,Lp1
**          write(6,'(a,i5,2f10.3)') 'NODE-1 i,z(i) =',i,z(i),zmi
         if(z(i).ge.zmi) goto 10
         enddo
 10      continue
         imi=i
**       write(6,'(a,i5)') '***** imi =',imi
         if(zma.lt.zmin) then
            ind=0
         else
            do i=Lp1,imi,-1
**             write(6,'(a,i5,2f10.3)') 'NODE-2 i,z(i) =',i,z(i),zma
            if(z(i).lt.zma) goto 20
            enddo
 20         continue
            ima=i
**          write(6,'(a,i5)') '***** ima =',ima
            do i=imi,ima
**             write(6,'(a,i5,2f10.3)') 'NODE-3 i,z(i) =',i,z(i),zmd
            if(z(i).ge.zmd) goto 30
            enddo
 30         continue
            imim=i
**          write(6,'(a,i5)') '***** imim =',imim
            do i=imim,imi,-1
**             write(6,'(a,i5,2f10.3)') 'NODE-4 i,z(i) =',i,z(i),zmd
            if(z(i).lt.zmd) goto 40
            enddo
 40         continue
            imam=i
**          write(6,'(a,i5)') '***** imam =',imam
            ind=1
         endif
      endif
      END