Function with the symmetry of a hyperbolic regular tiling I'm looking for a meromorphic complex function

$$f: \mathbb{D} \to \mathbb{C}$$

where $\mathbb{D} = \{z\in \mathbb{C} : \vert z \vert \le 1\}$ is the unit disk identified with the Poincaré disk model of the hyperbolic plane, which has a zero of multiplicity 7 at each vertex of the order-7 triangular tiling of the disk:



(with one of the vertices at the origin and one of the central edges in $\mathbb{R}^+$, for example, so that the function is unique up to scale), and a pole of multiplicity 3 at each vertex of its dual tiling.

Has this function been studied somewhere? Is there an explicit expression for it, maybe in terms of special functions?

I know that similar functions exist for regular spherical tilings in the Riemann sphere (see for example this question), and powers of Weierstrass P functions give the corresponding functions for Euclidean tilings of the complex plane.
First of all, the proper description of your map is as a holomorphic map ${\mathbb D}\to CP^1= {\mathbb C}\cup \{\infty\}$ (this is what it means to have poles) which is automorphic with respect to a certain group $\Gamma$ of linear fractional transformations of ${\mathbb D}$ preserving the order order 7 tiling:
$$
f(\gamma z)= f(z), \gamma\in \Gamma, z\in {\mathbb D}. 
$$ 
Such $\Gamma$ is not unique and neither is $f$ (even after your normalization), unless you impose further restrictions. Firstly, you want to have $\Gamma$ which is the $(2,3,7)$ von Dyck group. To describe this group geometrically, pick one of the triangles $T$ in your tiling. Let $v$ be one of its vertices, $o$ the center of $T$ and $p$ the midpoint of an $e$ edge containing $v$. The group $\Gamma$ is discrete, generated by order 7 primitive elliptic rotation around $v$, order 2 elliptic rotation around $p$ and order $3$ rotation around $o$. (Rotations are understood here as linear fractional transformations conjugate to the euclidean rotations.) The quotient $X={\mathbb D}/\Gamma$ is topologically a 2-sphere with three distinguished points $[v], [o], [p]$, projections of $v, o, p$. The complex structure from  ${\mathbb D}$ descends to a complex structure on $X$ making it a Riemann surface. As any compact Riemann surface of genus $0$, it is biholomorphic to the extended complex plane, $CP^1$. Such a biholomorphic map $h$ is not unique, but it becomes unique if you require 
$$
h([o])=\infty, h([v])=0, h([p])=1. 
$$
The composition of the (holomorphic) projection $\pi: {\mathbb D}\to X$ with $h$ is your automorphic function $f$. It will have a triple pole at $o$, order 7 zero at $v$ and also will send the edge $e$ to  the unit segment $01$ in the real line. By automorphicity, the same will be true for all other vertices are centers of triangles. 

This description is not explicit of course. However, the function $f$ (or its multivalued inverse) can be found as a solution of some ODE of the form $S(F)=\varphi$, where $\varphi$ is a certain meromorphic quadratic differential and $S$ is the Schwarzian derivative. Equivalently, $F$ is the ratio of two independent solutions of Fuchsian linear 2nd order ODE, of the type: $y'' + \varphi y=0$, where $y$ is a function of one complex variable.  

This ODE can be made very explicit if you look at the target space (the extended complex plane); I can write it down if you are interested and when I have time. The differential $\varphi$ can be made very explicit on the target space but less so in the unit disk (one can find it as an infinite sum, a Poincare series). All of what I describe was known to Poincare and can be found in books on automorphic functions, like Ford's "Automorphic functions". 

Edit. The inverse to the function $f$ is the ratio $F_1/F_2$, where $F_1(z), F_2(z)$ are two linearly independent solutions of the following holomorphic 2nd order ODE
$$
z(z-1) F'' + ((a+b+1)z-c)F' + abF=0
$$
on the complex plane, where
$$
|1-c|=1/7, |c-a-b|=1/2, |a-b|=1/3.
$$
One has to further fiddle with the choice of $F_1, F_2$ to make sure that $f^{-1}$ sends the upper half plane $H=\{z: Im(z)&gt;0\}$ in the target complex plane to your favorite triangle $T$ (as above) in the unit disk.

See for instance Notes on differential equations and hypergeometric functions page 29. I do not have access to any math books right now and I got this reference just by googling "monodromy of hypergeometric functions". 

Once you know $f: T\to H$, you get $f: {\mathbb D}\to CP^1$ by using the Schwarz reflection principle (to define $f$ on the image of $T$ under the hyperbolic reflection in the edge $e$) and then the automorphic equation $f(\gamma z)= f(z)$ on the rest of the unit disk.
I have not heard of this.  

A lot of people who play with hyperbolic tilings check in with me to seek my opinion on this or that.

The ring structures grow very fast, this series https://oeis.org/A265041 shows the first few rings of {3,3,3,5} as found by my process.  

The equiangular or steriographic projection of hyperbolic geometry is bounded.  The only unbounded one i know is the equiperimeter version, which has rings that run as the cosh(nr).

But i can't ever recall anyone ever putting the h2 or e2 onto the complex plane.  A number of people have tried to suggest some sort of coordinate system, but i point out these turn very dense very quickly.