a better approach than tracing path numerically could be to analytically calculate Fourier series of inverse Boettcher map (Laurent series) applied to w = r exp(i theta) and take the limit as r -> 1 from above

https://en.wikipedia.org/wiki/External_ray#Uniformization_2**edit** I think it works out quite easily:

\[

\mathcal{F}_t(\Psi_{M}(r e^{i t}))(\omega) =r \delta(\omega + 1) - {\frac{1}{2}}\delta(\omega)+{\frac{1}{8r}}\delta(\omega-1)-{\frac{1}{4r^{2}}}\delta(\omega-2)+{\frac{15}{128r^{3}}}\delta(\omega - 3)+\ldots

\]

**edit2** see also

https://www.mrob.com/pub/muency/laurentseries.htmlJust tried it, converges rather slowly.

Haskell code using the data-memocombinators package:

`import qualified Data.MemoCombinators as Memo`

import Control.Monad (forM_)

import Data.Complex (Complex((:+)), cis)

import Data.List (intercalate)

-- https://www.mrob.com/pub/muency/laurentseries.html#coef

b :: Integer -> Rational

b = Memo.integral b'

where

b' 0 = -1/2

b' n = -(w n (n + 1)) / fromInteger n

w :: Integer -> Integer -> Rational

w = Memo.memo2 Memo.integral Memo.integral w'

where

w' 0 _ = 0

w' n m = a (m - 1) + w (n - 1) m + sum [ a j * w (n - 1) (m - j - 1) | j <- [0 .. m - 2] ]

a :: Integer -> Rational

a = Memo.integral a'

where

a' j = u 0 (j + 1)

u :: Integer -> Integer -> Rational

u = Memo.memo2 Memo.integral Memo.integral u'

where

u' n k

| 2 ^ n - 1 == k = 1

| 2 ^ n - 1 > k = sum [ u (n - 1) j * u (n - 1) (k - j) | j <- [ 0 .. k ] ]

| 2 ^ (n + 1) - 1 > k = 0

| otherwise = (u (n + 1) k - sum [ u n j * u n (k - j) | j <- [ 1 .. k - 1 ] ]) / 2

psi :: Double -> [Double]

psi r = r : [ fromRational (b n) / r^n | n <- [(0 :: Integer) .. ] ]

point :: Int -> Double -> Double -> Complex Double

point n r t = sum [ fmap (p *) $ cis t ^^ m | (p, m) <- take n $ psi r `zip` [ 1, 0 .. ] ]

curve :: Int -> Int -> Double -> [Complex Double]

curve m n r = [ point n r (2 * pi * t) | i <- [ 0 .. m ], let t = fromIntegral i / fromIntegral m ]

plot :: [Complex Double] -> String

plot ps = unlines (map plot1 ps)

where

plot1 (x :+ y) = unwords (map show [x, y])

main :: IO ()

main = do

forM_ [ 0 .. 9 ] $ \l -> do

forM_ [ 1 .. 100 ] $ \n -> do

let r = 1 + 0.5 ^ l

writeFile (show l ++ "_" ++ show3 n ++ ".dat") (plot $ curve 1024 n r)

show3 :: Int -> String

show3 = reverse . take 3 . (++ "000") . reverse . show

**edit** there was a bug, code is fixed now...

image is with r = 1 + 0.5^9, 100 harmonics