# some quick python startup notes for linux #enter a linix bash shell with bash # enter the python environment python # to exit python use cntrl-D # do a little arithmetic 4/2 5/2 # this illustrates integer versus real arithmetic 5/2. # this illustrates integer versus real arithmetic 5*2 # multiplicatio 5**2 # exponents 6.99e11**4 # scientif notation and exponents x=2. y=4. z=x*y # assign variabels and do math print(z) # or just use z # to print the variable value print 'The value of z is', z # formatted print # try exponents 9**0.5 # try string variables name = 'Chip Kobulnicky' name # make a list of numbers x=[1,2,3,4,5] x x[0] # print first element x[1] # print second element x[1:4] # print a range of elements 1,2,3 len(x) # find the length of x l=len(x) # find the length of x and assign it to a variable int(3.54159) # take integer (does not round!) # to do serious math on arrays we need to import some # math modules import numpy as np x=np.array([1,2,3,4,5]) # define a true array of numbers x[1:3] # print elements 2 and 3 x[0:3] # print elements 1, 2 and 3 y=x*5 y z=np.arange(10,100,2) # make array from 10 to 100 every 2 units z np.log10(1000)# take a base 10 log np.sqrt(100) # take a square root np.sin(0) # sine function in radians!! np.sin(90) # sine function in radians!! np.exp(2) # e^2 exponentiation # make a plot # first import matplotlib import matplotlib.pyplot as plt x=np.array([1,2,3,4,5,6]) y=x**2 + 5*np.sqrt(x) # note use of np.sqrt for arrays! z=y*2 # note to do math on arrays you need to use np.sqrt(x) # or np.log10(x) # since the native python math. library does not handle arrays plt.plot(x,y,'bo',linestyle='-',label='My Plot') plt.plot(x,z,'r*',linestyle='--') plt.xlabel('X value') plt.ylabel('velocity (km/s)') plt.title('Some Title') plt.legend(loc='lower left') # adjust the axis range [xmin,xmax,ymin,ymax] plt.axis([0.1,5.5,0,100]) # save as a pdf file plt.savefig('test.pdf') # better yet run all of this as a script # put the above 10 lines in a file called # plot1.py and say # > import plot1 # at the python prompt # or at the shell level python plot1.py # put multiple plots on a page (x,y,#) plt.subplot(3,1,2) plt.subplots_adjust(hspace=0.1) # Here is how to make pretty formatted print statements area=3.4561e2 print ('Area is %5.2f' %(area) ) pc=2.9861e18 ly=3.2 print ('A parsec is %5.2e cm and about %i l.y.' %(pc, ly)) # or use the ;format method print ('A parsec is {0:5.2e} cm and about {1:3.1f} l.y.'.format(pc, ly)) # END OF BASIC STUFF # NOW FOR MORE ADVANCED STUFF # read in a multi coloumn data file from pylab import * mjd,mv,vel,v2=loadtxt('/d/www/chip/public_html/Classes/ASTR4610/Data/binary.dat',unpack=True,skiprows=6) # or do a formatted read n,name,mag=loadtxt('somefile.dat',unpack=True,skiprows=4,dtype='i,str,f') # for integer,string,float #fancier example of reading a csv files N,l,b,name,ra,dec,I1,I2,I3,I4,M24,W1,W2,W3,W4,H70,H160,Peak70,R,H,PeakJy=np.loadtxt('../../Bowshocks_master.csv', delimiter=',',unpack=True, skiprows=0, usecols=(0,1,2,3,4,5,6,7,8),dtype={'names' : ('N','l','b','name','ra','dec','I1'),'formats' : ('i','f','f','U19','U13','U13','f')} ) # and a formatted write ascii.write([N[ts],sname[ts],ra[ts],dec[ts],R[ts],H[ts],I1[ts]],'Tab1.dat',names=['N','sname','ra','dec','R','H','I1'],formats={'N':'%3i','sname':'s','ra':'s','dec':'s','R':'%2i','H':'%2i','I1':'%4f'}) # or data=np.array([a,b]) data=data.T np.savetxt('file.dat',data,fmt='%6.3f') # write a 2 or 3-column data file to values.dat from astropy.io import ascii x = np.array([1, 2, 3]) y = x ** 2 z = y + 2 ascii.write([x, y,z], 'values.dat', names=['x', 'y', 'z']) # or do a formatted write # fiest read data mjd,mv,vel,v2=loadtxt('/d/www/chip/public_html/Classes/ASTR4610/Data/binary.dat',unpack=True,skiprows=6) # then write data ascii.write([mjd, mv,vel], 'values.dat',names=['mjd', 'mv', 'vel'],formats={'MJD':'%12.2f','Mag(V)':'%4.2f','Vel(km/s)':'%5.2f'}) # plot a histogram # make a histogram # plot histogram import matplotlib.pyplot as plt plt.hist( vel, bins = 12 ) # or use the histogram method nv=np.histogram(v,bins=varray) np.shape(nv) # check the shape of array vals=nv[0] bincens=nv[1] plt.plot(bincens[0:-1],vals,'bo',fillstyle='full') # working with 2D arrays import numpy as np ; create a 3x3 array x=np.ndarray(shape=(3,3)) set second row to have value 4 x[1]=4 # set individual elements to specific values x[0,0]=1 x x[1,0]=2 x # make a 5x5 array x=np.ndarray(shape=(5,5)) y[2][3]=7 # extract a subarray t=y[1:3,2:4] # dealing with FITS files # TO READ: # import numpy as np from astropy.io import fits hdu = fits.open('myfile.fits') data=hdu[0].data hdr = hdul[0].header crval1=hdr['CRVAL1'] cdelt1=hdr['CDELT1'] crpix1=hdr['CRPIX1'] # to use header parameters to maek a wavelength array for a spectrum; CDELT1 specifies the increment per pixel and CRVAL1 gives the wavelength of the reference pixel, CRPIX1 which is usually 1, but not always. wavl=crval1 + cdelt1*(np.arange(0,len(data))-crpix1 ) TO WRITE the whole data structure into a new file: fits.writeto('out.fits', data, hdr) # # see http://docs.astropy.org/en/stable/io/fits/ # for more fits file options # # fit a Gaussian to a set of data # def gauss(x, a,b,c,d): # a, b, c, d = p y = a*np.exp(-np.power((x - b), 2.)/(2. * c**2.)) + d return y # set initial guesses for fit p_initial = np.array([0.02, -100, 30, 0.0]) # do the fit to x array v and observed y values ydata popt, pcov = curve_fit(gauss, v,ydata, p0=p_initial) # print best fit values a,b,c,d print(popt[0],popt[1],popt[2],popt[3] ) # print uncertainties on a,b,c print(np.sqrt(pcov[0,0]),np.sqrt(pcov[1,1 ]),np.sqrt(pcov[2,2 ]) ) ----- more on fits from time import time import sys import os import numpy as np from math import * from astrometry.util.starutil_numpy import * from astrometry.libkd import spherematch import random import pyfits from pyfits import * from numpy import radians, degrees, sin, cos, arctan2, hypot from astropy.io import fits if __name__ == '__main__': time0 = time() filename = 'test.fits' a = pyfits.open(filename) hdr = getheader(filename, 1) # header d = a[1].data ra = d.field('ra') dec = d.field('dec') z = d.field('z') col1 = fits.Column(name='RA(deg)', format='D', array=ra) col2 = fits.Column(name='DEC(deg)', format='D', array=dec) col3 = fits.Column(name='Redshift(zWARBNING0)', format='D', array=z) cols = fits.ColDefs([col1, col2, col3]) tbhdu = fits.BinTableHDU.from_columns(cols) outfile='test_out.fits' # clobber = True (it overwrites/updates) # tbhdu.writeto(outfile, clobber=True) # print 'Took ', time()-time0,' sec' # #