subroutine initia(m,phi,dcut,tol,xn,yn)
*     *****************************************************************
*     initialize the curve by giving a set of points and reparametrize
*     so they are equidistanced between every consecutive points.
*     the fluctuation part of tangent angle (phi) is also generated.
*     The initial length of the curve is sa*p2.
*     In this code, a closed curve based on
*        ro(th) = 1 + eps*(sin(5*al)+cos(4*al))    0 < al < 2pi
*     is used.
*     ******************************************************************
      integer m,j
      double precision xn(0:8192),yn(0:8192),sxn(0:8192),syn(0:8192)
      double precision slo(0:8192),x(0:8192),y(0:8192)
      double precision sx(0:8192),sy(0:8192),phi(0:8192)
      double precision alpha(0:8192)
      double precision pi,p2,h,hp,sa,ax,bx,dcut,tol,ta,pert,xr,yr,
     +                 tnext,tprev,tj,ra,rad,pt,tant,angle
      parameter( ax=1.5D0, bx=1.D0, pert=0.4D0 )
      common/inicur/rad,pt
      common/salpha/sa
      common/gridsz/h,hp
      common/valupi/pi,p2

      do 10 j=0,m
         ta = dfloat(j)*h*p2
         x(j) = xr(ta)
         y(j) = yr(ta)
 10   continue

      call fd1(m,x,sx,h,dcut)
      call fd1(m,y,sy,h,dcut)

      do 20 j=0,m-1
         slo(j) = dsqrt(sx(j)**2 + sy(j)**2)
 20   continue
      slo(m) = slo(0)

      call arc(m,h,slo,alpha,tol)

*     lay down the new points and translate to a coordinate system
*     where (0.5,0.5) is the center of the box.
      do 30 j=0,m
         ta    = p2*alpha(j)
         xn(j) = xr(ta) + 0.5D0
         yn(j) = yr(ta) + 0.5D0
 30   continue

      call fd1(m,xn,sxn,h,dcut)
      call fd1(m,yn,syn,h,dcut)

      sa = 0.D0
      do 50 j=0,m-1
         sa = sa + dsqrt(sxn(j)**2 + syn(j)**2)
 50   continue
      sa = sa/dble(m)/p2
      print *,'sa=',sa

*     generate phi from the arc tangent of y'/x':
*     subtract by alpha to make phi a periodic function (the correction
*     part).

      phi(0) = angle(sxn(0),syn(0))
      tj = 0.D0
      tprev = phi(0)
      do 60 j=1,m
         tnext = datan(syn(j)/sxn(j))
         if ( dabs(tnext-tprev) .ge. 1.D0 ) tj = tj + pi
         phi(j) = tnext + tj - p2*h*dfloat(j)
         tprev = tnext
 60   continue

      call kfilter(m,dcut,phi)

 101  format(2(f12.6,2x))
      end

      function angle(sx,sy)
      double precision sx,sy,tant,angle,pi,p2
      common/valupi/pi,p2
      tant = sy/sx
      angle = datan(tant)
      if ( tant .ge. 0.D0 .and. sy .lt. 0.D0 ) 
     +    angle = angle + pi
      if ( tant .lt. 0.D0 .and. sy .gt. 0.D0 )
     +    angle = angle + pi
      end

      function xr(al)
      double precision al,xr,rad,pt
      common/inicur/rad,pt
      xr = (rad + pt*(dsin(5.D0*al) + dcos(4.D0*al)))*dcos(al)
      end

      function yr(al)
      double precision al,yr,rad,pt
      common/inicur/rad,pt
      yr = (rad + pt*(dsin(5.D0*al) + dcos(4.D0*al)))*dsin(al)
      end