' Source in Uniform Stream    Milne-Thompson p202     12/5/2008
     ' 2-dimensional case ...      CFD by Supercool
     ' Uniform stream flow over a ridge at Esperance
     ' Not to scale in both axes
     ' Code is Quick Basic
     DIM psi(30), x(30, 500), y(30, 500)      ' set up arrays
     SCREEN 9, , 1, 1: pi = 4 * ATN(1): rtd = 180 / pi: dtr = pi / 180
     ' rtd = radians to degrees: dtr = degrees to radians conversion
     scl = .6: WINDOW (-.8 * scl, -.1 * scl)-(.4 * scl, .5 * scl)
     LINE (-.6, .005)-(.3, 0), 3, BF
        U = 50        ' velocity of uniform stream     metres per seoond
        m = 4         ' source strength                cubic metres per radian
    ' find values of psi (flux) for streamlines at -infinity
    ' these streamlines are for the uniform flow with zero source strength
    ' they are chosen to give even vertical spacing of the initial streamlines
          dely = .02                ' spacing of initial streamlines in metres
        FOR n = 0 TO 30
          psi(n) = -U * n * dely    ' flux for streamline n
        NEXT n                      ' n * dely is height above water, metres
      ' plot free-stream streamlines
      ' delete ' below to see uniform flow streamlines
         ' COLOR 2: FOR n = 0 TO 20: LINE (-10, n * dely)-(10, n * dely): NEXT n
      ' plot source streamlines
         ' FOR theta = pi / 12 TO 11 * pi / 12 STEP pi / 12
         '      y = n * dely
         '      x = y / TAN(theta)
         '      COLOR 13: LINE (0, 0)-(x, y)
         ' NEXT theta
        FOR n = 0 TO 30                   ' each n is a new streamline
            k = 0: stepp = .01
      ' now run thru source angles to form streamline coordinates for each n
          FOR theta = pi - stepp TO pi / 1600 STEP -stepp   ' radians
      ' find y for psi source = m * theta
      ' do by brute force, trying all y's until we get close enough
      ' this works because psi (the flux) is constant in a streamline
              FOR y = .001 * dely TO 1000 STEP .001
                psi = -U * y - m * theta
                IF psi > psi(n) THEN GOTO 5 ELSE GOTO 10
5             NEXT y
10            k = k + 1
              x(n, k) = y / TAN(theta)    ' cordinates of streamline
              y(n, k) = y                 ' by trigonometry
          NEXT theta                      ' next location on chosen streamline
            endk = k - 1                  ' counter for points on streamline
        NEXT n                            ' next streamline
     ' plot streamlines on screen
           FOR n = 1 TO 30: FOR k = 2 TO endk
              COLOR 2:  LINE (x(n, k - 1), y(n, (k - 1)))-(x(n, k), y(n, k))
             NEXT k: NEXT n
      ' dividing streamline which separates source and uniform flows
      ' at y=0, speed of unifrom flow equals speed of source flow
      ' ie stagnation point at base of the cliff!
            n = 1: k = 0
          FOR theta = .001 TO pi STEP stepp
            y = -(m * theta - m * pi) / U     ' stagnation point at y=0
            x(n, k) = y / TAN(theta)
            y(n, k) = y
            k = k + 1
          NEXT theta
            endk = k - 1
      ' plot dividing streamline which separates source and uniform flows
             FOR k = 2 TO endk
              COLOR 14:  LINE (x(n, k - 1), y(n, (k - 1)))-(x(n, k), y(n, k))
             NEXT k
             FOR k = 1 TO endk
              COLOR 6: LINE (x(n, k), y(n, k))-(x(n, k + 1), 0), 6, BF
             NEXT k
           LOCATE 2, 1: COLOR 14:         PRINT " Symbol code:                 "
           LOCATE 3, 3: COLOR 3: PRINT "      Southern Ocean         "
           LOCATE 4, 3: COLOR 6: PRINT "      Ridge at Esperance     "
           LOCATE 5, 3: COLOR 2: PRINT "      Streamlines in the air "
           LOCATE 20, 47: PRINT " Airflow streamlines over a ridge "
           COLOR 0: INPUT q      ' hide prompt
           END