cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c THIS SUBROUTINE, GIVEN A PARAMETRIZATION BETA(J)=j*h, GIVES A c SET OF EQUALLY SPACED POINTS IN ARCLENGTH c THIS IS DONE VIA NEWTON ITERATION c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine arc(m,h,ds,alpha,tol) implicit double precision (a-h,o-z) dimension alpha(0:8192),dsi(0:8192) dimension ds(0:8192),f(0:8192),fp(0:8192) dimension a(0:8192),b(0:8192),c(0:8192),d(0:8192) m1=m+1 pi=4.d0*datan(1.d0) p2=2.d0*pi c integrate ds and set-up functions for interpolation and iteration call fin(m,h,p2,dsi,ds) ave=dsi(m) do j=0,m f(j)=dble(j)*h*ave-dsi(j) fp(j)=-ds(j) end do do j=0,m a(j)=f(j)/dble(m) b(j)=0.d0 c(j)=fp(j)/dble(m) d(j)=0.d0 end do call fft842(0,m,a,b) call fft842(0,m,c,d) c for each point, iterate. alpha(0)=0.d0 do j=1,m ao=alpha(j-1)+h*(ave/ds(j)) c begin iterations do l=0,100 aoo=p2*ao call synth(f,v,aoo,m1,a,b) call synth(fp,vp,aoo,m1,c,d) v=v+(dble(j)*h-ao)*ave dalp=v/vp an=ao-dalp err=dabs(dalp) c write(6,*)dalp if (err .le. tol) goto 1000 ao=an end do 1000 continue alpha(j)=an end do return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c GIVEN AN ARCLENGTH SL, AND A TAN ANGLE THETA, THIS CODE c RECONSTRUCTS THE CURVE AND GIVES X,Y c c NEED EQUAL ARCLENGTH TO DO c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine recon (m,h,p2,sl,theta,x,y,avx,avy,dcut) implicit double precision (a-h,o-z) dimension theta(0:8192),x(0:8192),y(0:8192) dimension temp(0:8192),temp1(0:8192) call kfilter(m,dcut,theta) do j=0,m-1 temp(j)=sl*cos(theta(j)) temp1(j)=sl*sin(theta(j)) end do call fin(m,h,p2,x,temp) call fin(m,h,p2,y,temp1) avx=x(m) avy=y(m) c avx=1.d0 c avy=0.d0 do j=0,m x(j)=x(j)-dble(j)*h*avx y(j)=y(j)-dble(j)*h*avy end do call kfilter(m,dcut,x) call kfilter(m,dcut,y) do j=0,m x(j)=x(j)+dble(j)*h y(j)=y(j) end do avy=0.d0 avx=1.d0 return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c THIS CODE FORMS THE INDEFINITE INTEGRAL OF A PERIODIC FUNCTION c c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine fin (n,h,p2,yi,y) implicit double precision (a-h,o-z) dimension yi(0:8192),y(0:8192),a(0:8192) dimension b(0:8192) do j=0,n-1 a(j)=y(j) end do a(n)=a(0) call fast(a,n) b(0)=0.d0 b(1)=0.d0 do j=1,n/2 k=2*j b(k)=a(k+1)/(p2*dble(j)) b(k+1)=-a(k)/(p2*dble(j)) end do call fsst(b,n) b(n)=b(0) do j=0,n yi(j)=(a(0)/dble(n))*dble(j)*h + b(j) - b(n) end do return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c THIS FILE CONTAINS kfilter.f kfilterx.f fd1.f acurv.f c kfiltern.f c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c fourier filter c with krasny filtering with 1 input vector c c USING FAST.f FOURIER TRANSFORM c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine kfilter(n,dcut,u) implicit double precision(a-h,o-z) dimension u(0:8192) dimension b(0:8192) do j=0,n-1 b(j)=u(j) end do b(n)=b(0) c filter fourier coeffs call fast(b,n) test=(b(0)/dble(n))**2 +(b(1)/dble(n))**2 test=dsqrt(test) if (test .le. dcut) then b(0)=0.d0 b(1)=0.d0 end if do j=1,n/2 k=2*j test=(b(k)/dble(n))**2 + (b(k+1)/dble(n))**2 test=dsqrt(test) if (test .le. dcut) then b(k)=0.d0 b(k+1)=0.d0 end if end do do k=n-3,n+1 b(k)=0.d0 end do C C C Applying Fourier smoothing of 25th order C do k=1,n b(k)=b(k)*exp(-10.d0*(dble(k)/dble(n))**25) end do C call fsst(b,n) do j=0,n-1 u(j)=b(j) end do u(n)=b(0) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c kfilter x c filters a vector = periodic + linear c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine kfilterx(m,dcut,x,h) implicit double precision (a-h,o-z) dimension x(0:8192) do i=0,m x(i)=x(i)-dble(i)*h end do call kfilter(m,dcut,x) do i=0,m x(i)=x(i)+dble(i)*h end do return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c kfiltern x c filters N/2 mode of a vector = periodic + linear c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine kfilternx(m,x,h) implicit double precision (a-h,o-z) dimension x(0:8192) do i=0,m x(i)=x(i)-dble(i)*h end do call kfiltern(m,x) do i=0,m x(i)=x(i)+dble(i)*h end do return end c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c this routine filters the n/2 mode only c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine kfiltern(n,u) implicit double precision (a-h,o-z) dimension u(0:8192) dimension b(0:8192) do j=0,n-1 b(j)=u(j) end do c remove n/2 coeff call fast(b,n) b(n)=0.d0 b(n+1)=0.d0 call fsst(b,n) do j=0,n-1 u(j)=b(j) end do u(n)=u(0) return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c subroutine to compute SPECTRAL derivatives in space c for periodic data c c c c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine fd1(n,x,sx,h,dcut) implicit double precision(a-h,o-z) dimension x(0:8192),sx(0:8192),x1(0:8192),sx1(0:8192) call kfilter(n,dcut,x) do i=0,n-1 x1(i)=x(i) end do x1(n)=x(0) call fdiff(x1,sx1,n,h) do i=0,n-1 sx(i)=sx1(i) end do sx(n)=sx(0) call kfilter(n,dcut,sx) return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c THIS ROUTINE COMPUTES THE CURVATURE FROM THE TANGENT c ANGLE FORMULATION c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine acurv(m,sl,theta,dkap,dcut,h) implicit double precision (a-h,o-z) dimension theta(0:8192),dkap(0:8192) dimension stheta(0:8192) call fd1(m,theta,stheta,h,dcut) call kfilter(m,dcut,stheta) do j=0,m-1 dkap(j)=stheta(j)/sl end do dkap(m)=dkap(0) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c THIS FILE CONTAINS ajacoby.f rtvel.f rnvel.f rlen.f vel.f c c c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c THIS ROUTINE COMPUTES THE TANGENT VELOCITY GIVEN THE NORMAL c VELOCITY. IT ASSUMES THAT THE TANGENTIAL VELOCITY AT ALPHA=0 c IS 0 c c IT IS TO BE USED IN CONJUNCTION WITH THE TANGENT ANGLE FORMULATION c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine rtvel(m,h,p2,theta,unorm,utan,dcut) implicit double precision (a-h,o-z) dimension theta(0:8192),unorm(0:8192),utan(0:8192) dimension stheta(0:8192),temp(0:8192) call fd1(m,theta,stheta,h,dcut) call kfilter(m,dcut,stheta) do j=0,m-1 temp(j)=stheta(j)*unorm(j) utan(j)=0.d0 end do utan(m)=0.d0 temp(m)=temp(0) call fin(m,h,p2,utan,temp) do j=0,m-1 utan(j)=utan(j)-utan(m)*dble(j)*h end do utan(m)=utan(0) call kfilter(m,dcut,utan) return end