Создание спиральной структуры в Python с использованием гиперболического тангенса

Question

Я пытаюсь создать спиральную структуру, как спиральные рукава галактики, в 2D-массиве на Python. Первый и простой способ я сделал это, использовал простую функцию логарифмической спирали, определяемую как на изображении: log spiral function

Значение x и y созданы с помощью

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

M и N являются размерами массива. Радиус координат прост, как и уравнение последнего изображения,

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

Создание яркости профиля из F (R), а ploting

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

дать мне Коммон спиральную структуру, как спираль галактика.

Другой способ сделать это - использовать другую функцию, hyperbolic tangent. В уравнениях последнего изображения, кроме r, которые определены, как и прежде, все остальные параметры, являются регулируемыми числами.

Для этой функции у меня есть проблемы с созданием спиральной структуры в 2D-массиве. Я не знаю, если мне нужно использовать гиперболический тангенс, чтобы сделать преобразование координат в массиве или искажение матрицы/массива, чтобы создать спиральную структуру. Я попробовал, но не мог.

Как я могу обработать эту spira/image, используя определения выше? Спасибо за помощь!

Более подробная информация о предмете, в ссылках:

1. Пэн, Ю. Чен и др; Детальное структурное разложение изображений галактики, 2002
2. Peng, Y. Chien et al; Подробное разложение изображений Галактики. II. За пределами осесимметричных моделей, 2009
3. Peng, Y. Chien, Galfit Руководство пользователя, 2003
4. Rowe, Barnaby et al; GALSIM: Модульная галактика изображений моделирования инструментарий, 2015

Отредактировано:

код, который я использую выглядит следующим образом:

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() 

Набор параметров может быть:

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 

Answer 1

Я не уверен, что это точно, но я думаю, что это близко, и дает результаты, похожие на бумагу:

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 

Единственное изменение состоит в том, что я использую

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

создать галактики сюжет.

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

создает этот сюжет:

И я спина его вокруг, добавляя k*Or

Используя это с n2 = 3 я получаю: