FORTRAN Program

Reference

      Program wedge
c 
c     Copyright (1996) belongs to Harris Wong.
c     Use this program at your own risk.
c
c     This program solves the self-similar film profile of a film
c        wedge by a fourth order Runge Kutta method.
c     Given a contact angle (cangle) and a position of the contact 
c        line (x0), the program iterates on the curvature (c0) at the
c        contact line until it finds the wedge angle alpha.
c
c     Input:	
c		n	: the number of uniform steps
c  		ds	: displacement of each step
c		cangle	: contact angle in degrees
c		x0	: the position of the contact line
c
c     	 Note: 	n*ds = the domain length
c
c     Output:
c		(x, y) location of the film profile
c
c
c
      implicit double precision (a-h, o-z)
      real*8 x(0:50000), y(0:50000)
      open (unit=8, file='y.in', status='unknown')
      open (unit=9, file='y.out', status='unknown')
      open (unit=10, file='shape.out', status='unknown')
      pi = dasin(1d0)*2d0
c
c     Read in the controlling parameters
c
      read(8,*) n, ds
      read(8,*) cangle, x0

      gtod = 180d0/pi
      beta = cangle/gtod
      write(9,990) n, ds
      write(9,*) cangle
c
c     Compute alpha iteratively
c
      c00 = -20d0
      c01 = 5d0
      do i=1, 60
         c0 = (c00 + c01)/2d0
         call profile(n,ds,beta,x0,c0,x,y,alpha,c1)
         if (c1.lt.0d0) then
            c00 = c0
         else
            c01 = c0
         end if
         write(9,*) i, c0, c1, alpha*gtod
         if (dabs(c1).lt.1d-6) goto 100
      end do
  100 continue
      write(9,*) i, x0, c0, alpha*gtod
c
c     Print the (x, y) location of the film profile
c
      do i=0, n, n/200
         write(10,*) x(i), y(i)
      end do
      write(10,*) 

  990 format(i6, 5f18.6)
  991 format(5e15.6)
      stop
      end


      subroutine profile(n, ds, beta, x00, c0, x, y, alpha, c1)
c
c     This subroutine uses a fourth order Runge Kutta method to
c        computes the film profile given the initial
c        x position (x00) and the curvature (c0).
c
      implicit double precision (a-h, o-z)
      real*8 k1,k2,k3,k4,l1,l2,l3,l4,m1,m2,m3,m4,n1,n2,n3,n4 
      real*8 q1,q2,q3,q4, a(4), x(0:50000), y(0:50000)

      s = 0d0
      x0 = x00
      y0 = 0d0
      p0 = beta
      u0 = c0
      v0 = 0d0
      x(0) = x0
      y(0) = y0
c
c     compute p(s), x(s), y(s), u(s), and v(s)
c
      do 10 i=1, n
         k1 = dcos(p0)*ds
         l1 = dsin(p0)*ds
         m1 = u0*ds
         n1 = v0*ds
         q1 = g(x0, y0, p0)*ds
         k2 = dcos(p0 + m1/2d0)*ds
         l2 = dsin(p0 + m1/2d0)*ds
         m2 = (u0 + n1/2d0)*ds
         n2 = (v0 + q1/2d0)*ds
         q2 = g(x0+k1/2d0, y0+l1/2d0, p0+m1/2d0)*ds
         k3 = dcos(p0 + m2/2d0)*ds
         l3 = dsin(p0 + m2/2d0)*ds
         m3 = (u0 + n2/2d0)*ds
         n3 = (v0 + q2/2d0)*ds
         q3 = g(x0+k2/2d0, y0+l2/2d0, p0+m2/2d0)*ds
         k4 = dcos(p0 + m3)*ds
         l4 = dsin(p0 + m3)*ds
         m4 = (u0 + n3)*ds
         n4 = (v0 + q3)*ds
         q4 = g(x0 + k3, y0 + l3, p0 + m3)*ds
         x1 = x0 + (k1 + 2d0*k2 + 2d0*k3 + k4)/6d0
         y1 = y0 + (l1 + 2d0*l2 + 2d0*l3 + l4)/6d0
         p1 = p0 + (m1 + 2d0*m2 + 2d0*m3 + m4)/6d0
         u1 = u0 + (n1 + 2d0*n2 + 2d0*n3 + n4)/6d0
         v1 = v0 + (q1 + 2d0*q2 + 2d0*q3 + q4)/6d0
         s = s + ds
         x0 = x1
         y0 = y1
         p0 = p1
         u0 = u1
         v0 = v1
         x(i) = x1
         y(i) = y1
   10 continue
      alpha = p1
      c1 = u1
      return
      end


c
c     g(x,y,p) = dv/ds
c
      function g(x,y,p)
      implicit double precision (a-h,o-z)
      g = -(y*dcos(p) - x*dsin(p))/4d0
      return
      end