      program curves
! Curvilinear motion along a curve
      implicit double precision (a-h,o-z)
      common /par/dt,stol,xm,aa,a2,bb,b2,y0,xnse,nm,nl,shape
      common /pos/x,y,xo,yo,xn,yn,vx,vy
      character*6 shape
!
      open(unit=7,file='curves.out',form='formatted'
     1 ,status='unknown')
      open(unit=77,file='curves.dat',form='formatted'
     1 ,status='unknown') ! for gnuplot
! Input deck:
      nt=2000;    ! no of time steps 
      nl=80;      ! max no of Lagrange iterations
      dt=0.01d0; dth=0.5d0*dt;
      dt2=dt*dt; dtsh=0.5d0*dt2;
      stol=1.d-12; ! sigma tolerance
      xm=1.d0;
! Diverse shape parameters:
      aa=1.d0; a2=aa*aa; ! for circle, ellipse, disorectangle, superellipse
      bb=2.d0; b2=bb*bb; ! for ellipse, disorectangle, superellipse
      y0=0.8d0*bb;       ! for discorectangle
      xnse=2.5d0;        ! for superellipse (Piet Hein: xnse=2.5)
!
      vx=0.d0; vy=0.5d0;
!
!      shape='circle'
       shape='ellips'
!      shape='dscrct'
!      shape='supell'
!      shape='beancv'
!
      call startconf
!

      vx=0.5d0*(xn-xo)/dt; vy=0.5d0*(yn-yo)/dt; 
      vm=dsqrt(vx*vx+vy*vy);
      ek=0.5d0*xm*(vx*vx+vy*vy); 
      write(*,100) it,xn,yn,ek,ep,et,am
!
      do 1 it=1,nt
!  Verlet move:
        xn=2.d0*x-xo; yn=2.d0*y-yo
!  Lagrange-1 procedure:
        call guiseppe 
!        write(*,199) x,y,sig,xn,yn
199     format(1x,'curves: x_y,xn_yn = ',5e14.5)
        vx=0.5d0*(xn-xo)/dt; vy=0.5d0*(yn-yo)/dt;
        vm=dsqrt(vx*vx+vy*vy);
        ax=(xn-2.d0*x+xo)/dt2; ay=(yn-2.d0*y+yo)/dt2;
        am=dsqrt(ax*ax+ay*ay);
        ek=0.5d0*xm*(vx*vx+vy*vy); 
        if(it.gt.4) then 
          write(*,100) it,xn,yn,ek,am
          write(7,100) it,xn,yn,ek,am
          write(77,100) it,xn,yn,ek,am
100       format(1x,'it,xn,yn,ek,am= ',i6,7e11.3)
        endif
        xo=x; yo=y;
        x=xn; y=yn;
!
1     continue   ! end time step
!
      stop
      end
!
!------------------------------------------------------
      subroutine startconf
      implicit double precision (a-h,o-z)
      common /par/dt,stol,xm,aa,a2,bb,b2,y0,xnse,nm,nl,shape
      common /pos/x,y,xo,yo,xn,yn,vx,vy
      character*6 shape
!
      x=aa; if(shape.eq.'beancv') x=1.d0; 
      y=0.d0;
! Determine r(-1):
      xn=x-dt*vx; yn=y-dt*vy; call guiseppe; 
      xo=xn; yo=yn;  ! r(-1) for Verlet
      xa=0.d0; ya=0.d0; ! just initialized; will be computed in the nt loop  
      return
      end
!
!------------------------------------------------------
      subroutine guiseppe
      implicit double precision (a-h,o-z)
      common /par/dt,stol,xm,aa,a2,bb,b2,y0,xnse,nm,nl,shape
      common /pos/x,y,xo,yo,xn,yn,vx,vy
      character*6 shape
      do is=1,nl   ! Loop over Lagrange iterations
        call sigma(sig,dsx,dsy,dsxn,dsyn)
        sr=dsxn*dsx+dsyn*dsy;
        cg=sig/sr/4.d0;
        xn  =xn  -cg*dsx;
        yn  =yn  -cg*dsy;
        if(dabs(sig).lt.stol) exit
      enddo       ! end Lagrange loop
!
      return
      end
!
!------------------------------------------------------
      subroutine sigma(sig,dsx,dsy,dsxn,dsyn)
      implicit double precision (a-h,o-z)
      common /par/dt,stol,xm,aa,a2,bb,b2,y0,xnse,nm,nl,shape
      common /pos/x,y,xo,yo,xn,yn,vx,vy
      character*6 shape
 !     print *,shape
      if(shape.eq.'circle') then      ! circle
        sig=(xn*xn+yn*yn)/xl2-1.d0;
        dsx=2.d0*x/xl2; dsy=2.d0*y/xl2;
        dsxn=2.d0*xn/xl2; dsyn=2.d0*yn/xl2;
      else if(shape.eq.'ellips') then ! ellipse
        sig=xn*xn/a2+yn*yn/b2-1.d0
        dsx=2.d0*x/a2; dsy=2.d0*y/b2; 
        dsxn=2.d0*xn/a2; dsyn=2.d0*yn/b2;
      else if(shape.eq.'supell') then ! superellipse
        termx=dabs(xn/aa)**xnse;
        termy=dabs(yn/bb)**xnse;
        sig=termx+termy-1.d0;
        sx=dsign(1.d0,x);   sy=dsign(1.d0,y);
        dsx=(xnse/aa)*sx*dabs(x/aa)**(xnse-1.d0);
        dsy=(xnse/bb)*sy*dabs(y/bb)**(xnse-1.d0);
        sxn=dsign(1.d0,xn);   syn=dsign(1.d0,yn);
        dsxn=(xnse/aa)*sxn*dabs(xn/aa)**(xnse-1.d0);
        dsyn=(xnse/bb)*syn*dabs(yn/bb)**(xnse-1.d0);
      else if(shape.eq.'dscrct') then ! discorectangle
        if(yn.le.y0) then
          sig=(xn*xn/a2)-1.d0;
          dsxn=2.d0*xn/a2; dsyn=0.d0;
        else
          sig=(xn*xn+(yn-y0)*(yn-y0))/a2-1.d0;
          dsxn=2.d0*xn/a2; dsyn=2.d0*sign(1.d0,yn)*(dabs(yn)-y0);
        endif
        if(y.le.y0) then
          dsx=2.d0*x/a2; dsy=0.d0;
        else
          dsx=2.d0*x/a2; dsy=2.d0*sign(1.d0,y)*(dabs(y)-y0);
        endif
      else if(shape.eq.'beancv') then ! bean curve
        xn2=xn*xn; yn2=yn*yn;
        sig=xn2*xn2+xn2*yn2+yn2*yn2-xn*(xn2+yn2);
        dsx=4.d0*x*x*x-3.d0*x*x+2.d0*x*y*y-y*y;
        dsy=4.d0*y*y*y+2.d0*y*(x*x-x);
        dsxn=4.d0*xn*xn*xn-3.d0*xn*xn+2.d0*xn*yn*yn-yn*yn;
        dsyn=4.d0*yn*yn*yn+2.d0*yn*(xn*xn-xn);
      endif
!
      return
      end
!--------------------------------------------------------

