Classical modular forms of weight $k\geq 1$ can be computed explicitly by exhibiting their $q$-expansion. When the weight is at least two, the most standard method uses modular symbols. When the weight is one, we can not use modular symbols directly because such modular forms are not cohomological, but one can still compute their $q$-expansion by using ideas of Buhler, Frey, Serre, Buzzard and others. 

Say now that we want to compute explicitly modular forms of weight $k\geq 1$ on a Shimura curve associated to an order in an indefinite quaternion algebra over $\mathbb{Q}$. I can think of at least two ways of exhibiting such modular forms: either by (1) computing their Taylor expansion about a CM point or by (2) invoking Cerednik-Drinfeld theory and computing the harmonic cocycle on the Bruhat-Tits tree associated to it.

Method (1) is in principle available for all weights $k\geq 1$ but I'm not sure how well it works systematically in practice, as there is no canonical choice of CM point.

Method (2) seems more canonical and robust to me, but it has the drawback that it can only be used directly for $k\geq 2$. 

My question is: how would you compute in practice quaternionic modular forms of weight $1$? 

I seek an answer which should allow us to play around with these modular forms and do basic operations like for instance multiplying two forms of weight one and recognize the output as an explicit linear combination of eigenforms of weight two. This is beyond what one can gather from Jacquet-Langlands, I would say.