|||
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
20140103
20131101
坐标转换(ec2local.py)
20131030
客户端程序:
20131029
>>> import serial
>>> ser = serial.Serial(0) # open first serial port
>>> print ser.portstr # check which port was really used
>>> ser.write("hello") # write a string
>>> ser.close() # close port
从串口读东西(必须设置timeout,否则就锁死了。)>>> import serial
>>> ser =serial.Serial("COM3",9600,timeout=1)
>>> ser.isOpen()
>>>ser.write('1')
20130906
20130905
20130904
20130903
20130813
===============================================
===============================================
import numpy as np
#from mayavi.mlab import *
from mayavi import mlab
def test_flow():
x, y, z = np.mgrid[0:5, 0:5, 0:5]
r = np.sqrt(x**2 + y**2 + z**4)
u = y*np.sin(r)/r
v = -x*np.sin(r)/r
w = np.zeros_like(z)
obj = mlab.flow(u, v, w)
return obj
test_flow()
mlab.savefig(filename='D:Python27script\test_flow.jpg')
mlab.show()
===============================================
这个问题不知道怎么处理。不过已经不影响简单使用了。已经可以运行如下代码
============================================
==============================================
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-20 02:00
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社