subroutine initialize(Nx,Ny,x,y,f)
      implicit none
      integer Nx,Ny
      real*8 x(Nx),y(Ny),f(Nx,Ny),new
      
      integer i,j
      real*8 r2

c

c      write(*,*) "dentro init xc1",xc(1),"xc2",xc(2)
c      do i=1,Nx
c      write(*,*)"dentrox",i,x(i)
c      enddo
c      do j=1,Ny
c            write(*,*)"dentroy",j,y(j)
c                  enddo

	new=21.0d0

c
        do j=1,Ny
        do i=1,Nx
c          r2 = (x(i)-xc(1))**2 + (y(j)-xc(2))**2
c          f(i,j) = exp(-r2/delta**2)
           f(i,j)=1.0d0
        enddo
      enddo
c	write(*,*)'texas new',new

      return
      end

      subroutine calc_center(time,arm,x0,y0,xc)
      implicit none
      real PI
      parameter (PI=3.1415926535897932385d0)
      real*8 time,arm,x0,y0,xc(*)
      xc(1) = x0
c+ arm*cos(2.d0*PI*time)
      xc(2) = y0
c + arm*sin(2.d0*PI*time)
c      write(*,*) "xc1",xc(1),"xc2",xc(2)
      return
      end


      subroutine newposition(Nx,Ny,pn,time,pposx,pposy,x,y)
      implicit none
c      real*8 x(*),y(*)
      integer Nx,Ny,pn,i,j
      real*8 x(Nx),y(Ny)
      real PI,count
      parameter (PI=3.1415926535897932385d0)
      real*8 time,pposx(*),pposy(*)
c      vel = 0.01d0
      do i=1,Nx
      pposx(2)=pposx(2)+x(i)
c      write(*,*)"pposx",i,x(i),pposx(1)
      enddo
      do j=1,Ny
      pposy(2)=pposy(2)+y(i)
c            write(*,*)"pposy",j,y(j),pposy(1)
                  enddo

c      do i=1,pn
c                pposx(i) = pposx(i)+time
c                pposy(i) = pposy(i)+time
c                write(*,*) out(i),"pn",pn
c      enddo
      return
      end

c-----------------------------------------------------------

      subroutine save_particles(time,savex,savey,Np,rank)

      implicit none
      integer i,rank,Np,tag1,tag2
      real*8 time
      real*8 savex(Np), savey(Np)
      real*8 scratch(4*Np)
      character*200 fname
      character*5 intoc,ran

      fname = "2d_"
      ran = intoc(rank)
      tag1 = index(fname,' ')
      tag2 = index(ran,' ')
      fname = fname(1:tag1-1)//ran(1:tag2-1)

      do i=1,Np
      scratch(i) = savex(i)
      scratch(Np+i) = savey(i)
      scratch(2*Np+i) = 0.d0
      scratch(3*Np+i) = 0.d0
      enddo

      call gft_out_part_xyzf(fname,time,scratch,Np,1)
      return
      end

c-------------------------------------------------------------
      subroutine make_step(dt,time,px,py,vtype,Nx,Ny,x,y,vmx,vmy)

      implicit none
      integer i,j,k,vtype,Nx,Ny
      real*8 time,px,py,pxx,pyy,vx,vy,pi,dt,vxx,vyy
      real*8 x(Nx), y(Ny), vmx(Ny), vmy(Ny)
      real*8 vxg,vyg
c       Define pi
      pi = 4.0d0*atan(1.0d0)
c      write(*,*)'vtype=',vtype

	if (vtype .eq. 1) then

c===========================================
c  three ways to calculate the velocity 
c  selected by vtype
c============================================
c     angular vel=2*pi

c     first RK step
      vx=2.0d0*pi*(-py)
      vy=2.0d0*pi*(px)

      pxx=px+vx*dt
      pyy=py+vy*dt

c     Second RK step

      vxx=2.0d0*pi*(-pyy)
      vyy=2.0d0*pi*(pxx)
      px=0.5d0*(px+pxx+vxx*dt)
      py=0.5d0*(py+pyy+vyy*dt)


	else if (vtype .eq. 2) then

c=============================================
c     angular velocity= 2*pi*r

      vx=2*pi*(-py)*sqrt(px*px+py*py)
      vy=2*pi*(px)*sqrt(px*px+py*py)

      pxx=px+vx*dt
      pyy=py+vy*dt

      vxx=2*pi*(-pyy)*sqrt(pxx*pxx+pyy*pyy)
      vyy=2*pi*(pxx)*sqrt(pxx*pxx+pyy*pyy)
      px=0.5d0*(px+pxx+vxx*dt)
      py=0.5d0*(py+pyy+vyy*dt)

      else

c===========================================================
c       Grid defined velocity: find i,j of my px,py

c-----------------------------------------------------------
c       First we determine i and j

	i=1
	do while (x(i) .lt. px)
	
	write(*,*) 'today Cx i=',i-1,'x=',x(i),'px=',px	
        
	i=i+1
	enddo
	write(*,*)'today i Cx=',i-2 
c       one minus starting from 1, one minus the extra 
c       i=i+1

        j=1
        do while (y(j) .lt. py)

        write(*,*) 'today Cy j=',j-1,'y=',y(j),'py=',py

        j=j+1
        enddo
        write(*,*)'today j Cy=',j-2
c       one minus starting from 1, one minus the extra
c       j=j+1

c-----------------------------------------------------------
c       First Linear interpolation
        i=i-1
	j=j-1
	
	vxg=vmx(j)+((vmx(j+1)-vmx(j))/(y(j+1)-y(j)))*(py-y(j))
        write(*,*)'j=',j,'vmx(j)=',vmx(j),'vmx(j+1)='
     &  ,vmx(j+1),'vxg=',vxg

        vyg=vmy(i)+((vmy(i+1)-vmy(i))/(x(i+1)-x(i)))*(px-x(i))
        write(*,*)'i=',i,'vmy(i)=',vmy(i),'vmy(i+1)='
     &  ,vmy(i+1),'vyg=',vyg

        vx=vxg
        vy=vyg


c     first RK step

      pxx=px+vx*dt
      pyy=py+vy*dt


c---------------------------------------------------------
c       First we determine i and j
        i=1
        do while (x(i) .lt. pxx)

        i=i+1
        enddo
c       one minus starting from 1, one minus the extra
c       i=i+1

        j=1
        do while (y(j) .lt. pyy)

        j=j+1
        enddo
c       one minus starting from 1, one minus the extra
c       j=j+1
c-------------------------------------------------------
c     Second Linear interpolation
	i=i-1
	j=j-1
	write(*,*)' second calc'
        vxg=vmx(j)+((vmx(j+1)-vmx(j))/(y(j+1)-y(j)))*(pyy-y(j))
        write(*,*)'j=',j,'vmx(j)=',vmx(j),'vmx(j+1)='
     &  ,vmx(j+1),'vxg=',vxg

        vyg=vmy(i)+((vmy(i+1)-vmy(i))/(x(i+1)-x(i)))*(pxx-x(i))
        write(*,*)'i=',i,'vmy(i)=',vmy(i),'vmy(i+1)='
     &  ,vmy(i+1),'vyg=',vyg

        vxx=vxg
        vyy=vyg





c     Second RK step

      px=0.5d0*(px+pxx+vxx*dt)
      py=0.5d0*(py+pyy+vyy*dt)



	endif
      return
      end