Hello and welcome to our community! Is this your first visit?
Register
Enjoy an ad free experience by logging in. Not a member yet? Register.
Results 1 to 1 of 1
  1. #1
    New to the CF scene
    Join Date
    Nov 2011
    Location
    durham
    Posts
    1
    Thanks
    0
    Thanked 0 Times in 0 Posts

    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()
    Last edited by Inigoesdr; 11-28-2011 at 08:38 PM.


 

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •