Transcript 计算天文讲义

16. IDL速成及其它可视化工具简介
P. F. Chen
Department of Astronomy
Nanjing University
July 2, 2000
IDL (Interactive Data Language)



A complete computing environment for
the interactive analysis and visualization
of data.
It integrates a powerful language with
numerous mathematical analysis and
graphical display techniques.
It is a time-saving alternative to
programming in Fortran or C.
IDL特点

可交互式使用，亦可生成复杂的函数或程序；

运算可以数组为单位：如x(100), y(100)满足 y(i)=sin(x(i)) ，在IDL中
可表示为y=sin(x)；

直接编译、运行，及时反馈；

快速的多维图象显示及动画，易于诊断结果；

很多计算、统计程序（含Numerical Recipes）；

与Fortran、C可相互调用。

重要的是它提供了一个控制、处理及绘图的软件发展平台

可生成独立的可执行文件
运行方式

通常在IDL环境下运行，也可单独运行
(点击IDL图标for Windows；点击图标或运行idl or idlde for Unix or Linux)

1. 交互式
.r
2. IDL程序
***.pro
可直接输入，也可先敲

pro exam1
;上一行可有可无
device,decomposed=0
!p.color=0
!p.background=255
x=findgen(51)/50.*3.14 & y=sin(x)
plot, x, y,xst=1,thick=5 & y=y/1.2
oplot, x, y,psym=1
& y=y/1.2
oplot, x, y,psym=2
& y=y/1.2
oplot, x, y,linestyle=2
& y=y/1.2
window,1,xs=600,ys=500
plot, x, y,psym=10,xr=[1.,2.5],yr=[0.3,0.5]
end
常用的运算及语句







+ - * / ^ (exponentiation) sqrt ##(矩阵相乘)，
invert()矩阵求逆
and not or xor
pro exam2
for i=0,10 do begin
eq ge gt le lt ne
if i ne 7 then goto,l2
循环
print,'真的很好用'
判断
l2: print,'Now i=',$
i
GOTO
endfor
续行符 $
end
赋值





a=1. & b=2.d & c=’ ’ （c为字符型变量）
a=fltarr(200,/nozero) & b=bytarr(512,512)
d=[1.,2.,3.]
e=[[1.,2.,3.],[4.,5.,6.],[7.,8.,9.]] 矩阵
read data from a file
寻找文件
……
fn=dialog_pickfile(/read, /multiple_files,filter = ['*.jpg', '*.tif'])
if fn eq '' then return
openr,1,fn
……
寻找几个文件，若1个则不用multiple_files
有格式文件



打开
存取
关闭
a=fltarr(100)
a=fltarr(2,3,2)
openr,2,’input’
openr,2,’D:\input.txt’
readf,2,a
readf,2,a
close,2
close,2
openw,3,’output’
printf,a,b
close,3
Fortran程序中无格式保存数据时
write(1)((a(i,j),i=1,200),j=1,200)
无格式文件

打开
实型
a=fltarr(200,200,/nozero)
openr,2,’data’,/f77_unformatted
readu,2,a

读取
close,2
字节数据

关闭
a=bytarr(512,512)
openr,2,’data’
readu,2,a
close,2
Fits格式文件读取
fitsdata,’filename’,data
其中fitsdata.pro为一个子程序
data为2维数组，通过命令help,data可以知道其维数。
或fits_read,’filename’,data,index
画图

device,decomposed=0
loadct,3
a=dist(500)
tvscl,a
Plot, oplot, errplot, contour, tvscl, surface等及各
自的交互式命令，如live_surface
均匀分布的二维数组

屏幕显示，可多图

产生图象文件
多图显示可采用两种方式：
1. 设定!p.multi=[0,2,2]
2. 在命令后加pos=[0.1,0.1,0.3,0.3]等
不均匀分布的二维数组
a=dist(500)
x=findgen(500) & y=findgen(500)^2/499.
contour,a,x,y,nlevels=51,/fill
contour,a,x,y,nlevels=11
图像诊断及去坏点






cursor,x,y
box_cursor,x,y,dx,dy
profiles,a
b=profile(a)
a(where(a eq 999))=0
a(where(strtrim(a,2) eq strtrim(!values.f_nan,2)))=0.
求解代数方程
pro exam3
x= [1.0, 5.0]
;Provide an initial guess as the algorithm's starting point.
result = newton(x, 'newtfunc')
;Compute the solution.
print, result
end
x1  x2  3  0
function newtfunc, x
return, [x[0] + x[1] -3.0, 2.*x[0] + x[1] - 4.0]
; return, [x[0] + x[1] -3.0, x[0]^2 + x[1]^2 - 9.0]
end
2 x1  x2  4  0
x1  x2  3  0
x12  x 22  9  0
拟合曲线
y  0.9  ae
bx
pro exam4
x=findgen(10)
;;; X values of the sample points
y=exp(-x)+randomU(seed,10)*0.01+1. ;;; Y values of the sample points
weights=1./y
a=[2.,-2.]
;;;; initial guess of the coefficients
yfit=curvefit(x,y,weights,a,sigma,function_name='gfunct')
print,a
plot,x,y,psym=6,thick=5
oplot,x,0.9+a(0)*exp(a(1)*x)
end
pro gfunct,x,a,f,pder
bx=exp(a[1]*x)
f=a[0]*bx+0.9
;;;;; 0.9 can be changed
if n_params() ge 4 then pder=[[bx],[a[0]*x*bx]]
;;; pder is the expression of the partial differentiation with respect to
;;;
the fitting coefficients: df/da, df/db, ...
end
生成GIF文件
x=findgen(101)/100.*3.1415
plot,x,sin(x),psym=1,xtitle='chen'
write_gif, 'image.gif',tvrd()
end
生成一系列GIF文件
x=findgen(101)/100.*3.1415
for iy=0,35 do begin
plot,x,sin(x),psym=1,xtitle='chen'
if iy lt 10 then name=strcompress('inflow0'+string(iy)+'.gif',/remove_all)
if iy ge 10 then name=strcompress('inflow'+string(iy)+'.gif',/remove_all)
write_gif,name,tvrd()
endfor
end
Example 5: 读取CCD数据
pro exam5
!p.background=255
!p.charsize=0.6
;device,/encapsulated,filename='t2.eps', xsize=12.,ysize=12., $
xoffset=3.2,yoffset=5.,/color, bits_per_pixel=8
w=bytarr(512,512)
openr,2,'c:\mydocu~1\sunflare.dat'
readu,2,w
close,2
window,xs=512,ys=512
tvscl,w
!p.background=1
xyouts,0.5,0.3,'Chen',alignment=0.5,size=0.8,orientation=0.,$
color=1
end
Example 6: 动态控制
pro exam6
i=findgen(10,10)
window,0,xs=300,ys=300
loadct,3
while 1 do begin
case get_kbrd(0) of
'a': tvscl,congrid(sin(i*2*3.14/50),300,300,-1)
's': surface,sin(i*2*3.14/50)
'z': return
else:
endcase
endwhile
end
Example 7: 生成EPS文件
pro exam7
!p.background=255
!p.color=0.
loadct,4
x=findgen(200) & y=findgen(200) & z=dist(200)
set_plot,'ps'
!p.charsize=1.2
device,filename='D:\t2.eps', /encapsul,xsize=12.,ysize=12., $
xoffset=3.2,yoffset=5.,/color, bits_per_pixel=8
contour,z,x,y,pos=[0.1,0.1,0.8,0.9],xst=1,yst=1, $
levels=findgen(30)/29.*150.,/fill
contour,z,x,y,pos=[0.1,0.53,0.3,0.9],xst=1,yst=1,$
levels=findgen(30)/29.*150.,/norma,/overplot,thick=2
arrow,208.,110.,208.,150.,/data,thick=3
xyouts,100.,205.,'!8t!3=130!4s!DA!3',alignment=0.5,size=2.8,$
orientation=0.,charthick=3
bar=indgen(30)/29.*150.
bar=rebin(bar,30,5)
bar=transpose(bar)
contour, bar,indgen(5),indgen(30),xst=5, yst=5,
pos=[0.82,0.1,0.87,0.4],$
levels=findgen(30)/29.*150.,/fill,/nor,/noera
xyouts,8.,-0.5,'0',alignment=0.5,size=0.8,orientation=0.
xyouts,8.,28.5,'150',alignment=0.5,size=0.8,orientation=0.
device,/close
end
Example 7: y-x plot
; plotdots.pro,v 1.0 2001/03/26 22:35:33 by P. F. Chen
; NAME: PLOTDOTS
; PURPOSE: Plot 1D profiles with the data points (x,y) read from a two-column ASCII file.
; CALLING SEQUENCE: PLOTDOTS
; OUTPUTS: A figure and a PS file.
; MODIFICATION HISTORY: March 27, 2001: modified for Windows
pro plotdots
set_plot,'win'
!p.background=255
!p.color=0.
fn=dialog_pickfile(/read)
i=0
if fn eq '' then return
openr,1,fn
desc=['0, float, 0., LABEL_LEFT=Xmin:, WIDTH=10, TAG=xmin', $
'0, float, 0., LABEL_LEFT=Ymin:, WIDTH=10, TAG=ymin', $
'0, float, 1., LABEL_LEFT=Xmax:, WIDTH=10, TAG=xmax', $
'0, float, 1., LABEL_LEFT=Ymax:, WIDTH=10, TAG=ymax', $
Continued
'0, BUTTON, dots|line, EXCLUSIVE,LABEL_TOP=Exclusive:,'+'COLUMN,TAG=bg2', $
'1, BASE,, ROW', '0, BUTTON, OK, QUIT,'+'TAG=OK', $
'2, BUTTON, Cancel, quit'+'Tag=tag7']
stru=cw_form(desc,/column)
if stru.tag7 eq 1 then return
window,xs=400,ys=400
!p.charsize=1.5
plot,[stru.xmin],[stru.ymin],psym=3,xst=1,yst=1,xrange=[stru.xmin,stru.xmax],yrange=[stru.ymin,str
u.ymax],xtitle='x',ytitle='y'
WHILE NOT EOF(1) DO BEGIN
readf,1,x,y
if stru.bg2 eq 0 then begin
oplot,[x],[y],psym=3
endif else begin
i=i+1
if i eq 1 then plots,x,y,/data
if i gt 1 then plots,x,y,/continue,/data;,psym=3
Continued
endelse
endwhile
close,1
desc2=['0, text, *.ps, LABEL_LEFT=Enter PS file name:,WIDTH=10,'+'tag=fname',$
'1, BASE,, ROW', '0, BUTTON, OK, QUIT,'+'TAG=OK', $
'2, BUTTON, Cancel, QUIT']
stru2=cw_form(desc2,/column,title='Save as a PS file')
if stru2.fname eq '*.ps' then return
tmpimg=tvrd()
set_plot,'ps'
!p.charsize=1.2
device,filename=stru2.fname, xsize=15.,ysize=23.,xoffset=3.2,yoffset=5.
tvscl,tmpimg
device,/close
end
有用的网址与有趣的命令

http://www.dfanning.com/documents/tips.html
http://cow.physics.wisc.edu/~craigm/idl
http://comp.lang.idl-pvwave/
IDL新闻组
http://idlastro.gsfc.nasa.gov/homepage.html 天文程序库
http://114.212.201.49 (南大校内网)

calendar 生成日历，如calendar,5,2003或calendar,2003







当前系统的时间: a=systime(0)或a=systime(1)
在线帮助：？
在线演示：demo
3维方案
方案之一：剖面（如左图）
方案之二：投影
方案之三：Volume Rendering
其它工具之一 — AVS
1. AVS(Advanced Visual Systems)公司
2. 高级可视化系统,多维可视化开发平台,处理海量数据
3. 广泛应用于: 工程分析、
航空航天、国防、石油工业、
地理信息与遥感、环境、
电信、有限元分析、流体力学
计算、医学、金融等方面。
4. 面向对象的可视化编程环境
5.有大量预制的可视化编程对象
组成部分
1.图形显示软件包
2.数据可视化软件包
3.图象处理软件包
4.数据库软件包
……
支持的平台：Digital UNIX, Hewlett-Packard HP-UX, IBM AIX, Silicon
Graphics IRIX（SGI）, SunOS, Sun Solaris, PC（Windows 9x and NT,
2000）, Solaris 7.0（C/C++编译器升级到了Sun Workshop 5.0）, HP-UX 11,
Compaq Tru64 UNIX 4.0E, Red Hat LINUX 6.0
其它工具之二 — Vis5D
1. Vis5D 是一个软件系统，用于将 numerical weather
models 及相似领域数据的可视化
2. Vis5D works on data in the form of a five-dimensional
rectangle. That is, the data are real numbers at each point
of a“grid” which spans three space dimensions (行、列、
层）, one time dimension and a dimension for
enumerating multiple physical variables.
3. Vis5D 第一版： by Bill Hibbard and Dave Santek （威
思康辛大学空间科学与工程中心受NASA Marshall 空
间飞行中心及法国气象局资助，1988）。
4. Vis5D is available by anonymous ftp from
www.ssec.wisc.edu/pub/vis5d
主页：http://www.ssec.wisc.edu/~billh/vis5d.html
5. 数据形式：v5d文件和comp5d文件(标准),亦可自定义数据形式
特点： (1) 不同变量可具有不同维数
(2) 提供数据转换程序：foo_to_v5d.f
foo_to_v5d.c
6. It is almost completely controlled using the mouse with a
graphical user interface. The best way to learn to use it is
to experiment.
试试
Pavo% vis5d file1.v5d
一个显示窗口对应
一个控制面板
也可直接控制显示窗
口(如旋转，zoom)
可任意输入文字、
animate、step显示动画
可截取任意部分
可用鼠标probe数据
可存成GIF(含动画）
及PS文件
可显示导入变量
其它工具之三 — VRML
Virtual Reality Modeling Language虚拟现实模拟语言
VRML is one way to share three dimensional data. The
primary advantage of VRML is that it can be transferred
and viewed using the internet and free browser plug-ins.
优点：
* 利用网络共享3维数据
* 浏览与平台无关
(PC, Mac, SGI, SUN)
* 免费浏览插件
* 第二版包含 interactive capabilities
* 文件可由纯文本编辑器产生
范例1：地球磁场
规则：
The Box Shape
范例2：太阳系
Box {
size 2.0 2.0 2.0
（注：须下载IE插件
Cosmoplayer后才能显示）
}
The size line format is:
size [width(x)] [height(y)] [depth(z)]
The Cone Shape
文件格式及编写：
Cone {
bottomRadius 1.0
height 2.0
* 纯文本，但第一行须为 #VRML V2.0 utf8 side TRUE
bottom TRUE
* 文件后缀为.wrl或.vrml
}
等等等等
•编写工具：
1. Jot, Notepad, Emacs ……
2. 三维模拟工具如cosmoplayer, 3D Studio Max
3. Fortran
例子
这是VRML文件的标志
#VRML V2.0 utf8
Group {
children [
先加入一个组节点，花括号
内的内容视为一个整体
定义组节点的孩子域
Shape {
geometry Cone {
bottomRadius 1.0
height 7.0
side TRUE
bottom TRUE
}
}
]
}
第一个孩子是一个形态节点
appearance Appearance {
material Material {diffuseColor 1 0 0 }
}