Creating a spiral structure in Python, using hyperbolic tangent

I'm trying to create a spiral structure, like spiral arms of a galaxy, in a 2D array in Python. The first and easy way I did it, was using a simple log-spiral function, defined as in the image:log spiral function

The x and y values are created by

x,y=meshgrid(arange(0,M=400,1), arange(0,N=400,1))

Mand N are the dimensions of the array. The radius coordinate is simple, like the equation of the last image,

r=(abs(x-gal_center[1])**(2.0)+((abs(y-gal_center[0]))/(q))**(2.0))**(0.5)

Creating the profile brightness of f(r), and ploting

plt.imshow((abs(galaxy_model))**0.2)

give me a commom spiral structure, like a spiral galaxy.

Another way to do this, is to use another function, the hyperbolic tangent. In the equations of the last image, unless r, that is defined like before, all the others parameters, are ajustable numbers.

For this function, I have problems to make a spiral structure in a 2D array. I don't know, if I need to use the hyperbolic tangent to make a coordinate transformation in the array, or a matrix/array distortion, to create a spiral structure. I tried it, but I could not.

How can I proced to make this spira/image, using the definitions above? Thanks for the help!

More information about the subject, in the references:

1. Peng, Y. Chien et al; Detailed Structural Decomposition of Galaxy Images, 2002
2. Peng, Y. Chien et al; Detailed Decomposition of Galaxy Images. II. Beyond Axisymmetric Models, 2009
3. Peng, Y. Chien, Galfit User's Manual, 2003
4. Rowe, Barnaby et al; GALSIM:The modular galaxy image simulation toolkit, 2015

Edited:

The code that I'm using is as follows:

from __future__ import division
import numpy as np
from numpy import*
import matplotlib.pyplot as pyplot
import scipy as sp
from scipy import*
import pylab as pl
from pylab import*
import math 
from math import*
import pyfits as pf
from pyfits import*

def exponential_profile(Io,ro,r):
    Iexp=0.5*Io*np.exp(-r/ro)
    return Iexp

def sersic_profile(Io,ro,r,n):
    Iser=Io*np.exp(-(r/ro)**(1/n))
    return Iser

def galaxy_model1(q,c,gal_center,Io,ro,n,M,N,xi,p,n1,n2,s1,s2,k):
    x,y=meshgrid(arange(-M/2,M/2,1), arange(-N/2,N/2,1))
    r=(abs(x-0*gal_center[1])**(c+2.0)+((abs(y-0*gal_center[0]))/(q))**(c+2.0))**(1.0/(c+2.0))
    power=2.0
    fr=(30-xi*np.log(1.0+r**power)+(1.0/p)*np.cos(n1*arctan2(x,y)+k*np.log(s1+r**power))+(1.0/p)*np.cos(n2*arctan2(x,y)+k*np.log(s2+r**power))  )
    I_exp=exponential_profile(Io,ro,r)
    I_ser=sersic_profile(Io,ro,r,n)
    galaxy_model_1=0.1*I_exp+0.1*I_ser+0.5*fr
    return galaxy_model_1

def galaxy_model2(q,c,Cb,rout,rin,Oout,a,M,N,Io,ro,n):
    gal_center=(M/2,N/2)
    x,y=meshgrid(arange(0,M,1), arange(0,N,1))
    r=(abs(x-0*gal_center[1])**(c+2.0)+((abs(y-0*gal_center[0]))/(q))**(c+2.0))**(1.0/(c+2.0))
    A=2*Cb/(abs(Oout)+Cb)-1.00001
    B=(2-np.arctanh(A))*((rout)/(rout-rin))
    T=0.5*(np.tanh(B*(r/rout-1)+2)+1)
    Or=Oout*T*(0.5*(r/rout+1))**a
    I_exp=exponential_profile(Io,ro,r)
    I_ser=sersic_profile(Io,ro,r,n)
    galaxy_model_2=0.1*I_exp+0.1*I_ser+0.5*Or
    return galaxy_model_2
galaxy_model_1=galaxy_model1(q,c,(M/2,N/2),Io,ro,n,M,N,xi,p,n1,n2,s1,s2,k)
galaxy_model_2=galaxy_model2(q,c,Cb,rout,rin,Oout,a,M,N,Io,ro,n)
fig=plt.figure()
ax1=fig.add_subplot(121)
ax1.imshow((abs(galaxy_model_1))**0.2)
pf.writeto('gal_1.fits', galaxy_model_1,  clobber=1)  
ax2=fig.add_subplot(122, axisbg='white')
ax2.imshow((abs(galaxy_model_2))**0.2)
plt.show()

A set of parameters can be:

M=400
N=400
q=0.8
c=0.0
Io=100.0
ro=10.0
n=3.0
xi=2.0
p=1.7
n1=3.0
n2=3.0
s1=0.05
s2=0.5
k=3.0
Cb=0.23
rout=100.0
rin=10.0
Oout=pi/2
a=0.0

I'm not sure this is exactly right but I think it is close, and produces results similar to the paper:

def galaxy_model2(q,c,Cb,rout,rin,Oout,a,M,N,Io,ro,n):
    gal_center=(0,0)
    x,y=meshgrid(arange(-M/2,M/2,1), arange(-N/2,N/2,1))
    r=(abs(x-gal_center[1])**(c+2.0)+((abs(y-gal_center[0]))/(q))**(c+2.0))**(1.0/(c+2.0))
    A=2*Cb/(abs(Oout)+Cb)-1.00001
    B=(2-np.arctanh(A))*((rout)/(rout-rin))
    T=0.5*(np.tanh(B*(r/rout-1)+2)+1)
    Or=Oout*T*(0.5*(r/rout+1))**a
    Or=30-np.log(1.0+r**2.0)+(2.0/p)*np.cos(n2*arctan2(x,y)+k*Or)
    I_exp=exponential_profile(Io,ro,r)
    I_ser=sersic_profile(Io,ro,r,n)
    #galaxy_model_2=0.5*Or
    return Or

The only change is that I use

Or=30-np.log(1.0+r**2.0)+(2.0/p)*np.cos(n2*arctan2(x,y)+k*Or)

to create a galaxy plot.

np.cos(n1*arctan2(x,y))

creates this plot:

And i spin it around by adding k*Or

Using this with n2=3 I get: