DECLARE SUB FUNCT (X!, Y!, Z!, T!)
DECLARE SUB VISIBPLOT (H(), L(), U1!, U2!, V1!, V2!, IC)
DECLARE SUB UPDATEHL (H(), L(), U1!, U2!, V1!, V2!)
DEF FNF (X, Y) = X * X * (1 - X) * Y ^ 3 * (1 - Y)
REM - GRAPHICS 3D
DIM H(640), L(640)
INPUT "Line of Sight   XE, YE, ZE : ", XE, YE, ZE
S1 = XE * XE + YE * YE
S2 = SQR(S1)
S3 = SQR(S1 + ZE * ZE)
S4 = 1 / (S2 * S3)
INPUT "# Grid Gradations : M, N : ", M, N
INPUT "# of Iter : ", NI
INPUT "Steps of : ", ST
DIM X(M), Y(N), R(M, N, 1), W(M, N, 1)
XL = 0
YL = 0
XH = 1
YH = 1
DX = (XH - XL) / M
DY = (YH - YL) / N
X0 = XH
IF XE < 0 THEN DX = -DX: X0 = XL
Y0 = YL
IF YE < 0 THEN DY = -DY: Y0 = YL
CX = 250
CY = 240
FOR I = 0 TO M
  X(I) = X0 - I * DX
NEXT
FOR J = 0 TO N
  Y(J) = Y0 + J * DY
NEXT
ZMAX = 1E-20
FOR I = 1 TO M - 1
  FOR J = 1 TO N - 1
    W(I, J, 0) = FNF(I / M, J / N)
    IF W(I, J, 0) > ZMAX THEN ZMAX = W(I, J, 0)
  NEXT
NEXT
FOR I = 1 TO M - 1
  FOR J = 1 TO N - 1
    W(I, J, 0) = 2.8 * W(I, J, 0) / ZMAX
  NEXT
NEXT

StartTime = TIMER
SCREEN 12
FOR L = 1 TO NI
MM = 1E+18
UH = -MM
UL = MM
VH = -MM
VL = MM
FOR I = 0 TO 639
  L(I) = 479
  H(I) = 0
NEXT

FOR I = 0 TO M
  FOR J = 0 TO N
    X = X(I)
    Y = Y(J)
    Z = W(I, J, 0)
    U = (XE * Y - YE * X) / S2
    V = (ZE * (X * XE + Y * YE) - S1 * Z) * S4
    R(I, J, 0) = U
    R(I, J, 1) = V
    IF U > UH THEN UH = U
    IF U < UL THEN UL = U
    IF V > VH THEN VH = V
    IF V < VL THEN VL = V
  NEXT
NEXT
REM (S IS SCALE FACTOR)
IF L = 1 THEN
  SA = 319 / ABS(UL)
  SB = 319 / ABS(UH)
  SC = 239 / ABS(VL)
  SD = 239 / ABS(VH)
  S = SA
  IF SB < S THEN S = SB
  IF SC < S THEN S = SC
  IF SD < S THEN S = SD
END IF
FOR I = 0 TO M
  FOR J = 0 TO N
    R(I, J, 0) = INT(S * R(I, J, 0) + CX)
    R(I, J, 1) = INT(S * R(I, J, 1) + CY)
  NEXT
NEXT
PRINT L, W(2 * M / 3, 2 * N / 3, 0)
  StartTime = TIMER
  TimePast = 0
  WHILE TimePast < 1
    TimePast = TIMER - StartTime
  WEND
  StartTime = TIME
  CLS
IF ABS(XE) >= ABS(YE) THEN
  FOR I = 0 TO M
    FOR J = 1 TO N
      U1 = R(I, J - 1, 0)
      V1 = R(I, J - 1, 1)
      U2 = R(I, J, 0)
      V2 = R(I, J, 1)
      IC = 8 + W(I, J, 0) * 7
      CALL VISIBPLOT(H(), L(), U1, U2, V1, V2, IC)
      CALL UPDATEHL(H(), L(), U1, U2, V1, V2)
    NEXT
    FOR J = 0 TO N
      IF I < M THEN
        U1 = R(I, J, 0)
        V1 = R(I, J, 1)
        U2 = R(I + 1, J, 0)
        V2 = R(I + 1, J, 1)
        IC = 8 + W(I, J, 0) * 7
        CALL VISIBPLOT(H(), L(), U1, U2, V1, V2, IC)
        CALL UPDATEHL(H(), L(), U1, U2, V1, V2)
      END IF
    NEXT
  NEXT
END IF
FOR LL = 1 TO ST
  FOR I = 1 TO M - 1
    FOR J = 1 TO N - 1
      W(I, J, 1) = (1 - .8) * W(I, J, 0) + .2 * (W(I - 1, J, 0) + W(I + 1, J, 0) + W(I, J - 1, 0) + W(I, J + 1, 0))
    NEXT
  NEXT
  FOR I = 1 TO M - 1
    FOR J = 1 TO N - 1
      W(I, J, 0) = W(I, J, 1)
    NEXT
  NEXT
NEXT
NEXT
DO
LOOP WHILE INKEY$ = ""

SUB FUNCT (X, Y, Z, T)
  IF X > 0 AND X < 3.141 AND Y > 0 AND Y < 3.141 THEN
    Z = X * (1 - X) * Y * (1 - Y) ^ 3
  END IF
END SUB

SUB UPDATEHL (H(), L(), U1, U2, V1, V2)
      IF V1 > H(U1) THEN H(U1) = V1
      IF V2 > H(U2) THEN H(U2) = V2
      IF V1 < L(U1) THEN L(U1) = V1
      IF V2 < L(U2) THEN L(U2) = V2
      IF ABS(U1 - U2) > 1 THEN
        KM = U1
        KX = U2
        VM = V1
        VX = V2
        IF KM > KX THEN
          KM = U2
          KX = U1
          VM = V2
          VX = V1
        END IF
        FOR K = KM + 1 TO KX - 1
          VK = VM + (VX - VM) * (K - KM) / (KX - KM)
          IF VK > H(K) THEN H(K) = VK
          IF VK < L(K) THEN L(K) = VK
        NEXT
      END IF

END SUB

SUB VISIBPLOT (H(), L(), U1, U2, V1, V2, IC)
      T1 = 0
      T2 = 0
      G1 = 0
      G2 = 0
      IF V1 >= H(U1) THEN T1 = 1
      IF V2 >= H(U2) THEN T2 = 1
      IF V1 <= L(U1) THEN G1 = 1
      IF V2 <= L(U2) THEN G2 = 1
      OK = 0
      IF T1 = 1 AND T2 = 1 THEN
        LINE (U1, V1)-(U2, V2), IC
        OK = 1
      END IF
      IF G1 = 1 AND G2 = 1 AND OK = 0 THEN
        LINE (U1, V1)-(U2, V2), IC
        OK = 1
      END IF
      IF T1 + T2 + G1 + G2 > 0 AND OK = 0 THEN
        KM = U1
        KX = U2
        VM = V1
        VX = V2
        IF KM > KX THEN
          KM = U2
          KX = U1
          VM = V2
          VX = V1
        END IF
        IF KM <> KX THEN
          F1 = 0
          F2 = 0
          FOR K = KM TO KX
            VK = VM + (VX - VM) * (K - KM) / (KX - KM)
            IF VK > H(K) OR VK < L(K) THEN
              U1 = K
              V1 = VK
              F1 = 1
              K = KX
            END IF
          NEXT
          FOR K = KX TO KM STEP -1
            VK = VM + (VX - VM) * (K - KM) / (KX - KM)
            IF VK > H(K) OR VK < L(K) THEN
              U2 = K
              V2 = VK
              F2 = 1
              K = KM
            END IF
          NEXT
          IF F1 = 1 AND F2 = 1 THEN LINE (U1, V1)-(U2, V2), IC
        ELSE
          IF VX > H(U1) THEN LINE (U1, H(U1))-(U1, VX), IC
          IF VM < L(U1) THEN LINE (U1, L(U1))-(U1, VM), IC
        END IF
      END IF

END SUB