Преобразование значений плотности из сферической системы координат в декартову / скриптовую проверку

Question

У меня есть следующая проблема:

• сферическая система координат, куб данных 50*31*89, полный координат со значениями плотности и температуры
• Мне нужно перенести значение плотности из каждой сферической "ячейки" в соответствующую "ячейку" в декартовой системе с размером 150*150*50.

Я написал сценарий, но не уверен, что он правильный. Был бы рад, если бы кто-нибудь смог это проверить.

Вывод, который мне нужен, это файл ascii со значениями плотности. Позже я преобразую эти значения в двоичный файл данных вокселей, чтобы использовать его в Blender.

import numpy as np

# some "fix wired" information from input file
nI = 50
nJ = 31
nK = 89
n = nI * nJ * nK
radius = 3e10 

# define grid data arrays
r = np.zeros( nI )                 # radius
theta = np.zeros( nJ )             # polar angle
phi = np.zeros( nK  )              # azimuth angle
D = np.zeros( ( nI, nJ, nK ) )     # density
T = np.zeros( ( nI, nJ, nK ) )     # temperature 

# read data from input file
file = open( "PLANET3DP.DATA", "r" )
text = file.readlines()
file.close()
i = 0
for line in text:
    if( line.find( "Grid: radial" ) > -1 ):
        I1 = i + 2
    if( line.find( "Grid: polar" ) > -1 ):
        J1 = i + 2
    if( line.find( "Grid: azimuthal" ) > -1 ):
        K1 = i + 2
    if ( line.find( "Data:" ) > -1 ):
        data1 = i + 1
    i = i + 1
for I in range( 0, nI ):
    w = text[ I1 + I ].split()
    r[ I ] = float( w[ 1 ] )
#print( r )
for J in range( 0, nJ ):
    w = text[ J1 + J ].split()
    theta[ J ] = float( w[ 1 ] )
#print( theta )
for K in range( 0, nK ):
    w = text[ K1 + K ].split()
    phi[ K ] = float( w[ 1 ] )
#print( phi )
for i in range( 0, n ):
    w = text[ data1 + i ].split()
    I = int( w[ 0 ] ) - 1
    J = int( w[ 1 ] ) - 1
    K = int( w[ 2 ] ) - 1
    D[ I, J, K ] = float( w[ 3 ] )
    T[ I, J, K ] = float( w[ 4 ] )

# define 3 vectors that describe a 3D cartesian coordinate system
x = np.linspace(-radius,radius,150)
y = np.linspace(-radius,radius,150)
z = np.linspace(-radius,radius,50)

# pick densities from spherical coordiates
D2 = np.zeros( ( 150, 150, 50 ) )
for i in range( 0, 150 ):
    for j in range( 0, 150 ):
        for k in range( 0, 50 ):
            r2 = np.sqrt( x[ i ] ** 2 + y[ j ] ** 2 + z[ k ] ** 2 )
            theta2 = np.arccos( z[ k ] / r2 )
            phi2 = np.arctan2( y[ j ], x[ i ] )
            ii = np.where( r > r2 )
            jj = np.where( theta > theta2 )
            kk = np.where( phi > phi2 )
            if ii[ 0 ].size > 0 and jj[ 0 ].size > 0 and kk[ 0 ].size > 0:
                I = np.min( ii )
                J = np.min( jj )
                K = np.min( kk )
                D2[ i, j, k ] = D[ I, J, K ]
            else:
                D2[ i, j, k ] = 0
#print( D2 )
file = open( "only_density.csv", "w" )
for I in range( 0, 150 ):
    for J in range( 0, 150 ):
        for K in range( 0, 50 ):
            file.write( str( D2[ I, J, K ] ) + "\n" )
file.close()