使用 IDL 8 图形功能生成椭球体
21969 为本文评分：
4.5
使用 IDL 8 图形功能生成椭球体
匿名 2014年4月10日，星期四
这篇博客文章描述了如何使用 IDL 来生成如下图所示的椭球体：
创建此图形使用了 MESH_OBJ、POLYGON 和 PLOT3D 例程，所使用的代码显示在文章底部 (dj\_ellipsoid\_ng)。
当使用 MESH_OBJ 例程创建球形对象时，您需要提供一个二维的半径值网格。该网格的行对应经度，范围从 0 到 360 度 (theta)。网格的列对应纬度，范围从 -90 到 90 度 (phi)。
由于 MESH_OBJ 需要半径值网格，而半径值取决于纬度和经度，因此必须首先生成纬度和经度网格。为了生成长度网格 (theta)，我先生成了一个包含 51 个元素的数组，范围从 0 到 360，并转换为弧度（以便与 COS 和 SIN 例程一起使用）：
theta_1d = (findgen(51)*360/50)*!dtor
然后，我结合使用了 CONGRID 和 TRANSPOSE 例程，将 theta\_1d 的值复制到一个 51 行 51 列的二维数组中。第一次调用 TRANSPOSE 处理 theta\_1d 是为了让 CONGRID 例程将其视为二维数组。
theta = congrid(transpose(theta_1d),51,51) theta = reform(transpose(theta))
我使用相同的方法创建了纬度网格 (phi)：
phi_1d = ((findgen(51)-25)*180/50)*!DTOR phi = congrid(TRANSPOSE(phi_1d),51,51)
下一步是确定每个经纬度点对应的半径值。为了计算半径，我使用了下面所示的椭球体方程：
我将这个方程转换为球坐标，然后解出 "r"：
利用这个结果，我编写了函数 dj\_get\_rad，它接受纬度和经度网格，并返回半径值网格。然后，我将此函数返回的半径值用于 MESH_OBJ，以获得椭球体的顶点和多边形数据。
最后一步是创建数据的可视化。为此，我使用 NODATA 关键字调用 PLOT3D 以生成空的 3D 数据空间。接着，我调用 POLYGON，将其 vertices 参数设置为 MESH_OBJ 返回的顶点，并将 CONNECTIVITY 关键字设置为 MESH_OBJ 返回的多边形数据。
dj\_ellispoid\_ng 例程展示了如何完成此操作。它需要 3 个参数，分别对应椭球体方程中的 "a"、"b" 和 "c" 参数。以下是一些使用不同参数值生成的示例图像：
IDL> dj_ellipsoid_ng, 2,1.5,1 ;a=2,b=1.5,c=1
IDL> dj_ellipsoid_ng, 2,2,1
IDL> dj_ellipsoid_ng, 1,1,2
function dj_get_rad, theta, phi, a,b,c
compile_opt idl2
;use the formula for the ellipsoid to
;determine the radius for various values
;of theta and phi.
radius = sqrt(((a*b*c)^2)/( $
b^2*c^2*(cos(theta)^2)*(cos(phi)^2)+ $
a^2*c^2*(sin(theta)^2)*(cos(phi)^2)+ $
a^2*b^2*(sin(phi)^2)))
return, radius
end
pro dj_ellipsoid_ng, a, b, c, data=data
compile_opt idl2
;create a grid of theta values
;these are the longitude of the of the
;ellipsoid and go from 0 to 360 degrees
theta_1d = (findgen(51)*360/50)*!dtor
theta = congrid(transpose(theta_1d),51,51)
theta = reform(transpose(theta))
;create a grid of phi values. these
;are the latitude of the ellipsoid and go
;from -90 to 90 degrees
phi_1d = ((findgen(51)-25)*180/50)*!dtor
phi = congrid(transpose(phi_1d),51,51)
;use "get_rad" function to determine the
;radius values at each point of the theta
;and phi grid
radius = dj_get_rad(theta,phi,a,b,c)
;radius = get_rad(theta,phi,a,b,c)
;create a mesh using the radius values
;this ouputs the vertices (vert) and polygons
;that will be used to generate the plot
mesh_obj,4,vert,pol,radius
;output the data
if (arg_present(data)) then begin
data=vert
endif
;determin the biggest value and create a scale from it
m = max([a,b,c])
scale = (findgen(10)-5)*m/5
;plot the scale and axis with no data. set clip=0, to
;prevent edges of ellipsoid from getting cut off.
p_scale =plot3d(scale,scale,scale,/nodata,clip=0,$
aspect_z=1, aspect_ratio=1)
;use the polygon to plot the mesh. use the
;vertices and polygons output from the mesh_obj
;to fill out the data argument and connectivity keyword.
;use bright green as the color.
p = polygon(vert,connectivity=pol,/data,$
clip=0,fill_color=[0,255,0])
end
;main level program
;for some reason, if 'a; does
;not equal 'b', the ellipsiod
;comes out messed up
dj_ellipsoid_ng,1,1.1,1.2,data=out_vert
end