c  a-two-dimensional electromagnetic plasma simulation code
c  *************  written by esam ahmad tawfiq  ***********   
c ***** program to solve the plasma parameters in the preasent of electric 
c             ***** and magnetic fields ****** 
c    ***************  input variables ***************


c    p ........ number of particles
c    dt ........ time step
c    nt ........ number of time step to run (  ending  time=nt*dt )
c    nth ....... number of time steps betwween history plots
c    mmax ...... maximum number of different fourier mods to plot
c    nspm ...... maximum number of species
c    iw ........ mover algorithm selector
c    psi ....... multiplier in poisson equation
c    wp ........ plasma frequency
c    wc ........ cyclotron frequency
c    qm ........ q/m  charge:mass ratio
c    lx ........ length of the system in the x-direction
c    ly ........ length of the system in the y-direction
c    ngx ....... number of grid points in the x-direction
c    ngy ....... number of grid points in the y-direction
c    vx0 ....... drift velocity in x-direction
c    vy0 ....... drift velocity in y-direction
c    x1 ........ for loading sinusoidal perturbation in the density 
c                 in the x-direction
c    y1 ......... for loading sinusoidal perturbation in the density
c                 in the y-direction
c    irho ...... plotting interval for  rho (charge density). =0 for no plot
c    irhos ..... plotting interval for smoothed density
c    iphi ......  plotting interval for phi (potential )
c    ie ........ plotting interval for e (electric field )
c    ixvx ...... plotting interval for x vs.vx phase space
c    ifvx ...... plotting  intervar for f(vx) distribution function
c    ivxvy ..... plottinginterval for  vx vs. vy  phase spase
c    mpolt ..... fourier mode numbers to plot
c    mode,x1,v1,thetax and thetay  are for loading 
c    a sinusoidal perturbation         
	 integer nspm,mmax,nth1,nth2,ih
	  parameter(mngx=65,mngy=64)
          parameter (p=25602,ih=4000)
	  common /param/ nsp,lx,ly,dx,dy,dt,nt
       common /cfield/ ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)   
     $,bz(34,34)
c        particle coordinate and velocities
	 common /xvp/ x(p),y(p),vx(p),vy(p)
	 dimension ek(34),ehk(34)
	 common /cntrl/ it,time,ithlx,ithly,iex,iey,ixvx,iyvy,
     $                ivxvy,ifvx,ifvy           
	parameter(nth=4000,mmax=10,nspm=3)
	parameter(nth1=nth+1,nth2=nth+2,nspm1=nspm+1)
	common ms(nspm),t(nspm),qs(nspm),ts(nspm),vts(nspm),nms(nspm)
C       nth= number of time step between history plots.
c       mmax= maximum number of different fourier modes to plot.
c       nsp= maximum number of species.
      real lx,ly,rho0,ek,ehk,kl,ml,wo,wpe,ld,io,te1,wcs(ns
     $pm),ms,qs,ts,yn,xn,nms,dt,ke1,ese1,wps(nspm),qe,me,lxr,lyr
     $,eyyn(34),eyyo(34),bz,bbz(34),eyy(34),EY1(34),DRE,TE,dcofc,dcof
     $,vthr,zz,ne,kw,kw1,te1o,ke1o,dreo,eseo
	  integer ins(4),nsp,ngy,it,nt,ith,ngx,ie,irho,ixy,ifvx,ihis  
     $    ,lat,npr,z
c     npr= number of laser periodesin the plasma
          character*15 nf1,nf2,nf3
c          record /xycoord/ xy
          open(1,file='em2db.dat')
	  read(1,1000)a
          read(1,*) io,kl,ml,wo,wpe,ld,lat,npr,z,ne
	  read(1,1000)a
	  read(1,*) nsp,lx,ly,ngy,ngx,nt,dt,ie,irho,ixy,ifvx,ihis
	  
c     **************  default input parameter  *****************
	ins(1)=1   
	pi=3.14
	rho0=0.    
	qe=1.6e-19   
	c=3e+8
	me=9.1e-31
	lxr=lx*ld
	lyr=ly*ld
	dtr=dt/wpe                    
c     lx= simul. length
c **** epo = dimensionless electric field *****   
c *******    dt======== real time step **********        
	   dx=lx/ngx
	   dy=ly/ngy
	   dxi=1/dx
	   dyi=1/dy
	   yn=ngy
	   xn=ngx
c   absorption coeff.
C*****************************************************************************
	   do 10 is=1,nsp
10       call init(ins(is),ins(is+1),ms(is),qs(is),wcs(is),
     $nms(is),vts(is),wps(is))
           te1=1.
           zz=26*sqrt(.001*te1/(1+.001*(26/z)**2*te1))
           ne=(wpe*2*3.14/5.64e4)**2
           kw=5.64e-11*ne**2*7/(wo*wpe)**2*1/sqrt(1-1/wo**2)
           epo=-(qe/(me*wpe*c)*ms(1)/qs(1)*wps(1))
          
          eo=epo*27.5*sqrt(iO)
          read(1,1000)a
          read(1,35) nf1,nf2,nf3
          write(*,35) nf1,nf2,nf3
          open(3,file=nf1)
	  open(4,file=nf2)
          open(5,file=nf3)
          open(6,file='rrr')
*****************************************************************************             
c            call plotxy(0,ins(1),ins(2),ngx,ngy,ixy)
******************************************************************************            
	    call field(0)     
****************************************************************************              
	   do is=1,nsp
          call setv(ins(is),ins(is+1)-1,ms(is),qs(is),wcs(is))
	  end do

*****************************************************************************
*****************   begain time step loop      *****************************
********* advance velocities from it-0.5 to it+0.5 **************
	 
100      continue          
	   
c###################################################################        
          zz=26*sqrt(.001*te1/(1+.001*(26/z)**2*te1))
          dcof=kw1*zz/(te1)**3/2*eo
          write(6,*) te1,dcof
          eo=eo-dcof
         call  fvxy(ins(1),ins(NSP+1)-1,vts(1)*dt/dx,it,dx/dt,ifvx)

          do 12 is=1,nsp
12       call accel (ins(is),ins(is+1)-1,ms(is),qs(is),
     $ke1,te1,DRE,dcof)
        
c#######################################################################
*****  advance positions from it to it+1 ***********
        
	 do 13 is=1,nsp
13        call move(ins(is),ins(is+1)-1,qs(is),xn,yn)
          call plotxy(it,ins(1),ins(nsp+1)-1,ngx,ngy,dx/dt,ixy)
c           if(it.eq.200)stop
          call plotorb(ngx,ngy)
          if(it.ge.nt) go to 30
	  it=it+1
	  time=it*dt
	  print*,'it',it
c***** get fields at time step it *******
          call fieldf(it)
	  call field(it)
          ESE=0.
	  DO I=1,NGY+1
	     DO J=1,NGX+1
	       ESE=ESE+EX(I,J)*EX(I,J)+EY(I,J)*EY(I,J)
     
	     END DO
	 END DO         
         ESE=0.5*MS(1)*WPS(1)*WPS(1)*LX*(INS(2)-INS(1))
     $    /(NGX*NGx*NGX)*ESE
                    
        TE=ESE+KE1+DRE
        if(it.eq.1) then
          te1o=te1
          ke1o=ke1
          eseo=ese
          dreo=dre
        end if
        print*,te1/te1o,ke1/ke1o,ese/eseo,dre/ke1o
c        te1=sqrt(2*te1/(ms(1)*ins(2)))*1
        write(4,40)it,ke1/ke1o,ese/eseo,te1/te1o,DRE/ke1o,TE/ke1o



c#############  adding laser electric field ###############          
C          if(it.ge.lat) go to 16
c	  call settextposition(2,1,xy) 
	 do j=1,ngx+1
c ***** eyy rpresents the electric field of the laser pulse ********
                eyyn(J)=eo*sin(wo*it*dt-(j-1)*2*npr*pi*dx/lx)
             do  i=1,ngy+1
	       ey(i,j)=ey(i,j)+eyyn(J)
	       bz(i,j)=eyyn(j)*qs(1)/ms(1)
	     end do
	  end do
	  do j=1,ngx+1
	     bbz(j)=bz(16,j)
	  end do 
16        do j=1,ngx
	   ek(j)=ey(16,j)*1/dx
	   
	   ehk(j)=0.
	end do
	do j=1,ngx
	eyy(j)=ey(16,j)-eyyn(j)
	EY1(J)=EY(16,J)
	end do
         
	call lasf(it,eyy,EY1,bbz,ngx) 
	 call plotex(it,ie)          
	 call plotey(it,ie)
	 go to 100
***** end of run ********
30       continue
 31      format(i4,8f8.4)
 32      format(12f8.4)
 33      format(10f8.4)
 35      format(7a15)
 40      format(i5,3x,5e10.3)
1000     format(80a1)
	 end  
							 
c************************************************************            
        subroutine init(il1,il2,m,q,wc,nm,vt,wp) 
c************************************************************        
	common /param/nsp,lx,ly,dx,dy,dt,nt
	common /cfield/ ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)   
	common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
******  replaces uniform distribution plasma *******              
	 real lx,ly,m,q,p,wp,vx0,vy0,x1,y1,vt,rho0,nm,yn,xn
     $    ,dny,dnx,dyy,dxx,gs,dx,dy,qc,qdxy,rand,vv,vxy,wc,to
	 real ngx,ngy,k,j1,i1,iw,ny,nx,a,seed
	 read(1,1000)a
	 read(1,*) p,wp,wc,qm,vx0,vy0,x1,y1,vt,iw
	 t=tan(-wc*dt/2)
	 seed=5773
	 il2=il1+p
	 q=ly*wp*wp/(p*qm)
	 m=q/qm
	nm=p*m
	ny=sqrt(p)*ly/lx
	nx=sqrt(p)*lx/ly
	dny=ly/ny
	dnx=lx/nx
	dxi=1/dx
	dyi=1/dy
	yn=ngy
	xn=ngx
	 do i=il1,il2-1
	   vx(i)=0.
	   vy(i)=0.
	end do
	
	if(vt.eq.0.) go to 45
	vv=0.
	vxy=0.
	do i=il1,il2-1
	 z=2*3.142*rand(seed)
	 v=sqrt(-2*ALOG(RAND(SEED)))
	 v=vt*v
	 vy(i)=v*cos(z)
	 vx(i)=v*sin(z)
	 vv=vv+v*v
	 vxy=vxy+vx(i)*vx(i)+vy(i)*vy(i)
	 end do
	 
 45    if(vx0.eq.0.) go to 47
	do 2 i=il1,il2-1
	vx(i)=vx0+vx(i)
2       continue
47      if(vy0.eq.0.) go to 50
	do 3 i=il1,il2-1
	vy(i)=vy0+vy(i)
3       continue 
50       k=il1-1
	do 4 i=1,ny
	yk=(i-0.5)*dny
	do 5 j=1,nx
	k=k+1
	x(k)=(j-0.5)*dnx
	y(k)=yk
 5       continue
 4       continue
	 k=il1-1
	if(y1.le.0.)go to 30
	do 6 i=1,ny
	dyy=y1*sin(3.14*2*y(k+1)/ly)
	do 7 j=1,nx
	 k=k+1
	 y(k)=y(k)+dyy
7       continue
6       continue
30      if(x1.le.0.)go to 60
	do 8 j=1,nx
	dxx=x1*sin(3.14*2*x(j)/lx)
	do 9 i=1,ny
	x((i-1)*nx+j)=x((i-1)*nx+j)+dxx
9       continue
8       continue
	il1=1
60      gs=p*q/(lx*ly)
	qdxy=q/(dx*dy)
	if(il1.ne.1)go to 122
	do 10 i=1,ngy+1
	do 10 j=1,ngx+1                  
10      g1(i,j)=0.
122     rho0=rho0-gs
************************************************************************
************************************************************************
	do 11 i=1,ngy
	do 11 j=1,ngx
11         g1(i,j)=g1(i,j)-gs
	go to (100,200),iw
*****  collect charge densities to the grid points  ********

******************** N G P ************************
100    continue
       do 12 i=il1,il2-1
       x(i)=x(i)*dxi
       y(i)=y(i)*dyi
       if(x(i).lt.0.) x(i)=x(i)+xn
       if(x(i).gt.xn) x(i)=x(i)-xn
       if(y(i).lt.yn) y(i)=y(i)+yn
       if(y(i).gt.yn) y(i)=y(i)-yn
       j1=x(i)+0.5
       i1=y(i)+0.5
       g1(i1+1,j1+1)=g1(i1+1,j1+1)+qdxy
12     continue
       return
******************** LINEAR *************************
200    continue
       qc=qdxy
       do 13 i=il1,il2-1
       x(i)=x(i)*dxi
       y(i)=y(i)*dyi
       if(x(i).lt.0.) x(i)=x(i)+xn
       if(x(i).gt.xn) x(i)=x(i)-xn
       if(y(i).lt.0.) y(i)=y(i)+yn
       if(y(i).gt.yn) y(i)=y(i)-yn
       j1=x(i)
       i1=y(i)
       g1(i1+2,j1+2)=qc*((y(i)-i1)*(x(i)-j1))+g1(i1+2,j1+2)
       g1(i1+1,j1+1)=qc*((i1+1-y(i))*(j1+1-x(i)))+g1(i1+1,j1+1)
       g1(i1+2,j1+1)=qc*((i1+1-y(i))*(x(i)-j1))+g1(i1+2,j1+1)
 13    g1(i1+1,j1+2)=qc*((y(i)-i1)*(j1+1-x(i)))+g1(i1+1,j1+2)
       return
1000   format(80a1)
	end
c######################################################################
      function rand(k)
      integer k,m,const1
      real rand,const2
      parameter(const1=2147483647,const2=.4656613e-9)
      save
      data m/0/
      if(m.eq.0) m=k
      m=m*65539
      if(m.lt.0) m=(m+1)+const1
      rand=m*const2
      end
c#####################################################################
**************************************************************************
       subroutine field (ith)
**************************************************************************      
	common/param/ nsp,lx,ly,dx,dy,dt,nt
	common/cfield/ ngx,ngy ,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)   
	common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
	common/cntrl/it,time,ithlx,ithly,iex,iey,ixvx,iyvy,
     $  ivxvy,ifvx,ifvy          
	integer ith,l 
	real lx,ly,pi,erun1,erun2
	 gx=ngx
	gy=ngy
	pi=2.*3.14
	ityp=2
	h=ityp-1
	mn=ngy
	nn=ngx
	m=mn-2
	n=nn-2
	mh=m/2
	nh=n/2
	m1=m+2
	n1=n+2
	m0=m
	np=n+1
	mp=m+1
	h1=1.
	h2=1.
	
c        do 5 i=1,ngy
c  5     g1(i,1)=g1(i,1)+g1(i,ngx+1)
c        Do 6 j=1,ngx
c  6     g1(1,j)=g1(1,j)+g1(ngy+1,j)
	      
	do 9 j=1,ngx
	erun2=0.
	do 10 k=1,ngy
	erun1=0.
	do 12 i=1,k
	l=i+1
	if(l.gt.ngy) l=1
  12    erun1=erun1+.5*dy*(g1(i,j)+g1(l,j))
  10    erun2=erun2+erun1
  9     ey(1,j)=-erun2/ngy
  
	do 19 i=1,ngy
	erun2=0.
	do 20 k=1,ngx
	erun1=0.
	do 22 j=1,k
	l=j+1
	if(l.gt.ngx) l=1
  22    erun1=erun1+.5*dx*(g1(i,j)+g1(i,l))
  20    erun2=erun2+erun1
  19    ex(i,1)=-erun2/ngx
	do 30 i=1,ngy 
	do 30 j=1,ngx-1
  30    ex(i,j+1)=ex(i,j)+0.5*dx*(g1(i,j+1)+g1(i,j))
	do 32 j=1,ngx 
	do 32 i=1,ngy-1
  32    ey(i+1,j)=ey(i,j)+0.5*dy*(g1(i+1,j)+g1(i,j))
	do 33 j=1,ngx
  33    ey(ngy+1,j)=ey(1,j)
	do 34 i=1,ngy
  34    ex(i,ngx+1)=ex(i,1)
**********************************************
****** centered difference across two cells ******
	do 27 i=1,ngy
	do 27 j=1,ngx
 27       g1(i,j)=rho0
	do  131 i=1,ngy+1
  131   g1(i,ngx+1)=0.
	do  132 j=1,ngx+1
  132   g1(ngy+1,j)=0.
	ael=1.
	return
	end


*****************************************************************************
**************************************************************************
       subroutine fieldf (ith)
**************************************************************************      
       include 'fgraph.fd'
	common/param/ nsp,lx,ly,dx,dy,dt,nt
	common/cfield/ ngx,ngy ,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)   
	common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
	common/cntrl/it,time,ithlx,ithly,iex,iey,ixvx,iyvy,
     $  ivxvy,ifvx,ifvy          
	integer*2 dummy
       record /wxycoord/ wxy  
	integer ith,l,ng2,hdx,hdy,ik,jk 
	real lx,ly,pi,erun1,erun2,kdx2,sm(34),ksq(34),eset,phif(34,34)
     $,rhof(34,34),sc(34),rhok(34),lxlyi
	data ng2,pi/0,3.128/
	if(ng2.ne.0) go to 2
	a1=0.
	a2=0.
	ng2=ngx/2
	do k=1,ng2
	   kdx2=(pi/ngx)*k
	   sm(k)=exp(a1*sin(kdx2)**2-a2*tan(kdx2)**4)
	   ksq(k)=1.0/((2.0*sin(kdx2)/dx)**2)*sm(k)**2
	   ksq(k)=1/ksq(k)
       end do
2       continue
	do 5 i=1,ngy
  5     g1(i,1)=g1(i,1)+g1(i,ngx+1)
	Do 6 j=1,ngx
  6     g1(1,j)=g1(1,j)+g1(ngy+1,j)
	hdx=.5*dx
	hdy=.5*dy
	lxlyi=1/(lx*ly)

	do i=1,ngy
	   do j=1,ngx
	      rhok(j)=g1(i,j)*hdx
	      sc(j)=0.
	   end do
	      call cpft(rhok,sc,ngx,1,1)
	      call rpft2(rhok,sc,ngx,1)
	      rhok(1)=0.
	   do j=1,ngx
	      rhof(i,j)=rhok(j)
	   end do
	end do
	do j=1,ngx
	   do i=1,ngy
	      rhok(i)=rhof(i,j)*hdy
	      sc(i)=0.
	   end do
	      call cpft(rhok,sc,ngx,1,1)
	      call rpft2(rhok,sc,ngx,1)
	      rhok(1)=0.
	   do i=1,ngy
	      rhof(i,j)=rhok(i)
	    end do
	 end do
	    
	phif(1,1)=0.
	ese=0.
	do i=2,ngy/2
	   ik=ngy+2-i
	   do j=2,ngx/2
	      jk=ngx+2-j
	      phif(i,j)=rhof(i,j)/(ksq(j-1)+ksq(i-1))
	      phif(ik,jk)=rhof(ik,jk)/(ksq(j-1)+ksq(i-1))
	      phif(i,jk)=rhof(i,jk)/(ksq(j-1)+ksq(i-1))
	      phif(ik,j)=rhof(ik,j)/(ksq(j-1)+ksq(i-1))
	      ese=ese+rhof(i,j)*phif(i,j)+rhof(ik,jk)*phif(ik,jk)
	   end do
	   phif(i,ngx/2)=rhof(i,ngx/2+1)/ksq(ngx/2)   
	   ese=ese+phif(i,ngx/2+1)*rhof(i,ngx/2+1)
	end do   
	eset=ese/(Lx*ly)
	do i=1,ngy
	   do j=1,ngx
	      phif(i,j)=phif(i,j)*rhof(i,j)*lxlyi
	   end do

	end do
        goto 30
c        call plotm(phif,ngx,ngy)
	call setviewport(300,250,600,350)
       dummy=setwindow(.true.,0,100,500,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       dummy=rectangle_w($gborder,0,100,500,-100)  
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=setcolor(9)
       ekm=abs(phif(2,2))
       do j=1,ngx    
	ekm=amax1(ekm,phif(1,j),phif(2,j),phif(3,j))
       end do
       call moveto_w(500.0d0/ngx,phif(1,1)*100*.8d0/ekm,wxy) 
       do j=1,ngx 
	   dummy=lineto_w(j*500.0d0/ngx,phif(1,j)*100*.8d0/ekm)    
       end do
       call moveto_w(500.0d0/ngx,phif(2,1)*100*.8d0/ekm,wxy) 
       do  j=1,ngx 
	   dummy=lineto_w(j*500.0d0/ngx,phif(2,j)*100*.8d0/ekm)    
       end do
       call moveto_w(500.0d0/ngx,phif(3,1)*100*.8d0/ekm,wxy) 
       do  j=1,ngx 
	   dummy=lineto_w(j*500.0d0/ngx,phif(3,j)*100*.8d0/ekm)    
       end do
  30   do j=1,ngx
	  rhok(j)=phif(2,j)
	end do
        write(3,20) ith,(rhok(j),j=1,8)
	

 10     format(8e8.4)
 20     format(i4,8e10.3)  
	return
	end

***************************************************************************
	 
	 
	 subroutine setv(il,iu,m,q,t)
***************************************************************************
******* converts particle velocities at t=0 to computer normalization

******** of t=-dt/2 *********
	 common/param/nsp,lx,ly,dx,dy,dt,nt
	 common /cfield/ ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)    
	 common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
	 integer iu,il
	 real q,m,t,ke,lx,ly,dx,dy,dt,DRE
	 dtdx=dt/dx
	 if(t.eq.0.) go to 20
	 c=1/sqrt(1+t*t)
	 s=c*t
	 do 1 i=il,iu-1
	 vxx=vx(i)
	 vx(i)=c*vxx+s*vy(i)
1        vy(i)=-s*vxx+c*vy(i)
20       do 2 i=il,iu-1
	 vx(i)=vx(i)*dtdx
2        vy(i)=vy(i)*dtdx
******************************************************************************
         call accel(il,iu,m,-.5*q,ke,Te,DRE,0.)            
******************************************************************************
	 return
	  end
**************************************************          
         subroutine accel(ilp,iup,m,q,ke,te,DRE,Einb)            
*************************************************
	  common /param/nsp,lx,ly,dx,dy,dt,nt
	  common/cfield/ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)     
	  common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
          real m ,q ,ke,ae,tem,ael,exx,exxx,eyy,eyyy,vxvy,lx
     $,ly,dx,dy,dt,ex,ey,dre,sumx,sumy,te,bzz,bzzz,einb,m2
	  integer j1,i1 
	  ael=1.
*******  advance velocity one time step
	 dxdt=dx/dt
	 ae=(q/m)*dt/dxdt
	 AE=0.5*ae
	 if(ae.eq.ael) go to 60
	  tem=ae/ael
	 do 1 i=1,ngy+1
	 do 1 j=1,ngx+1
	 ex(i,j)=ex(i,j)*tem
	 ey(i,j)=ey(i,j)*tem 
1        continue
	 ael=ae
60       continue
******** select acceleration weighting *******
c         goto 300
        do 3 i=ilp,iup-1
	i1=y(i)
	j1=x(i)
	exx=ex(i1+2,j1+2)*(y(i)-i1)*(x(i)-j1)+ex(i1+2,j1+1)*(y(i)-i1)*
     $    (j1+1-x(i))+ex(i1+1,j1+2)*(i1+1-y(i))*(x(i)-i1)                                                 
	 exxx=exx+ex(i1+1,j1+1)*(i1+1-y(i))*(j1+1-x(i)) 
	 eyy=ey(i1+2,j1+2)*(y(i)-i1)*(x(i)-j1)+ey(i1+2,j1+1)*(y(i)-i1)*
     $    (j1+1-x(i))+ey(i1+1,j1+2)*(i1+1-y(i))*(x(i)-i1)
	 eyyy=eyy+ey(i1+1,j1+1)*(i1+1-y(i))*(j1+1-x(i))
	bzz=bz(i1+2,j1+2)*(y(i)-i1)*(x(i)-j1)+bz(i1+2,j1+1)*(y(i)-i1)*
     $    (j1+1-x(i))+bz(i1+1,j1+2)*(i1+1-y(i))*(x(i)-i1)                                                 
	bzzz=bzz+bz(i1+1,j1+1)*(i1+1-y(i))*(j1+1-x(i)) 
	c=cos(bzzz*dt)
	s=sin(bzzz*dt)
	vy(i)=vy(i)+eyyy
	 vx(i)=vx(i)+exxx
	 vyy=-vx(i)*s+vy(i)*c
	 vxx=vx(i)*c+vy(i)*s
	 vy(i)=vyy+eyyy
3        vx(i)=vxx+exxx    
	 vxvy=0.
c         am=0.5*wp2m*(lx/ngx)**2
         Einb=Einb*.5/(iup-1)*1/(dxdt*dxdt)
         m2=m/2
         do i=ilp,iup-1
            if(vx(i).lt.0.) then
               vx(i)=-sqrt(vx(i)*vx(i)+einb/m2)
            else
               vx(i)=sqrt(vx(i)*vx(i)+einb/m2)
            end if
            if(vy(i).lt.0.) then
               vy(i)=-sqrt(vy(i)*vy(i)+einb/m2)
            else
               vy(i)=sqrt(vy(i)*vy(i)+einb/m2)
            end if
         end do
          do 4 i=ilp,iup-1
    4    vxvy=vxvy+(vx(i)*vx(i)+vy(i)*vy(i))
	 ke=.5*m*vxvy*dxdt*dxdt 
	 
c****** te= electron plasma teparature**************         
	 te=vxvy*dxdt*dxdt*m*.5/(iup-ilp)
c  ***** to find particle drift energy ******         
       sumx=0.
       sumy=0.
       do i=ilp,iup-1
	if(vx(i).gt.0) then 
	 sumx=sumx+vx(i)*vx(i) 
	else 
         sumx=sumx-vx(i)*vx(i)
	end if
	if(vy(i).gt.0) then 
	 sumy=sumy+vy(i)*vy(i)
	else
	 sumy=sumy-vy(i)*vy(i)
	end if
       end do
       dre=0.5*m*abs((sumx+sumy))*dxdt*dxdt
       return
       end
**************************************************************************                    
	 subroutine move (ilp,iup,q,xn,yn)
************************************************************************
	  common/param/nsp,lx,ly,dx,dy,dt,nt
	  common/cfield/ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)     
	  common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
******  advances position one time step and accumulate charge density*****
	  integer j1,i1,ilp,iup
	  real qc,qdxy
	  qdxy=q/(dx*dy)
	  go to (100,200),iw
*************** N G P ****************
100      continue
	 do 1 i=ilp,iup 
	 x(i)=x(i)+vx(i)
	 y(i)=y(i)+vy(i) 
	 if(x(i).lt.0.) x(i)=x(i)+xn
	 if(x(i).gt.xn) x(i)=x(i)-xn
	 if(y(i).lt.0.) y(i)=y(i)+yn
	 if(y(i).gt.yn) y(i)=y(i)-yn
	 j1=x(i)+0.5
	 i1=y(i)+0.5
 1       g1(i1+1,j1+1)=g1(i1+1,j1+1)+qdxy
 
	 return
********** LINEAR ***********
200      continue    
	 qc=qdxy
       do i=ilp,iup
	 x(i)=x(i)+vx(i)
	 y(i)=y(i)+vy(i)
	 if(x(i).lt.0.) x(i)=x(i)+xn
	 if(x(i).gt.xn) x(i)=x(i)-xn 
	 if(y(i).lt.0.) y(i)=y(i)+yn
	 if(y(i).gt.yn) y(i)=y(i)-yn
	 j1=x(i)
	 i1=y(i)
	 g1(i1+2,j1+2)=qc*((y(i)-i1)*(x(i)-j1))+g1(i1+2,j1+2)
	 g1(i1+1,j1+1)=qc*((i1+1-y(i))*(j1+1-x(i)))+g1(i1+1,j1+1)
	 g1(i1+2,j1+1)=qc*((i1+1-y(i))*(x(i)-j1))+g1(i1+2,j1+1)
	 g1(i1+1,j1+2)=qc*((y(i)-i1)*(j1+1-x(i)))+g1(i1+1,j1+2)

       end do   
	
	return
	 end

C#####################################################
*****************************************************************************
       subroutine plot(it,irho)         
*****************************************************************************       
       include 'fgraph.fd'
	common /cfield/ ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)  
       integer*2 dummy
       integer it,irho
       real gmax
c       record /xycoord/ xy
	record /wxycoord/ wxy
	if((it/irho)*irho.ne.it) return
	call setviewport(0,0,300,200)
       dummy=setwindow(.true.,0,100,600,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=lineto_w(500,0)
       dummy=lineto_w(600,50)  
       dummy=lineto_w(100,50)  
       dummy=lineto_w(0,0)  
c       call _fq_settextposition(20,10,xy)
c       call moveto_w(0,g1(1,1)*20-25,wxy)
       gmax=200. 
       dummy=setcolor(9)
       do 10 j=1,ngx  
       call moveto_w((j-1)*500.0d0/(ngx-1),dble(g1(1,j))*gmax,wxy)
       do  11 i=1,ngy
       
   11 dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100.0d0/(ngy-1),(i-1)
     $*50/(ngy-1)+dble(g1(i,j))*gmax)
10      continue  
	do 12 i=1,ngy
      call moveto_w((i-1)*100.0d0/(ngy-1),(i-1)*50.0d0/(ngy-1)+
     $dble(g1(i,1))*gmax,wxy)   
       do  13 j=1,ngx
      
   13  dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100/(ngy-1),(i-1)
     $*50.0d0/(ngy-1)+g1(i,j)*gmax)
   12     continue  
	
	return
       end
*****************************************************************************       
*****************************************************************************       
*****************************************************************************       

       subroutine plotm(ex,ngx,ngy)         
******************************************************************************
       include 'fgraph.fd'
       integer*2 dummy
c       record /xycoord/ xy
	record /wxycoord/ wxy
	integer ie 
	real emax,ex(34,34)
	emax=ex(1,1)
c       if(it.le.1) dummy=setvideomode($vres16color)
c        if((it/ie)*ie.ne.it) return
c       do i=1,ngy
c          do j=1,ngx
c            emax=amax1(emax,abs(ex(i,j)))
c          end do
c       end do
       emax=1e-4
       
       if(emax.le.0) emax=.0001
       call setviewport(100,0,400,200)
       dummy=setwindow(.true.,0,100,600,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=lineto_w(500,0)
       dummy=lineto_w(600,50)  
       dummy=lineto_w(100,50)  
       dummy=lineto_w(0,0)  
c      call _fq_settextposition(20,10,xy)
c      call moveto_w(0,ex(1,1)*20-25,wxy)
       dummy=setcolor(9)
       do 10 j=1,ngx  
       call moveto_w((j-1)*500.0d0/ngx,dble(ex(1,j))*50.0d0/emax,wxy)
       do  11 i=1,ngy,2
11     dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100.0d0/(ngy-1),(i-1)
     $*50.0d0/(ngy-1)+ex(i,j)*50/emax)
10     continue  
       do 12 i=1,ngy,2
      call moveto_w((i-1)*100.0d0/(ngy-1),(i-1)*50.0d0/(ngy-1)+ex(i,1)*
     $50/emax,wxy)   
       do  13 j=1,ngx
13     dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100/(ngy-1),(i-1)
     $*50.0d0/(ngy-1)+ex(i,j)*50/emax)
12     continue  
       return
       end

******************************************************************************
       
       
       subroutine plotex(it,ie)         
******************************************************************************
       include 'fgraph.fd'
	common /cfield/ ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)  
       integer*2 dummy
c       record /xycoord/ xy
	record /wxycoord/ wxy
	integer ie 
	real emax
	emax=ex(1,1)
c       if(it.le.1) dummy=setvideomode($vres16color)
	if((it/ie)*ie.ne.it) return
       do i=1,ngy
	  do j=1,ngx
	    emax=amax1(emax,abs(ex(i,j)))
	  end do
       end do
       emax=1.5*emax
       
       if(emax.le.0) emax=1
       call setviewport(0,0,300,200)
       dummy=setwindow(.true.,0,100,600,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=lineto_w(500,0)
       dummy=lineto_w(600,50)  
       dummy=lineto_w(100,50)  
       dummy=lineto_w(0,0)  
c      call _fq_settextposition(20,10,xy)
c      call moveto_w(0,ex(1,1)*20-25,wxy)
       dummy=setcolor(9)
       do 10 j=1,ngx  
       call moveto_w((j-1)*500.0d0/ngx,dble(ex(1,j))*50.0d0/emax,wxy)
       do  11 i=1,ngy,2
11     dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100.0d0/(ngy-1),(i-1)
     $*50.0d0/(ngy-1)+ex(i,j)*50.0d0/emax)
10     continue  
       do 12 i=1,ngy,2
      call moveto_w((i-1)*100.0d0/(ngy-1),(i-1)*50.0d0/(ngy-1)+ex(i,1)* 
     $50.0d0/emax,wxy)   
       do  13 j=1,ngx
13     dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100.0d0/(ngy-1),(i-1)
     $*50.0d0/(ngy-1)+ex(i,j)*50.0d0/emax)
12     continue  
       return
       end
*****************************************************************************       
******************************************************************************
       subroutine plotey(it,ie)         
******************************************************************************
       include 'fgraph.fd'
	common /cfield/ ngx,ngy,iw,rho0,g1(34,34),ex(34,34),ey(34,34)
     $,bz(34,34)  
       integer*2 dummy
c       record /xycoord/ xy
	record /wxycoord/ wxy
	integer ie
	real emax
	emax=ey(1,1)
c      if(it.le.1) dummy=setvideomode($vres16color)
       if((it/ie)*ie.ne.it) return
       do i=1,ngy
	  do j=1,ngx
	    emax=amax1(emax,abs(ey(i,j)))
	  end do
       end do
       emax=1.5*emax
       
       call setviewport(300,0,600,200)
       dummy=setwindow(.true.,0,100,600,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=lineto_w(500,0)
       dummy=lineto_w(600,50)  
       dummy=lineto_w(100,50)  
       dummy=lineto_w(0,0)  
c      call _fq_settextposition(20,10,xy)
c      call moveto_w(0,eY(1,1)*20-25,wxy)
       dummy=setcolor(9)
       do 10 j=1,ngx  
       call moveto_w((j-1)*500.0d0/ngx,eY(1,j)*50.0d0/emax,wxy)
       do  11 i=1,ngy,2
11     dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100.0d0/(ngy-1),(i-1)
     $*50.0d0/(ngy-1)+ey(i,j)*50.0d0/emax)
10     continue  
       do 12 i=1,ngy,2
      call moveto_w((i-1)*100.0d0/(ngy-1),(i-1)*50.0d0/(ngy-1)+ey(i,1)* 
     $50.0d0/emax,wxy)   
       do  13 j=1,ngx
13    dummy=lineto_w((j-1)*500.0d0/(ngx-1)+(i-1)*100.0d0/(ngy-1),(i-1)
     $*50.0d0/(ngy-1)+ey(i,j)*50.0d0/emax)
12    continue  
       return
       end
c*****************************************************************************       
       subroutine plotxy(it,io,in,ngx,ngy,dxdt,ixy)     
c******************************************************************************
       include 'fgraph.fd'
       common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
       integer*2 dummy
       real dxdt
       integer it,io,in,ngx,ngy,ixy
       record /xycoord/xy
       record /wxycoord/wxy
       if((it/ixy)*ixy.ne.it) return
	if(it.le.1) dummy=setvideomode($vres16color)
       call setviewport(300,100,600,200)
       dummy=setwindow(.true.,0,100,600,-100)
c       call _fq_settextposition(5,5,xy)
       print*,'it=',it
       call clearscreen($gviewport)
       dummy=setcolor(9)
       dummy=rectangle_w($gborder,0,100,600,-100)  
       dummy=setcolor(9)
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=lineto_w(600,0)
c      call _fq_settextposition(20,10,xy)
       do 10 i=io,in 
10     dummy=setpixel_w(x(i)*600.0d0/ngx,VY(i)*dxdt*100.0d0/10.0d0)    
       return 
       end
c****************************************************************************      
	subroutine plotorb(ngx,ngy)     
c******************************************************************************
       include 'fgraph.fd'
       common /xvp/ x(25602),y(25602),vx(25602),vy(25602)
       integer*2 dummy
       integer ngx,ngy
       record /xycoord/ xy
c       if((it/ixy)*ixy.ne.it) return
c        if(it.le.1) dummy=setvideomode($vres16color)
       call setviewport(300,350,600,450)
       dummy=setwindow(.true.,0,100,600,0)
c       call _fq_settextposition(5,5,xy)
       print*,'it=',it
c       call clearscreen($gviewport)
       dummy=setcolor(9)
       dummy=rectangle_w($gborder,0,200,600,0)  
c      call _fq_settextposition(20,10,xy)
       do 10 i=10000,10000 
       if(x(i).lt.0.or.x(i).gt.ngx) goto 10
       if(y(i).lt.0.or.y(i).gt.ngy) goto 10
       dummy=setpixel_w(x(i)*600.0d0/ngx,Y(i)*100.0d0/ngy)    
 10    continue      
       return 
       end

c#########################################################
       subroutine  fvxy(io,in,vt,it,dxdt,ifvx)
c##########################################################
       include 'fgraph.fd' 
       common /xvp/ x(25602),y(25602),vx(25602),vy(25602) 
       integer*2 dummy
       record /wxycoord/ wxy
       record/xycoord/ xy
       integer  fvx(101),fvy(101),io,in,jjx,jjy,it
       real dv,vmax,vmin,vt,xj,yj,dxdt
       if((it/ifvx)*ifvx.ne.it) return
       vmax=10.0/dxdt
       vmin=-10.0/dxdt
       dv=(vmax-vmin)/100
       
       if(dv.eq.0.) return

       do  j=1,100
	fvx(j)=0
	fvy(j)=0
      end  do
      do i=io,in
	 xj=vx(i)/dv
	 jjx=xj+50
         if(jjx.gt.100.or.jjx.le.0) goto 5
	 fvx(jjx)=fvx(jjx)+1
   5     yj=vy(i)/dv
	 jjy=yj+50
         if(jjy.gt.100.or.jjy.le.0) goto 6
         fvy(jjy)=fvy(jjy)+1
   6     continue
       end  do
       write(5,*) it
       do j=1,100
          write(5,10) (j-50)*dv,fvx(j),fvy(j)
       end do
c       mfv=fvx(1)
c       do j=1,100
c       mfv=max0(fvx(j),fvy(j),mfv) 
c       end do
       mfv=2500
       call setviewport(300,150,600,250)
       dummy=setwindow(.true.,0,100,600,0)
       call clearscreen($gviewport)
       
       dummy=setcolor(9)
       dummy=rectangle_w($gborder,0,100,600,0)  
c       call _fq_settextposition(20,10,xy)
       call  moveto_w(0.0d0,fvx(1)*100.0d0/mfv,wxy)
	do j=1,100
	dummy=lineto_w((j-1)*600.0d0/100.0d0,fvx(j)*100.0d0/mfv)
	end do
         dummy=setcolor(9)      
	call  moveto_w(0.0d0,fvy(1)*100.0d0/mfv,wxy)
	do  j=1,100
	  dummy=lineto_w((j-1)*600.0d0/100.0d0,fvy(j)*100.0d0/mfv)
	end  do
c       DO i=io,in
c         vx(i)=vx(i)/dxdt
c         vy(i)=vy(i)/dxdt
c       end do
  10    format(f8.4,2i5)
        return
	end
c##########################################################       
	
	
c###############################################################        
*****************************************************************************       
c       subroutine his(it,ihis,ke,ese,te)
****************************************************************************
c       include 'fgraph.fd'
c       dimension ke(4000),ese(4000),te(4000)
c       integer*2 dummy
c       record /xycoord/ xy  
c       record /wxycoord/ wxy
c       integer it,ihis
c       real ke,ese,kem,esem,TEM
c       if((it/ihis)*ihis.ne.it.or.it.eq.0) return
c       call setviewport(0,250,300,350)
c       dummy=setwindow(.true.,0,100,500,0)
c       call clearscreen($gviewport)
c       dummy=setcolor(9)
c       dummy=rectangle_w($gborder,0,100,500,0)  
c       call moveto_w(0,0,wxy)
c       dummy=setcolor(12)
c       kem=ke(1)
c       esem=ese(1)
C ****       TEM=MAXIMUM ELECTRON TEMPERATURE 
c       TEM=te(1)
c       do 2 i=1,it
c       if(kem.le.ke(i)) kem=ke(i)
c       if(TEM.le.te(i)) tem=tE(I)
c   2   if(esem.le.ese(i)) esem=ese(i)
c       call _fq_settextposition(10,20,xy)
c       do 10 i=1,it 
c 10    dummy=lineto_w(i*500/4000,ke(i)*100*.9/kem)    
c       call moveto_w(0,0,wxy)
c       dummy=setcolor(7)
c       do 11 i=1,it 
c11     dummy=lineto_w(i*500/4000,ese(i)*100*.9/esem)
c       call setviewport(0,350,600,450)
c       dummy=setwindow(.true.,0,100,500,0)
c       call clearscreen($gviewport)
c       DUMMY=setcolor(9)
c       dummy=rectangle_w($gborder,0,100,500,0)  
c       call moveto_w(1*500/4000,te(1)*100*.9/tem,wxy)
c       dummy=setcolor(11)
c       do i=1,it 
c       dummy=lineto_w(i*500/4000,Te(i)*100*.9/TEM)
c       END DO
cc       return 
c      end
*****************************************************************************       
       subroutine lasf(it,eyyn,EY1,bbz,ngx)
****************************************************************************
       include 'fgraph.fd'
       integer*2 dummy
c       record /xycoord/ xy  
       record /wxycoord/ wxy
       integer it,ngx
       real eyyn(34),bbz(34),bbzm,eyynm,EY1(34)
c       if((it/ihis)*ihis.ne.it.or.it.eq.0) return
C       call setviewport(0,0,600,450)
C       call clearscreen($gviewport)
       call setviewport(0,200,300,300)  
       dummy=setwindow(.true.,0,100,500,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       dummy=rectangle_w($gborder,0,100,500,-100)  
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=lineto_w(500,0)
       dummy=setcolor(9)
       eyynm=abs(eyyn(1))
       bbzm=(bbz(1))
       do 2 j=1,ngx
       eyynm=Amax1(eyynm,abs(eyyn(j)))
  2    bbzm=Amax1(bbzM,abs(bbz(j)))
       CALL MOVETO_W(500.0d0/NGX,EYYN(1)*100.0d0*.9/EYYNM,WXY)
       do 10 j=1,ngx 
 10    dummy=lineto_w(j*500.0d0/ngx,eyyn(j)*100.0d0*.9d0/eyynm)    
       CALL MOVETO_W(500.0d0/NGX,EY1(1)*100.0d0*.9d0/EYYNM,WXY)
C       dummy=setcolor(14)
C       do 11 j=1,ngx 
C 11    dummy=lineto_w(j*500/ngx,ey1(j)*100*.9/eyynm)    
       if(bbzm.le.0) bbzm=1.
       call setviewport(0,300,300,400)
       dummy=setwindow(.true.,0,100,500,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       dummy=rectangle_w($gborder,0,100,500,-100)  
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=lineto_w(500,0)
       dummy=setcolor(9)
       CALL MOVEto_w(1d0*500.0d0/ngx,bbz(1)*100.0d0*.9d0/bbzm,WXY)  
       do 12 j=1,ngx 
 12    dummy=lineto_w(j*500.0d0/ngx,bbz(j)*100.0d0*.9d0/bbzm)    
       
       return 
       end
*****************************************************************************       
       
       subroutine pfor(ek,ngx)
****************************************************************************
       include 'fgraph.fd'
       dimension ek(ngx)
       integer*2 dummy
       record /wxycoord/ wxy  
       integer ngx,j
       real ek,ekm
c       if((it/ihis)*ihis.ne.it) return
	call setviewport(300,250,600,350)
       dummy=setwindow(.true.,0,100,500,-100)
       call clearscreen($gviewport)
       dummy=setcolor(9)
       dummy=rectangle_w($gborder,0,100,500,-100)  
       call moveto_w(0.0d0,0.0d0,wxy)
       dummy=setcolor(12)
       ekm=abs(ek(1))
	 print*,'it',it  
       do j=1,ngx
	 if(ekm.le.abs(ek(j))) ekm=abs(ek(j))
       end do
       if(ekm.le.0) return 
       do  j=1,ngx 
	   dummy=lineto_w(j*500.0d0/ngx,ek(j)*100.0d0*.8d0/ekm)    
       end do
       return 
       end
c########################################################
	subroutine cpft(r,i,n,incp,signp)
************************************************************             
	     real r(1), i(1)
	     integer signp,span,rc
	     real sines(15),i0,i1
	     data sines(1)/0./
	     if(sines(1).eq.1) go to 1
	     sines(1)=1.
	     t=atan(1.)
	     do 2 is=2,15
	     sines(is)=sin(t)
   2         t=t/2.
   1         continue
	   if(n.eq.1) return
	   inc=incp
	   sgn=signp
	   ninc=n*inc
	   span=ninc
	   it=n/2
	   do 3 is= 1,15
	   if(it.eq.1) go to 12
   3       it=it/2
 10        t=s+(s0*c-c0*s)
	   c=c-(c0*c+s0*s)
	   s=t
 11        k1=k0+span
	   r0=r(k0+1)
	   r1=r(k1+1)
	   i0=I(k0+1)
	   i1=i(k1+1)
	   r(k0+1)=r0+r1
	   i(k0+1)=i0+i1
	   r0=r0-r1
	   i0=i0-i1
	   r(k1+1)=c*r0-s*i0
	   i(k1+1)=s*r0+c*i0
	   k0=k1+span
	   if(k0.lt.ninc) go to 11
	   k1=k0-ninc
	   c=-c
	   k0=span-k1
	   if(k1.lt.k0) go to 11
	   k0=k0+inc
	   if(k0.lt.k1) go to 10
   12      continue 
	   span=span/2
	   
	   k0=0
  13        k1=k0+span
	    r0=r(k0+1)
	    r1=r(k1+1)
	    i0=i(k0+1)
	    i1=i(k1+1)
	    r(k0+1)=r0+r1
	    i(k0+1)=i0+i1
	    r(k1+1)=r0-r1
	    i(k1+1)=i0-i1
	    k0=k1+span
	    if(k0.lt.ninc)go to 13
	    if(span.eq.inc)go to 20
	    c0=2.*sines(is)**2
	    is=is-1
	    s=sign(sines(is),sgn)
	    s0=s
	    c=1.-c0
	    k0=inc
	    go to 11
   20    n1=ninc-inc
	n2=ninc/2
	   rc=0
	  ij=0
	  ji=0
	  if(n2.eq.inc)return
	 go to 22
  21     ij=n1-ij
	 ji=n1-ji
	 t=r(ij+1)
	 r(ij+1)=r(ji+1)
	  r(ji+1)=t
	  t=i(ij+1)
	  i(ij+1)=i(ji+1)
	  i(ji+1)=t
	  if(ij.gt.n2)go to 21
  22      ij=ij+inc
	  ji=ji+n2
	  t=r(ij+1)
	  r(ij+1)=r(ji+1)
	  r(ji+1)=t
	  t=i(ij+1)
	  i(ij+1)=i(ji+1)
	  i(ji+1)=t
	  it=n2
 23       it=it/2
	   rc=rc-it
	if(rc.ge.0) go to 23
	rc=rc+2*it
	ji=rc
	ij=ij+inc
	if(ij.lt.ji)go to 21
	if(ij.lt.n2)go to 22
	return
	end 
c################ end of cpft ############################
****************************************************************        
	subroutine rpft2(a,b,n,incp)
*****************************************************************      
	real a(1),b(1)
	  real ip,im
	inc=incp
	ninc=n*inc
	a(1)=a(1)+A(1)
	b(1)=b(1)+b(1)
	lp=inc
	lm=ninc-lp
	if(lp.ge.lm) go to 2
 1       rp=a(lp+1)
	 rm=a(lm+1)
	 ip=b(lp+1)
	 im=b(lm+1)
	 a(lp+1)=rm+rp
	 b(lm+1)=rm-rp
	 b(lp+1)=ip+im
	 a(lm+1)=ip-im
	 lp=lp+inc
	 lm=ninc-lp
	 if(lp.lt.lm)go to 1
  2       if(lp.gt.ninc)return
	  a(lp+1)=a(lp+1)+a(lp+1)
	  b(lp+1)=b(lp+1)+b(lp+1)
	  return
	  end
c############# end of rpft2 ###########################                 
************************************************************
	subroutine rpfti2 (a,b,n,incp)
*************************************************************        
	real a(1),b(1)
	inc=incp
	ninc=n*inc
	lp=inc
	lm=ninc-lp
	if(lp.ge.lm)return
 3       ca=a(lp+1)
	sb=b(lm+1)
	cb=b(lp+1)
	sa=a(lm+1)
	a(lp+1)=ca-sb
	a(lm+1)=ca+sb
	b(lp+1)=cb+sa
	b(lm+1)=cb-sa
	lp=lp+inc
	lm=ninc-lp
	if(lp.lt.lm)go to 3
	return
	end
c######### eof rpfti2 ##################################3