20170517

使用starlink-pyast

import starlink.Grf as Grf
import starlink.Ast as Ast
import starlink.Atl as Atl
import matplotlib.pyplot as plt

#  Open the FITS file using (for instance) pyfits. A list of the HDUs
#  in the FITS file is returned.
hdu_list = pyfits.open( 't12_new.fits' )

#  Create a FitsChan and fill it with the FITS header cards in the
#  primary HDU (element zero of the list returned by pyfits.open).
fitschan = Ast.FitsChan( Atl.PyFITSAdapter(hdu_list[ 0 ]) )

#  Note which encoding has been used for the WCS information. This is only
#  necessary if you later wish to write modified headers out using the
#  same encoding. See section "...Write a Modified WCS Calibration to a
#  Dataset" below.
encoding = fitschan.Encoding

#  Read WCS information from the FitsChan. Additionally, this removes all
#  WCS information from the FitsChan.
wcsinfo = fitschan.read()
#print wcsinfo
#wcsinfo.show()

xpixel=[10,20]
ypixel=[20,40]
zpixel=[20,40]
print xpixel
print ypixel
print zpixel

(xworld,yworld,zworld)=wcsinfo.tran([xpixel,ypixel,zpixel])

print xworld
print yworld
print zworld

(xpixel2,ypixel2,zpixel2)=wcsinfo.tran([xworld,yworld,zworld],False)
print xpixel2
print ypixel2
print zpixel2

20170423

画天图（skymap）

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import pylab

lonsh,lonsm,latsd,latsm = np.loadtxt('FRB_radec.txt',unpack=True,usecols=[2,3,4,5])

lons=(lonsh+lonsm/60.0)*15.0
lats=(latsd+latsm/60.0)

m = Basemap(projection='hammer',lon_0=180)

parallels=np.arange(-90,90,30)
meridians=np.arange(0,360,60)

m.drawparallels(parallels,labels=[0,1,0,0],labelstyle="+/-",fmt='%g')
m.drawmeridians(meridians,labels=[1,0,0,1],labelstyle="+/-",fmt='%g')


rabins = np.arange(0.,360.,6)
decbins = np.arange(-90.,90.,6)
coverage = np.zeros([len(decbins),len(rabins)],"int")
for i in range(len(decbins)):
   for j in range(len(rabins)):
       coverage[i,j]=i+j
       
x,y= m(lons,lats)
m.scatter(x,y,3)
rabins, decbins = np.meshgrid(rabins, decbins)
x,y = m(rabins,decbins)
m.pcolor(x,y,coverage,cmap=pylab.cm.hot_r)


coordx,coordy=m(0,0)
plt.text(coordx,coordy,'360$^\circ$')
coordx,coordy=m(180,0)
plt.text(coordx,coordy,'180$^\circ$')

plt.show()

20170422

data3 = np.zeros(nrow,dtype=[('TSUBINT','float64'),('OFFS_SUB','float64'),('LST_SUB','float64'),('RA_SUB','float64'),('DEC_SUB','float64'),('GLON_SUB','float64'),('GLAT_SUB','float64'),('FD_ANG','float32'),('POS_ANG','float32'),('PAR_ANG','float32'),('TEL_AZ','float32'),('TEL_ZEN','float32'),('DAT_FREQ','float32',(nchan)),('DAT_WTS','float32',(nchan)),('DAT_OFFS','float32',(2*nchan)),('DAT_SCL','float32',(2*nchan)),('DATA','uint8',(nsblk,npol,nchan,1))])

注意，uint8表示8比特无符号整形变量，不能写成int8（有符号整形变量）。

20170401

pyephem中内置的主要行星和卫星

• ephem.Mercury

• ephem.Venus

• ephem.Mars, ephem.Phobos, ephem.Deimos

• ephem.Jupiter, ephem.Io, ephem.Europa, ephem.Ganymede,
  ephem.Callisto

• ephem.Saturn, ephem.Mimas, ephem.Enceladus, ephem.Tethys,
  ephem.Dione, ephem.Rhea, ephem.Titan, ephem.Iapetus, ephem.Hyperion,

• ephem.Uranus, ephem.Miranda, ephem.Ariel, ephem.Umbriel,
  ephem.Titania, ephem.Oberon

• ephem.Neptune

• ephem.Pluto

20161012

vc for python 2.7

http://origin.www.ms.akadns.NET/en-us/download/details.aspx?id=44266

matplotlib总是出现segment fault （core dump）问题，

pip uninstall matplotlib

pip install matplotlib

解决

20161011

easy_install ez-setup后成功安装FiPy

python解偏微分方程

http://ipython-books.github.io/featured-05/

安装了sympy和gravipy，可以进行广义相对论的符号计算。

无题

import numpy as np
import matplotlib.pyplot as plt
from pylab import *
from random import randint
#from scipy.optimize import curve_fit

def powerlaw(x,a,b,c):
   return a*x**b+c
def powerlaw2(x,a,c):
   return a*x**0.5+c

def Sll(xx,yy,zz,uu,vv,ww,nbin):
   ndim=len(xx)
   x=np.zeros([ndim,ndim])
   xt=np.zeros([ndim,ndim])
   y=np.zeros([ndim,ndim])
   yt=np.zeros([ndim,ndim])
   z=np.zeros([ndim,ndim])
   zt=np.zeros([ndim,ndim])
   u=np.zeros([ndim,ndim])
   ut=np.zeros([ndim,ndim])
   v=np.zeros([ndim,ndim])
   vt=np.zeros([ndim,ndim])
   w=np.zeros([ndim,ndim])
   wt=np.zeros([ndim,ndim])
   for i in xrange(ndim):
       x[i,:]=xx
       xt[:,i]=xx
       y[i,:]=yy
       yt[:,i]=yy
       z[i,:]=zz
       zt[:,i]=zz
       u[i,:]=uu
       ut[:,i]=uu
       v[i,:]=vv
       vt[:,i]=vv
       w[i,:]=ww
       wt[:,i]=ww
   dx=x-xt
   dy=y-yt
   dz=z-zt
   du=u-ut
   dv=v-vt
   dw=w-wt
   lengthtemp=np.sqrt(dx*dx+dy*dy+dz*dz)
   sflltemp=(du*dx+dv*dy+dw*dz)/lengthtemp
   length=lengthtemp[0,1:]
   sfll=sflltemp[0,1:]
   for i in xrange(1,ndim):
       length=np.append(length,lengthtemp[i,i+1:])
       sfll=np.append(sfll,sflltemp[i,i+1:])
   length2=np.arange(min(length),max(length),(max(length)-min(length))/nbin)
   sfll2=np.zeros(len(length2))
   for i in xrange(0,len(length2)):
       index=length>(length2[i]-(max(length)-min(length))/nbin/2)
       templ=length[index]
       temps=sfll[index]
       index2=templ<(length2[i]+(max(length)-min(length))/nbin/2)
       templ2=templ[index2]
       temps2=temps[index2]
       #sfll2[i]=np.sqrt(np.mean((temps2-np.mean(temps2))**2))
       sfll2[i]=np.mean((temps2-np.mean(temps2))**2)
   return np.array([length2,sfll2])

def CVD(xx,yy,uu,nbin):
   ndim=len(xx)
   x=np.zeros([ndim,ndim])
   xt=np.zeros([ndim,ndim])
   y=np.zeros([ndim,ndim])
   yt=np.zeros([ndim,ndim])
   u=np.zeros([ndim,ndim])
   ut=np.zeros([ndim,ndim])
   for i in xrange(ndim):
       x[i,:]=xx
       xt[:,i]=xx
       y[i,:]=yy
       yt[:,i]=yy
       u[i,:]=uu
       ut[:,i]=uu
   dx=x-xt
   dy=y-yt
   du=u-ut
   lengthtemp=np.sqrt(dx*dx+dy*dy)
   cvdtemp=du
   length=lengthtemp[0,1:]
   cvd=cvdtemp[0,1:]
   for i in xrange(1,ndim):
       length=np.append(length,lengthtemp[i,i+1:])
       cvd=np.append(cvd,cvdtemp[i,i+1:])
### </Note>: include both v_i-v_j and v_j-v_i
#        length=np.append(length,lengthtemp[i,:i-1])
#        cvd=np.append(cvd,cvdtemp[i,:i-1])
### <Note>
   length2=np.arange(min(length),max(length),(max(length)-min(length))/nbin)
   cvd2=np.zeros(len(length2))
   meandiff=np.zeros(len(length2))
   stt=np.zeros(len(length2))
   cvd=np.abs(cvd)
   for i in xrange(0,len(length2)):
       index=length>(length2[i]-(max(length)-min(length))/nbin/2)
       templ=length[index]
       temps=cvd[index]
       index2=templ<(length2[i]+(max(length)-min(length))/nbin/2)
       templ2=templ[index2]
       temps2=temps[index2]
       cvd2[i]=np.sqrt(np.mean((temps2-np.mean(temps2))**2))
   #    cvd2[i]=np.sqrt(np.mean((temps2)**2))
       meandiff[i]=np.mean(np.abs(temps2))
   #    meandiff[i]=np.median(np.abs(temps2))
       stt[i]=np.sqrt(np.mean((temps2)**2))
   #    stt[i]=np.mean((temps2)**2)
   length2=length2/57.3*140.0
   return np.array([length2,cvd2,meandiff,stt])

def gensample(zdimtemp,bet):
#    xdim=128
#    ydim=128
#    zdim=128
   xdim=8
   ydim=8
   zdim=8
   delta=(bet+2.0)/2.0
   v_fieldx=np.random.randn(xdim,ydim,zdim)
   v_fieldy=np.random.randn(xdim,ydim,zdim)
   v_fieldz=np.random.randn(xdim,ydim,zdim)

   vk_fieldx=np.fft.fftn(v_fieldx)
   vk_fieldy=np.fft.fftn(v_fieldy)
   vk_fieldz=np.fft.fftn(v_fieldz)
   for i in range(1,xdim+1):
       if(i>xdim/2):
         helpi=xdim-i
       else:
         helpi=i
       helpi=helpi*1.0
       for j in range(1,ydim+1):
           if(j>ydim/2):
             helpj=ydim-j
           else:
             helpj=j
           helpj=helpj*1.0
           for k in range(1,zdim+1):
               if(k>zdim/2):
                 helpk=zdim-k
               else:
                 helpk=k
               helpk=helpk*1.0
               knorm=np.sqrt(helpi**2+helpj**2+helpk**2)
               if(knorm>0):
                 vk_fieldx[i-1,j-1,k-1]=vk_fieldx[i-1,j-1,k-1]*knorm**(-delta)
                 vk_fieldy[i-1,j-1,k-1]=vk_fieldy[i-1,j-1,k-1]*knorm**(-delta)
                 vk_fieldz[i-1,j-1,k-1]=vk_fieldz[i-1,j-1,k-1]*knorm**(-delta)
   v_fieldpx=np.fft.ifftn(vk_fieldx).real
   v_fieldpx=v_fieldpx/np.sqrt(np.var(v_fieldpx))
   v_fieldpy=np.fft.ifftn(vk_fieldy).real
   v_fieldpy=v_fieldpy/np.sqrt(np.var(v_fieldpy))
   v_fieldpz=np.fft.ifftn(vk_fieldz).real
   v_fieldpz=v_fieldpz/np.sqrt(np.var(v_fieldpz))
   dim=xdim*ydim*zdim
   cx=np.zeros(dim)
   cy=np.zeros(dim)
   cz=np.zeros(dim)
   velox=np.zeros(dim)
   veloy=np.zeros(dim)
   veloz=np.zeros(dim)
   index=1
   for i in range(1,xdim):
       for j in range(1,ydim):
           for k in range(1,zdimtemp):
               if(np.random.rand()*xdim*ydim*zdimtemp < 4096):
                 cx[index-1]=i*1.0
                 cy[index-1]=j*1.0
                 cz[index-1]=k*1.0
                 velox[index-1]=v_fieldpx[i-1,j-1,k-1]
                 veloy[index-1]=v_fieldpy[i-1,j-1,k-1]
                 veloz[index-1]=v_fieldpz[i-1,j-1,k-1]
                 index=index+1
   cx=cx[0:index-1]
   cy=cy[0:index-1]
   cz=cz[0:index-1]
   velox=velox[0:index-1]
   veloy=veloy[0:index-1]
   veloz=veloz[0:index-1]
   return cx,cy,cz,velox,veloy,veloz
####
   


#xx,yy,zz,uu,vv,ww = np.loadtxt('taurus.txt',unpack=True,usecols=[0,1,2,3,6,9])
#xx,yy,zz,uu,vv,ww = np.loadtxt('twa.txt',unpack=True,usecols=[0,1,2,3,6,9])
#uu,xx,yy = np.loadtxt('core_taurus.txt',unpack=True,usecols=[6,16,17])
#uu,xx,yy = np.loadtxt('core_ophiuchus.txt',unpack=True,usecols=[6,16,17])
#uu,xx,yy = np.loadtxt('core_perseus.txt',unpack=True,usecols=[6,16,17])
#uu=uu/1000.0
zdimtemp=4
bet=1.0
xx,yy,zz,ux,uy,uz=gensample(zdimtemp,bet)


#result=CVD(xx,yy,ux,50)
result=Sll(xx,yy,zz,ux,uy,uz,50)
print result


lengthtemp=result[0,:]
upperscale=18.0
index=lengthtemp<=upperscale

#popt, pcov = curve_fit(powerlaw,result[0,index],result[1,index],[1,5,1])
#popt2, pcov2 = curve_fit(powerlaw,result[0,index],result[2,index],[1,5,1])
#popt3, pcov3 = curve_fit(powerlaw,result[0,index],result[3,index],[1,5,1])
#popt, pcov = curve_fit(powerlaw2,result[0,index],result[1,index],[1,1])
#popt2, pcov2 = curve_fit(powerlaw2,result[0,index],result[2,index],[1,1])
#popt3, pcov3 = curve_fit(powerlaw2,result[0,index],result[3,index],[1,1])
#print popt
#print popt3

plt.plot(result[0,:],result[1,:],marker='D',lw=0,mew=1.0,mec='#FFA500',ms=6.0,alpha=1.0,mfc='None')
#plt.plot(result[0,:],result[2,:],marker='p',lw=0,mew=1.0,mec='#FF0000',ms=6.0,alpha=1.0,mfc='None')
#plt.plot(result[0,:],result[3,:],marker='h',lw=0,mew=1.0,mec='#0000FF',ms=6.0,alpha=1.0,mfc='None')
xn=np.arange(0,upperscale,0.1)
#yn=powerlaw(xn,popt[0],popt[1],popt[2])
#yn2=powerlaw(xn,popt2[0],popt2[1],popt2[2])
#yn3=powerlaw(xn,popt3[0],popt3[1],popt3[2])
#yn=powerlaw2(xn,popt[0],popt[1])
#yn2=powerlaw2(xn,popt2[0],popt2[1])
#yn3=powerlaw2(xn,popt3[0],popt3[1])
#plt.plot(xn,yn,'k')
#plt.plot(xn,yn2,'k')
#plt.plot(xn,yn3,'k')
#plt.plot(rayso,decyso,marker='D',lw=0,mew=0.1,mec='#FFA500',ms=3.0,alpha=1.0,mfc='#FFA500')
xlabel('$L_{2D}$ (pc)')
#ylabel('$S_{ll}$(km/s)$^2$')
ylabel('CVD, CVD$_1$(km/s)')
#xlim(0,60)
#ylim(0,20)
plt.show()

20160115

读类psrfits文件

https://github.com/siemion/guppi_daq/blob/master/src/rawRead.py

20151226

easy_install astropy

总是报错，找不到vc9.0编译器。于是去Windows下python的民间扩展库（http://www.lfd.uci.edu/~gohlke/pythonlibs/）下载了一个astropy，解压缩以后把路径添加到python的路径（PYTHONPATH）里就可以了。

20151203

J2000->Jnow

20150821

mayavi只对python2.6管用。我到机器上python默认是2.7版的。2.6版也有，但是需要用
python2.6
启动。启动后发现没有pyfits，easy_install不管用，因为是针对python2.7的。怒了，直接用yum装了一个pyfits，世界清净了。

关于easy_install，2.6版的python可以用easy_install-2.6

20150423

python  读二进制表FITS文件

通常数据在第2页中（标号为1）

hdulist=pyfits.open(fitsname)

header=hdulist[1].header

20141225

画三色图，叠加外流、壳层结构的程序

import pyfits
from pylab import *

import math
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
#import matplotlib.colors as colormap
import mpl_toolkits
from matplotlib.patches import Ellipse
from matplotlib.patches import Wedge
from mpl_toolkits.axes_grid1.axes_rgb import make_rgb_axes, RGBAxes

hdulist24 = pyfits.open('/home/lqian/Work/Programming/idl/Observation/taurus24_reshape.fits')

image24 = hdulist24[0].data
nx24 = hdulist24[0].header['naxis1']
ny24 = hdulist24[0].header['naxis2']
crvalx24 = hdulist24[0].header['crval1']
cdeltax24 = hdulist24[0].header['cdelt1']
crpixx24 = hdulist24[0].header['crpix1']
crvaly24 = hdulist24[0].header['crval2']
cdeltay24 = hdulist24[0].header['cdelt2']
crpixy24 = hdulist24[0].header['crpix2']

x24 = np.arange(-crpixx24*cdeltax24+crvalx24,(nx24-1-crpixx24)*cdeltax24+crvalx24,cdeltax24)
y24 = np.arange(-crpixy24*cdeltay24+crvaly24,(ny24-1-crpixy24)*cdeltay24+crvaly24,cdeltay24)

name24='image24.npy'
name70='image70.npy'
name160='image160.npy'

image24=np.load(name24)
image70=np.load(name70)
image160=np.load(name160)

edge=400
image=np.zeros([ny24,nx24,3])
#image2=np.zeros([ny24,nx24,3])
image=np.zeros([ny24+2*edge,nx24+2*edge,3])
#image[:,:,0]=image160*0+1
#image[:,:,1]=image70*0+1
#image[:,:,2]=image24*0+1
#image[:,:,1]=image160*1.0/255
#image[:,:,0]=image70*1.0/255
#image[:,:,2]=image24*1.0/255

image[:,:,1]=image160
image[:,:,0]=image70
image[:,:,2]=image24


rad,ram,ras,decd,decm,decs,blueflag,blueangle,blueram,bluedecm,redflag,redangle,redram,reddecm = np.loadtxt('/home/lqian/Work/Programming/idl/Observation/bigmap_input.txt',unpack=True,usecols=[0,1,2,3,4,5,6,7,8,9,10,11,12,13])


ra=(rad+ram/60.0+ras/3600.0)*15.0
dec=decd+decm/60.0+decs/3600.0

convert=np.pi/180.0

cosc24=sin(crvaly24*convert)*sin(dec*convert)+cos(crvaly24*convert)*cos(dec*convert)*cos((ra-crvalx24)*convert)
xr24=crvalx24+cos(dec*convert)*sin((ra-crvalx24)*convert)/cosc24/convert
yd24=crvaly24+(cos(crvaly24*convert)*sin(dec*convert)- sin(crvaly24*convert)*cos(dec*convert)*cos((ra-crvalx24)*convert))/cosc24/convert

ra=xr24
dec=yd24


shellra,shelldec,shellr1,shellr2,shellangle1,shellangle2 = np.loadtxt('input_shell.txt',unpack=True,usecols=[0,1,2,3,4,5])


cosc24=sin(crvaly24*convert)*sin(shelldec*convert)+cos(crvaly24*convert)*cos(shelldec*convert)*cos((shellra-crvalx24)*convert)
xr24=crvalx24+cos(shelldec*convert)*sin((shellra-crvalx24)*convert)/cosc24/convert
yd24=crvaly24+(cos(crvaly24*convert)*sin(shelldec*convert)- sin(crvaly24*convert)*cos(shelldec*convert)*cos((shellra-crvalx24)*convert))/cosc24/convert

shellra=xr24
shelldec=yd24



plt.ion()
fig=plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)

minx=x24[nx24-1]+edge*1.0/nx24*(x24[nx24-1]-x24[0])
maxx=x24[0]-edge*1.0/nx24*(x24[nx24-1]-x24[0])
miny=y24[0]-edge*1.0/ny24*(y24[ny24-1]-y24[0])
maxy=y24[ny24-1]+edge*1.0/ny24*(y24[ny24-1]-y24[0])

#im=plt.imshow(image,extent=[min(x24),max(x24),min(y24),max(y24)],interpolation='none')
#plt.imshow(image,extent=[min(x24),max(x24),min(y24),max(y24)],interpolation='none')
plt.imshow(image,origin='lower',extent=[maxx,minx,miny,maxy],interpolation='none')
#plt.imshow(image,extent=[min(x24),max(x24),min(y24),max(y24)],interpolation='gaussian')
#plt.imshow(image,extent=[min(x24),max(x24),min(y24),max(y24)],interpolation='bicubic')
#ax=RGBAxes(fig,[0.1,0.1,0.8,0.8],pad=0.0)
#ax.imshow_rgb(image160,image70,image24)
xlabel('Ra')
ylabel('Dec')
#plt.colorbar(im,orientation='vertical')


for i in range(len(ra)):
   if(blueflag[i]==1):
      minor=np.min(blueram[i]/60.0,bluedecm[i]/60.0)
      bluelen=np.sqrt(blueram[i]**2+bluedecm[i]**2)/60.0
      blueopen=np.arctan(minor/2.0/bluelen)/convert
      blueangle[i]=90-blueangle[i]
      sector=Wedge((ra[i],dec[i]),bluelen,blueangle[i]-blueopen,blueangle[i]+blueopen,color='b',alpha=0.4,linewidth=0.0)
      ax.add_artist(sector)
   if(redflag[i]==1):
      minor=np.min(redram[i]/60.0,reddecm[i]/60.0)
      redlen=np.sqrt(redram[i]**2+reddecm[i]**2)/60.0
      redopen=np.arctan(minor/2.0/redlen)/convert
      redangle[i]=90-redangle[i]
      sector=Wedge((ra[i],dec[i]),redlen,redangle[i]-redopen,redangle[i]+redopen,color='r',alpha=0.4,linewidth=0.0)
      ax.add_artist(sector)

for i in range(len(shellra)):
#    sector=Wedge((shellra[i],shelldec[i]),(shellr2[i]+shellr1[i])/2.0/60.0,90-shellangle1[i],90.0-shellangle2[i],width=(shellr2[i]-shellr1[i])/60.0,color='w',alpha=0.4,linewidth=0.0)
   sector=Wedge((shellra[i],shelldec[i]),shellr2[i]/60.0,shellangle1[i],shellangle2[i],width=(shellr2[i]-shellr1[i])/60.0,color='w',alpha=0.4,linewidth=0.0)
   ax.add_artist(sector)


#sector=Wedge((65,20),5.0,0.0,90.0,width=0.2,color='w',alpha=0.4,linewidth=0.0)
#ax.add_artist(sector)

plt.show()
#savefig('taurus_outflow.eps')
savefig('taurus_outflow.pdf')
plt.close()

20141128
安装ppgplot需要用最新版本，1.0和1.1版都会报错。
http://code.google.com/p/ppgplot/downloads/detail?name=ppgplot-1.4.tar.gz&can=2&q=
http://ppgplot.googlecode.com/files/ppgplot-1.4.tar.gz

20140807

python读取二进制表fits文件

==============================================

#!/usr/bin/env python

import pyfits
import os,sys
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
import time

fitsname_array = os.popen('cd .../Data;ls .../Data/*.fit').readlines()
fitsname = fitsname_array[0].split()[0]


hdulist = pyfits.open(fitsname)
#print hdulist
imageheader = hdulist[0].header
#print imageheader
imagedata = hdulist[0].data
#print imagedata
tableheader = hdulist[1].header
print tableheader
tabledata = hdulist[1].data
#print tabledata[0]
for index in range(len(tabledata)):
   print index,'Source',tabledata[index][0],'RA',tabledata[index][1],'Dec',tabledata[index][2]

==============================================

20140718

==============================================

import os
import pyfits
import matplotlib.pyplot as plt
import time

#currentpath_array = os.popen('pwd').readlines()
#currentpath = currentpath_array[0].split()[0]
#print currentpath
datafiles_array = os.popen('cd /data;ls /data/*.fits').readlines()
#print datafiles_array
#commandstr='cd '
#commandstr+=currentpath
#print commandstr

plt.ion()
for i in range(len(datafiles_array)):
 print i
 datafiles = datafiles_array[i].split()[0]
 #print datafiles
 hdulist = pyfits.open(datafiles)
 hdu0 = hdulist[0]
 hdu1 = hdulist[1]
 data = hdu1.data
 spec = data.field(1)
 spec0 = spec[0,:]

 x = range(len(spec0))
 fig = plt.figure(figsize=(8,4))
 if (i%2 == 0):
   plt.semilogy(x,spec0,color='green')
 if (i%2 == 1):
   plt.semilogy(x,spec0,color='red')
 plt.draw()
 time.sleep(1)
 plt.close()

==============================================

20140628

看到一篇文章，讲用python实现访问google translate

http://www.phpfans.net/article/htmls/201104/MzM2Nzg2.html

20140430

Orrery

===================================

from socket import *
import time
import datetime
import struct,os,sys
import ephem
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np


getlocal = ephem.Observer()
getlocal.lon, getlocal.lat = '116.97345','40.55666' # at Miyun


# First set up the figure, the axis, and the plot element we want to animate  
fig = plt.figure()  
ax = plt.axes(xlim=(0, 360), ylim=(-90, 90),axisbg='#000000')  
line, = ax.plot([], [], 'yo')  
line2, = ax.plot([], [], 'r>')  
line3, = ax.plot([], [], color='green',ls='*',marker='.',markersize=3)  
line4, = ax.plot([], [], color='red',marker='*',markersize=20)  
line5, = ax.plot([], [], 'r-')  
line6, = ax.plot([], [], color='yellow',ls='*',marker='.',markersize=3)  

 
# initialization function: plot the background of each frame  
def init():  
   line.set_data([], [])  
   return line,
def azalt(ra,dec,ct3):
   getlocal.date = ct3 # UT
   body = ephem.FixedBody()
   body._ra = ra
   body._dec = dec
   body._epoch = ephem.J2000
   body.compute(getlocal)
   x= float(body.az)*180/np.pi
   y= float(body.alt)*180/np.pi
   return x,y

def azalt_of_mw(ct3):
   mwx = range(0,360,5)
   mwy = range(0,360,5)
   j=0
   for i in mwx:
       body_gl=i
       body_gb=0
       galactic_coord = ephem.Galactic(body_gl, body_gb)
       eq = ephem.Equatorial(galactic_coord)
 
       getlocal = ephem.Observer()
       getlocal.lon, getlocal.lat = '116.97345','40.55666' # at Miyun
       getlocal.date = ct3 # UT
 
       body = ephem.FixedBody()
       body._ra = eq.ra
       body._dec = eq.dec
       body._epoch = ephem.J2000
       body.compute(getlocal)
       mwx[j]=float(body.az)*180/np.pi
       mwy[j]=float(body.alt)*180/np.pi
       j=j+1
   return mwx,mwy

def azalt_of_eq(ct3):
   eqx = range(0,360,5)
   eqy = range(0,360,5)
   j=0
   for i in eqx:
       body = ephem.FixedBody()
       body._ra = i
       body._dec = 0.0
       body._epoch = ephem.J2000
       body.compute(getlocal)
       eqx[j]=float(body.az)*180/np.pi
       eqy[j]=float(body.alt)*180/np.pi
       j=j+1
   return eqx,eqy

def gc2azalt(gl,gb,ct3):
   body_gl=gl
   body_gb=gb
   galactic_coord = ephem.Galactic(body_gl, body_gb)
   eq = ephem.Equatorial(galactic_coord)
 
   getlocal = ephem.Observer()
   getlocal.lon, getlocal.lat = '116.97345','40.55666' # at Miyun
   getlocal.date = ct3 # UT
 
   body = ephem.FixedBody()
   body._ra = eq.ra
   body._dec = eq.dec
   body._epoch = ephem.J2000
   body.compute(getlocal)
   return float(body.az)*180/np.pi,float(body.alt)*180/np.pi
 
# animation function.  This is called sequentially  
# note: i is framenumber  
def animate(i):  
   x = np.random.rand(8)
   y = np.random.rand(8)
   ct = time.gmtime()
   ct2 = time.strftime("%Y/%m/%d",ct)
   ct3 = time.strftime("%Y/%m/%d %H:%M:%S",ct)

   sat = ephem.Saturn(ct2)
   jup = ephem.Jupiter(ct2)
   mar = ephem.Mars(ct2)
   ven = ephem.Venus(ct2)
   sun = ephem.Sun(ct3)
   ura = ephem.Uranus(ct3)
   nep = ephem.Neptune(ct3)
   mon = ephem.Moon(ct3)

   getlocal.date = ct3 # UT
   body = ephem.FixedBody()
   body._ra = sun.ra
   body._dec = sun.dec
   body._epoch = ephem.J2000
   body.compute(getlocal)
   x[0]= float(body.az)*180/np.pi
   y[0]= float(body.alt)*180/np.pi
#    x[0],y[0]=azalt(sun.ra,sun.dec,ct3)
   print ct3,sun.name, x[0],y[0]
   x[1],y[1]=azalt(sat.ra,sat.dec,ct3)
   print ct3,sat.name, x[1],y[1]
   x[2],y[2]=azalt(jup.ra,jup.dec,ct3)
   print ct3,jup.name, x[2],y[2]
   x[3],y[3]=azalt(mar.ra,mar.dec,ct3)
   print ct3,mar.name, x[3],y[3]
   x[4],y[4]=azalt(ven.ra,ven.dec,ct3)
   print ct3,ven.name, x[4],y[4]
   x[5],y[5]=azalt(ura.ra,ura.dec,ct3)
   print ct3,ura.name, x[5],y[5]
   x[6],y[6]=azalt(nep.ra,nep.dec,ct3)
   print ct3,nep.name, x[6],y[6]
   x[7],y[7]=azalt(mon.ra,mon.dec,ct3)
   print ct3,mon.name, x[7],y[7]
   
   line.set_data(x, y)  
   line2.set_data(x, y+10)  
   mwx,mwy = azalt_of_mw(ct3)
   line3.set_data(mwx,mwy)
   gcx,gcy=gc2azalt(0,0,ct3)
   line4.set_data(gcx,gcy)
   groundx=[0,360]
   groundy=[0,0]
   line5.set_data(groundx,groundy)
   eqx,eqy = azalt_of_eq(ct3)
   line6.set_data(eqx,eqy)
   time_text = ax.text(x[0],y[0]+5,'Sun',color='green')
   time_text2 = ax.text(x[1],y[1]+5,'Saturn',color='green')
   time_text3 = ax.text(x[2],y[2]+5,'Jupiter',color='green')
   time_text4 = ax.text(x[3],y[3]+5,'Mars',color='green')
   time_text5 = ax.text(x[4],y[4]+5,'Venus',color='green')
   time_text6 = ax.text(x[5],y[5]+5,'Uranus',color='green')
   time_text7 = ax.text(x[6],y[6]+5,'Neptune',color='green')
   time_text8 = ax.text(x[7],y[7]+5,'Moon',color='green')
   time_text9 = ax.text(gcx,gcy+5,'GC',color='green')
   return line,line2,line3,line4,line5,line6,time_text,time_text2,time_text3,time_text4,time_text5,time_text6,time_text7,time_text8,time_text9

 
# call the animator.  blit=True means only re-draw the parts that have changed.  
#anim = animation.FuncAnimation(fig, animate, init_func=init,  
#                               frames=200, interval=2000, blit=True)  
anim = animation.FuncAnimation(fig, animate, init_func=init,  
                               interval=2000, blit=True)  
 
#anim.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])  
 
plt.show()

20140417