calculate :: Mask -> XY -> Int -> Mask
calculate msk@(A2D.Array2d cols rows _) c r =
    let size = (cols, rows) in
    let dirs = [ (i, j) | i <- [-1..1], j <- [-1..1], i /= 0 || j /= 0 ] in
    foldl (\msk (i, j) -> let l1 = castLight c size 1 1.0 0.0 0 i j 0 False msk in
                          let l2 = castLight c size 1 1.0 0.0 i 0 0 j False l1 in
                          l2
                      ) msk dirs

 where castLight :: XY -> (Int,Int) -> Int -> Double -> Double -> Int -> Int -> Int -> Int -> Bool -> Mask -> Mask
       castLight c@(sx, sy) size@(width,height) row start end xx xy yx yy blocked mask =
           if start < end
           then mask
           else loop1 row start 0.0 blocked mask
        where loop1 :: Int -> Double -> Double -> Bool -> Mask -> Mask
              loop1 distance start newStart blocked mask =
                  let deltaY = -distance in
                  if distance > r || blocked
                  then mask
                  else loop2 (-distance) (-distance) start newStart blocked mask
               where loop2 :: Int -> Int -> Double -> Double -> Bool -> Mask -> Mask
                     loop2 deltaY deltaX start newStart blocked mask =
                         if deltaX > 0
                         then mask
                         else let currentX = sx + deltaX * xx + deltaY * xy in
                              let currentY = sy + deltaX * yx + deltaY * yy in
                              let leftSlope = fromIntegral deltaX - 0.5 / fromIntegral deltaY + 0.5 in
                              let rightSlope = fromIntegral deltaX + 0.5 / fromIntegral deltaY - 0.5 in
                              if not (currentX >= 0
                                   && currentY >= 0
                                   && currentX < width
                                   && currentY < height)
                                || start < rightSlope
                              then loop2 deltaY (deltaX + 1) start newStart blocked mask
                              else if end > leftSlope
                              then mask
                              else let mask' = if inRadius (deltaX, deltaY) r
                                               then A2D.put mask (currentX, currentY) True
                                               else mask
                                   in
                                   if blocked
                                   then if any not (A2D.get mask' (currentX, currentY) )
                                        then loop2 deltaY (deltaX + 1) start rightSlope blocked mask'
                                        else loop2 deltaY (deltaX + 1) newStart newStart False mask'
                                   else if any not (A2D.get mask' (currentX, currentY))
                                        then let mask'' = castLight c size (distance + 1) start leftSlope xx xy yx yy False mask' in
                                             loop2 deltaY (deltaX + 1) start rightSlope True mask''
                                        else loop2 deltaY (deltaX + 1) start newStart blocked mask'