myprojects-hg-reborn

Fork

#! /usr/bin/env python
from scipy import *
import pylab

nt=129
nr=88
nz=129
read3d=0

if (read3d!=0):
   vel3d=pylab.load("vsection")
   print 'the shape of vel3d is', shape(vel3d)
   vel3d=vel3d.reshape(nt,nr,nz)
   vel_section=vel3d[:,:,65]
   pylab.save("vel_cross_section",vel_section)

elif (read3d==0):
   vel_section=pylab.load("vel_cross_section")

rg2=pylab.load("rg2.out")
r=rg2[:,0]


theta=linspace(0.,2.*pi,nt)
#print "theta is ", theta
sin_t=sin(theta)
cos_t=cos(theta)
rsin_t=r[newaxis,:]*sin_t[:,newaxis]
rcos_t=r[newaxis,:]*cos_t[:,newaxis]
rsin_t=ravel(rsin_t)
rcos_t=ravel(rcos_t)

rsin_t.shape=(nt,nr)
rcos_t.shape=(nt,nr)


#pylab.colorbar()
#pylab.clf()
pylab.figure()
X = rsin_t.transpose()
Y = rcos_t.transpose()
Z = vel_section.transpose()
velmin = vel_section.min()
velmax = vel_section.max()
print velmin, velmax
levels = arange(velmin, velmax+0.01, 0.01)
pylab.contourf(X, Y, Z, levels, cmap=pylab.cm.jet)
pylab.colorbar()

#pylab.show()

pylab.savefig("velocity_on_section_DNS.png")
#pylab.hold(False)