[daniel.corry]

I have written a code which writes the output data i desire to a .txt file on my desktop. I would no like to plot that data using pylab but can not work out how to extract the values and plot them all on the same graph.

I am attempting to plot dist on the x axis and ivrv, ivrg on the y axis. I have called the text file dist_GA_VA.txt.

Any help would be deeply appreciated.
Thanks

import numpy as n
import pylab
import math
import coord
import sys

def ezp(x):
    return math.exp(max(-20.0,x))

    if x==0 and y==0 and z==0:
        x=0.001
    
def sbf2flow(x,y,z):
    
    # Cosmology parameters
    h0 = 78.4
    wx = -55
    wy = 143
    wz =  -8
    omega = 0.2
    
    # Common parameters
    rcore = 2
    thermal = 187
 
    # VA parameters
    xv = -3.1
    yv = 16.6
    zv = -1.6
    deltav = 0.974
    gammav = 1.5
    rcutv = 11.7
    sigv = 650

    # GA parameters
    xg = -36.7
    yg =  14.6
    zg = -17.1
    deltag = 0.781
    gammag = 2.0
    rcutg = 49.5
    sigg = 500
    
    # Fornax parameters
    xf =  -1.9
    yf = -15.0
    zf = -13.4
    sigf = 235

    # Quadrupole parameters
    rquad = 50
    qxx =  2.57
    qxy =  2.04
    qxz = -7.87
    qyz = -3.63
    qzz =  8.55

    ######################################################

    # Peculiar velocity
    vx = wx
    vy = wy
    vz = wz
    ivrw = int((vx*x + vy*y + vz*z) / math.sqrt(x*x+y*y+z*z))
    

    # Hubble flow contribution
    vx = vx + h0*x
    vy = vy + h0*y
    vz = vz + h0*z
    ivrh = int(h0*math.sqrt(x*x+y*y+z*z))

    # Virgo Attractor
    distv = max(0.001, math.sqrt(xv**2+yv**2+zv**2))
    d3 = (distv/rcore)*(distv/rcore)*(distv/rcore)
    delta0 = deltav*d3/(ezp(-distv/rcutv)*((1+d3)**(1-gammav/3)-1))
    ratv = max(0.001, n.sqrt((x-xv)**2+(y-yv)**2+(z-zv)**2))
    d3 = (ratv/rcore)*(ratv/rcore)*(ratv/rcore)
    delta = delta0/d3*(ezp(-ratv/rcutv)*((1+d3)**(1-gammav/3)-1))

 
    uinfall = h0/3 * ratv * omega**0.6 * delta / (1+delta)**0.25
    ux = -uinfall * (x-xv)/ratv
    uy = -uinfall * (y-yv)/ratv
    uz = -uinfall * (z-zv)/ratv
    vx = vx + ux
    vy = vy + uy
    vz = vz + uz
    ivrv = int((ux*x + uy*y + uz*z) / n.sqrt(x*x+y*y+z*z))

    # Great Attractor
    distg = max(0.001, n.sqrt(xg**2+yg**2+zg**2))
    d3 = (distg/rcore)*(distg/rcore)*(distg/rcore)
    delta0 = deltag*d3/(ezp(-distg/rcutg)*((1+d3)**(1-gammag/3)-1))
    ratg = max(0.001, n.sqrt((x-xg)**2+(y-yg)**2+(z-zg)**2))
    d3 = (ratg/rcore)*(ratg/rcore)*(ratg/rcore)
    delta = delta0/d3*(ezp(-ratg/rcutg)*((1+d3)**(1-gammag/3)-1))
    uinfall = h0/3 * ratg * omega**0.6 * delta / (1+delta)**0.25
    ux = -uinfall * (x-xg)/ratg
    uy = -uinfall * (y-yg)/ratg
    uz = -uinfall * (z-zg)/ratg
    vx = vx + ux
    vy = vy + uy
    vz = vz + uz
    ivrg = int((ux*x + uy*y + uz*z) / n.sqrt(x*x+y*y+z*z))

    # Quadrupole
    ratq = n.sqrt(x*x+y*y+z*z)
    cut = ezp(-0.5*ratq*ratq/(rquad*rquad))
    qyy = -qxx - qzz
    ux = cut*(x*qxx + y*qxy + z*qxz)
    uy = cut*(x*qxy + y*qyy + z*qyz)
    uz = cut*(x*qxz + y*qyz + z*qzz)
    vx = vx + ux
    vy = vy + uy
    vz = vz + uz
    ivrq = int((ux*x + uy*y + uz*z) / n.sqrt(x*x+y*y+z*z))

    # Radial component

    vr = (vx*x + vy*y + vz*z) / n.sqrt(x*x+y*y+z*z)

    # Thermal velocity
    vsig = thermal*thermal
    vsig = vsig + sigv*sigv*ezp(-ratv*ratv/(rcore*rcore))
    vsig = vsig + sigg*sigg*ezp(-ratg*ratg/(rcore*rcore))
    # Fornax
    ratf = max(0.001, n.sqrt((x-xf)**2+(y-yf)**2+(z-zf)**2))
    vsig = vsig + sigf*sigf*ezp(-ratf*ratf/(rcore*rcore))
    vsig = n.sqrt(vsig)

    return (vx,vy,vz,vr,vsig,ivrq,ivrw,ivrg,ivrv,ivrh)

l_deg = input('galactic l:')
b_deg = input('galactic b:')

h0 = 78.4
l = (l_deg*math.pi)/180 # rad
b = (b_deg*math.pi)/180  # rad
(sgal_l, sgal_b ) = coord.gal2sgal(l, b)
sgx = math.cos(sgal_l)*math.cos(sgal_b)
sgy = math.sin(sgal_l)*math.cos(sgal_b)
sgz = math.sin(sgal_b)

out = open('dist_'+'GA_'+'VA'+'.txt','w')
for i in range(1,100,1):
    x = float(i)*sgx
    y = float(i)*sgy
    z = float(i)*sgz
    
    (vx,vy,vz,vr,vsig,ivrq,ivrw,ivrg,ivrv,ivrh) = sbf2flow(x,y,z)
    #print x*h0,y*h0,z*h0,vx,vy,vz
    
    dist= float(math.sqrt(x**2+y**2))

#    print dist
    out.write(str(dist) + ' '+str(ivrv)+' '+str(ivrg)+'\n')

out.close()


#pylab.plot(dist,ivrv)
#pylab.show()