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 !--------------------------------------------------------