function f1(t,x,a,b,c,d)
implicit none
real*8::v1,vb,x0,x,t
complex*16::f1,a,b,c,d

real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi

pi=3.141592653d0
vb=10.0d0
v1=x**2+vb*(cos(pi*x/5.2d0))**2
f1=(0.0d0,1.0d0)*((a+b)/(2.0d0*xh**2)-(1.0d0/xh**2+v1+g1/pi*c*conjg(c)+g/pi*d*conjg(d))*c)
end function


function f2(t,x,a,b,c,d)
implicit none
real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi
real*8::v2,x,t
complex*16::f2,a,b,c,d
pi=3.141592653d0
v2=2.0d0*x**2-3.0d0
f2=(0.0d0,1.0d0)*((a+b)/(2.0d0*xh**2)-(1.0d0/xh**2+v2+g2/pi*c*conjg(c)+g/pi*d*conjg(d))*c)
end function

    
program main
implicit none
integer,parameter::m=500
complex*16::y1(-m-1:m+1),y2(-m-1:m+1)
complex*16::K1(-m-1:m+1),L1(-m-1:m+1),M1(-m-1:m+1),N1(-m-1:m+1)
complex*16::K2(-m-1:m+1),L2(-m-1:m+1),M2(-m-1:m+1),N2(-m-1:m+1)
complex*16::a,b,c,d
real*8::x,t,th
complex*16,external::f1,f2
real*8::limitt
integer::l,j
real*8::p(-m-3:m+3),im(-m-3:m+3),re(-m-3:m+3)
real*8::ilm,rle,irm,rre,NL,NR,N,cta1

real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi

xh=0.01d0
th=0.0001d0
g1=0.1d0
g2=0.0d0
g=0.0d0
limitt=100.0d0
pi=3.141592653d0
l=0
open(unit=1,file='1.txt')
open(unit=2,file='2.txt')

do j=-m-1,m+1,1
  y1(j)=(exp(-(j*xh+2.0d0)**2))*(0.55d0*pi)**(-0.25)+(exp(-(j*xh-2.0d0)**2))*(220.0d0*pi)**(-0.25)
  y2(j)=exp(-0.5d0*(j*xh)**2)*pi**(-0.25)
enddo


do t=0.0d0,limitt,th
   do j=-m,m,1
   x=j*xh
   K1(j)=f1(t,j*xh,y1(j+1),y1(j-1),y1(j),y2(j))
   K2(j)=f2(t,j*xh,y2(j+1),y2(j-1),y2(j),y1(j))
   enddo

   do j=-m,m,1
    x=j*xh
    L1(j)=f1(t+0.5d0*th,j*xh,y1(j+1)+0.5d0*th*K1(j+1),y1(j-1)+0.5d0*th*K1(j-1),y1(j)+0.5d0*th*K1(j),y2(j)+0.5d0*th*K2(j))
    L2(j)=f2(t+0.5d0*th,j*xh,y2(j+1)+0.5d0*th*K2(j+1),y2(j-1)+0.5d0*th*K2(j-1),y2(j)+0.5d0*th*K2(j),y1(j)+0.5d0*th*K1(j))
   enddo

   do j=-m,m,1
    x=j*xh
    M1(j)=f1(t+0.5d0*th,j*xh,y1(j+1)+L1(j+1)*0.5d0*th,y1(j-1)+L1(j-1)*0.5d0*th,y1(j)+L1(j)*0.5d0*th,y2(j)+L2(j)*0.5d0*th)
    M2(j)=f2(t+0.5d0*th,j*xh,y2(j+1)+L2(j+1)*0.5d0*th,y2(j-1)+L2(j-1)*0.5d0*th,y2(j)+L2(j)*0.5d0*th,y1(j)+L1(j)*0.5d0*th)
   enddo

   do j=-m,m,1
    x=j*xh
   N1(j)=f1(t+th,j*xh,y1(j+1)+th*M1(j+1),y1(j-1)+th*M1(j-1),y1(j)+th*M1(j),y2(j)+th*M2(j))
   N2(j)=f2(t+th,j*xh,y2(j+1)+th*M2(j+1),y2(j-1)+th*M2(j-1),y2(j)+th*M2(j),y1(j)+th*M1(j))
   enddo

   do j=-m,m,1
   x=j*xh
   y1(j)=y1(j)+th/6.0d0*(K1(j)+2.0d0*L1(j)+2.0d0*M1(j)+N1(j))
   y2(j)=y2(j)+th/6.0d0*(K2(j)+2.0d0*L2(j)+2.0d0*M2(j)+N2(j))
    p(j)=y1(j)*conjg(y1(j))
   im(j)=p(j)*aimag(y1(j))
   re(j)=p(j)*real(y1(j))
   enddo
   
     do j=-m,-1,1  
     ilm=ilm+0.5d0*xh*(im(j)+im(j+1))
     rle=rle+0.5d0*xh*(re(j)+re(j+1))
     NL=NL+0.5d0*xh*(p(j)+p(j+1))
     enddo
   
     do j=0,m,1
     irm=irm+0.5d0*xh*(im(j)+im(j+1))
     rre=rre+0.5d0*xh*(re(j)+re(j+1))
     NR=NR+0.5d0*xh*(p(j)+p(j+1))
    end do
   
    N=NL-NR
   cta1=(datan2(ilm,rle)-datan2(irm,rre))/pi  
   if (cta1>=0.9d0)then
     cta1=cta1-2.0d0
   endif
   if (cta1<=-0.9d0)then
     cta1=cta1+2.0d0
   endif

  if (mod(l,100)==0) then
   write(1,*)t,cta1
   write(2,*)t,N
   endif

  l=l+1
 
  ilm=0.0d0
  rle=0.0d0
  irm=0.0d0
  rre=0.0d0
  NL=0.0d0
  NR=0.0d0

  end do
 stop
 end