# Plotting from a .txt file using pylab

Printable View

• 11-28-2011, 05:37 PM
daniel.corry
Plotting from a .txt file using pylab
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

Code:

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