Planet

SCREEN 13

' Earth landscape
FOR i = 1 TO 63: pal i, i, 31 + i / 2, i: NEXT

' Earth water
FOR i = 0 TO 7: pal 64 + i, i * 8, 16 + 6 * i, 63: NEXT

' Starfeld
FOR i = 0 TO 63: pal 72 + i, i, i, i: NEXT

TYPE vec3
  x AS SINGLE
  y AS SINGLE
  z AS SINGLE
END TYPE

RANDOMIZE 1

DIM c AS vec3, o AS vec3, sun AS vec3

' Planet center
o.x = 0
o.y = 0
o.z = 1.5

' Sun position
sun.x = 1
sun.y = 0
sun.z = -.5

' Details
dt = 32

normalize sun

FOR i = 0 TO 255: PSET (RND * 320, RND * 200), RND * 63 + 72: NEXT

DO

  FOR y = 100 TO -100 STEP -1
  FOR x = -160 TO 160

    ' Initial
    c.x = x / 100
    c.y = y / 100
    c.z = 1

    ' Intersection
    m = sphere(c, o, 1)

    ' Ray intersect
    IF m > 0 THEN

      r = RND

      ' Normal vector
      c.x = c.x * m - o.x
      c.y = c.y * m - o.y
      c.z = c.z * m - o.z
      normalize c

      ' UV-calc
      u = atan2(c.x, c.z)
      v = atan2(c.y, c.z)
      u = u + rot

      ' Get fractional part
      u = u - INT(u)
      v = v - INT(v)
      m = fbm(dt * u, dt * v) * 63

      ' Directional Light
      dl = c.x * sun.x + c.y * sun.y + c.z * sun.z
      IF dl < 0 THEN dl = 0

      ' Water or surface
      IF m < 32 THEN
        m = 64
      ELSE
        m = 2 * (m - 32)
      END IF

      IF r > dl THEN m = 72 + dl * 63
      PSET (160 + x, 100 - y), m

    END IF

  NEXT
  NEXT

  rot = rot + .025

LOOP WHILE INKEY$ = ""

FUNCTION atan2 (x, y)

  atan2 = 0
  pi = 3.141592
  IF x <> 0 OR y <> 0 THEN

    m = ATN(y / x)
    IF x >= 0 THEN m = m + pi
    atan2 = m / (2 * pi) + 1 / 4

  END IF

END FUNCTION

FUNCTION fbm (x, y)

  value = 0
  amp = .5
  freq = 0

  FOR i = 0 TO 5

    value = value + amp * noise(x, y)
    x = x * 2
    y = y * 2
    amp = amp * .5

  NEXT

  fbm = value

END FUNCTION

FUNCTION noise (x, y)

  ix = INT(x): fx = x - ix
  iy = INT(y): fy = y - iy

  a = rand(ix, iy)
  b = rand(ix + 1, iy)
  c = rand(ix, iy + 1)
  d = rand(ix + 1, iy + 1)

  ux = fx ^ 2 * (3 - 2 * fx)
  uy = fy ^ 2 * (3 - 2 * fy)

  noise = a * (1 - ux) + b * ux + (c - a) * uy * (1 - ux) + (d - b) * ux * uy

END FUNCTION

SUB normalize (c AS vec3)

  d = SQR(c.x ^ 2 + c.y ^ 2 + c.z ^ 2)
  c.x = c.x / d
  c.y = c.y / d
  c.z = c.z / d

END SUB

SUB pal (a, r, g, b)

  OUT 968, a
  OUT 969, r
  OUT 969, g
  OUT 969, b

END SUB

FUNCTION rand (x AS SINGLE, y AS SINGLE)

  m = SIN(x * 12.9898 + y * 78.233) * 43758.54531229988#
  rand = m - INT(m)

END FUNCTION

FUNCTION sphere (d AS vec3, o AS vec3, r AS SINGLE)

  sphere = -1

  a = d.x * d.x + d.y * d.y + d.z * d.z
  b = -2 * (d.x * o.x + d.y * o.y + d.z * o.z)
  c = o.x * o.x + o.y * o.y + o.z * o.z - r
  det = b ^ 2 - 4 * a * c

  IF det >= 0 THEN

    det = SQR(det)
    x1 = (-b - det) / (2 * a)
    x2 = (-b + det) / (2 * a)

    IF x1 < 0 AND x2 < 0 THEN sphere = -1
    IF x1 < 0 AND x2 > 0 THEN sphere = x2
    IF x1 > 0 AND x2 < 0 THEN sphere = x1
    IF x1 > 0 AND x2 > 0 AND x1 < x2 THEN sphere = x1
    IF x1 > 0 AND x2 > 0 AND x1 >= x2 THEN sphere = x2

  END IF

END FUNCTION