Lensing

'Lensing by @_sorceress in April 2020
'Run in QB64
'-------------------------------------------------------

CurvatureLeft = 10
CurvatureRight = 10
LensThickness = 2

SCREEN 13
PRINT "use numeric pad to change lens"
SLEEP

DO
  k$ = INKEY$
  IF k$ = "1" THEN CurvatureLeft = CurvatureLeft + .5
  IF k$ = "4" THEN CurvatureLeft = CurvatureLeft - .5
  IF k$ = "2" THEN LensThickness = LensThickness + .5
  IF k$ = "5" THEN LensThickness = LensThickness - .5: IF LensThickness < 0 THEN LensThickness = 0
  IF k$ = "3" THEN CurvatureRight = CurvatureRight + .5
  IF k$ = "6" THEN CurvatureRight = CurvatureRight - .5

  'clear
  CLS

  'draw lens
  LINE (0, 100)-(320, 100), 8
  FOR t = -1 TO 1 STEP .01
    GOSUB getLens
    LINE (xL, y)-(xR, y), 11
  NEXT

  'draw ray
  sy = sy + 1: IF sy >= 200 THEN sy = sy - 200

  FOR col = 40 TO 55 STEP .5 'rainbow
    ok = 0: oldinLens = 0
    rayX = 0: rayY = sy: rayA = 0
    PSET (rayX, rayY)
    eta = 1.5 + .08 * (col - 40) / 15 '1.50 to 1.58 for glass
    FOR i = 0 TO 400
      oldinLens = inLens
      rayX = rayX + COS(rayA)
      rayY = rayY + SIN(rayA)
      t = (rayY - 100) / 50: GOSUB getLens
      IF rayX > xL AND rayX < xR AND ABS(t) < 1 THEN inLens = 1 ELSE inLens = 0
      dLens = inLens - oldinLens
      IF dLens <> 0 THEN 'lens transition
        IF dLens > 0 THEN e = 1 / eta 'air to glass
        IF dLens < 0 THEN e = eta 'glass to air
        ok = 1
        t = (rayY - 100) / 50
        GOSUB getLens
        IF ABS(rayX - xL) < ABS(rayX - xR) THEN nA = ATN2(nyL, nxL) ELSE nA = ATN2(nyR, nxR)
        IF dLens = -1 THEN nA = nA + _PI
        g = SIN(rayA - nA) * e
        IF ABS(g) < 1 THEN
          g = g / SQR(1 - g * g): rayA = ATN(g) + nA 'refraction!
        ELSE
          rayA = 2 * nA - rayA 'reflection!
        END IF
      END IF
      IF ok = 1 THEN c = col ELSE c = 15
      LINE -(rayX, rayY), c
    NEXT
  NEXT

  _DISPLAY
  _LIMIT 60

LOOP
END

getLens:
y = 100 + 50 * t
xL = 80 - LensThickness - (1 - t * t) * CurvatureLeft
xR = 80 + LensThickness + (1 - t * t) * CurvatureRight
nxL = 50: nyL = -2 * t * CurvatureLeft
nxR = -50: nyR = -2 * t * CurvatureRight
RETURN

FUNCTION ATN2 (dy, dx)
  IF dx <> 0 THEN
    ATN2 = ATN(dy / dx)
    IF dx < 0 THEN
      IF dy >= 0 THEN ATN2 = ATN2 + _PI ELSE ATN2 = ATN2 - _PI
    END IF
  ELSE
    ATN2 = SGN(dy) * _PI / 2
  END IF
END FUNCTION