科学网

 找回密码
  注册

tag 标签: 画图

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

谈谈gnuplot(十八):图例
热度 1 yusufma 2011-11-1 02:33
在同一图像中包含多组数据或函数时,图例是必要的。我们这一次谈一谈图例的微调。 这次来画前 3 阶的第一类贝塞尔函数 J n (x) 。在 gnuplot 里,0 阶和 1 阶贝塞尔函数已经定义了,分别为 besj0(x) 和 besj1(x) ,而 2 阶贝塞尔函数可以通过递推关系构造出来。下面是例子: gnuplot set term wxt enhanced gnuplot besj2(x) = besj1(x)*2/x - besj0(x) gnuplot set xrange gnuplot set xtics 2 gnuplot set xlabel "X" gnuplot set ylabel "Y" gnuplot set title "Bessel Functions of the First Kind" gnuplot set grid gnuplot set style line 1 lw 2 lc rgb "#F62217" gnuplot set style line 2 lw 2 lc rgb "#D4A017" gnuplot set style line 3 lw 2 lc rgb "#2B60DE" gnuplot plot besj0(x) ls 1 t "J_0(x)", besj1(x) ls 2 t "J_1(x)", besj2(x) ls 3 t "J_2(x)" 之前我们讲过, plot 命令后面可以跟随一些参数(例如 linewidth , linecolor 等)来改变点线风格。在上面的例子中,我们把这些参数单独拿出来放到了 set style 命令里,定义了三个 linestyle ,然后在 plot 命令里再调用这些 linestyle 。这样子做和我们之前的做法效果上没什么不同,唯一的区别是让 plot 命令短了一些。另外,改变风格可能容易一点。 上面是默认的图例,下面让我们进行微调。 为图例加上边框 gnuplot set key box gnuplot replot 改变图例显示位置 gnuplot set key center at 10,0.7 gnuplot replot 把图例的 title 和图线示例调换位置 gnuplot set key reverse gnuplot replot 调整图例边框宽度 width (或高度 height ) gnuplot set key width 1 gnuplot replot 调整 title 文字对齐方式( Left 或者 Right ,注意首字母大写) gnuplot set key Left gnuplot replot 调整图例行间隔 gnuplot set key spacing 1.2 gnuplot replot 调整图线示例长度 gnuplot set key samplen 2 gnuplot replot 这些并不是 set key 的全部参数。在 gnuplot 里,如果想深入了解任何命令的详细用法,不要忘记使用 help 命令。
个人分类: 开源软件|21528 次阅读|1 个评论
谈谈gnuplot(十七):边框和坐标轴
yusufma 2011-10-31 12:02
我们现在所有绘图的坐标刻度均标在图像边框上,无论上下左右。这样做的好处是函数或数据图线清楚,不会和坐标标注混在一起。其实,我们小时候数学课上最早学习坐标系的时候,都是让 X 轴和 Y 轴正交于原点,而刻度标注在坐标轴上。这样的图像在定性表现函数关系,尤其有一定对称性的函数关系时,比较一目了然。 让我们来看看怎样用 gnuplot 得到这样的效果。 用 unset border 命令把边框去掉; 用 set zeroaxis 命令画出正交于原点的坐标轴; 在设定坐标刻度时加上 axis 参数,这样刻度会出现在坐标轴上面,而不是边框上。 为了避免审美疲劳,这次我们拿高斯函数举个例子: gnuplot set term wxt enhanced font "Times New Roman,16" gnuplot gauss(x) = exp(-pi*x*x) gnuplot set title "函数 e^{-πx^2}" gnuplot set samples 500 gnuplot set xrange gnuplot set yrange gnuplot unset key gnuplot unset border gnuplot set zeroaxis lt -1 lw 2 gnuplot set xtics axis -2,1,2 gnuplot set ytics axis 0,1,1 gnuplot plot gauss(x) lw 3 例子中的参数前面都介绍过,如果不记得了,可以复习一下“ 坐标取值范围及刻度 ”和“ 点线风格 ”等章节。这里的图像已经很像模像样了,除了标签位置还不那么理想,而且坐标轴没有箭头。幸好,我们上一讲刚刚谈到过箭头,下面来试试做一下微调: gnuplot set title "函数 e^{-πx^2}" offset 12,-5 gnuplot set xtics axis -2,1,2 offset 0.4,0 gnuplot set ytics axis 0,1,1 offset 0,0.4 gnuplot set arrow 1 from 2,0 to 3.2,0 filled size 0.2,15,60 lw 2 gnuplot set arrow 2 from 0,1 to 0,1.22 filled size 0.2,15,60 lw 2 gnuplot set rmargin 4 gnuplot set label 1 "X" at 3.0,-0.1 gnuplot set label 2 "Y" at -0.3,1.2 gnuplot replot 这里有几个命令同时用到了新的参数: offset 。它的作用就是把命令里提到的标签文字平移一段距离。在这里, offset 默认的坐标系统是 character 。我们慢慢会体会到这种做法的好处,它使得我们很多时候改变字体大小,而不必重新设置 offset 。 另外, set rmargin 命令用于设置图像右边空白宽度,单位也是 character 。一般情况下,四边空白宽度都是自动设置的。现在我们在右边增加了箭头,而绘图显示区域不会因此自动扩大,这样会导致箭头无法完整显示,所以要手动改一下设置。相应的,上、左、下边的空白宽度,分别由 tmargin , lmargin , bmargin 参数控制。
个人分类: 开源软件|21028 次阅读|0 个评论
谈谈gnuplot(十六):箭头
yusufma 2011-10-29 02:32
有了坐标系的知识打底,其他很多东西很好谈了。我们的图上除了标签之外,还有一个常用的标志:箭头。关于箭头的命令是 set arrow ,语法和 label 有些类似,包括以下这些常用参数: from ... to ... 箭头的起点和终点坐标。如果把 to 换成 rto ,第二个坐标就表示相对位置而不是绝对坐标。 nohead, head, backhead, heads 分别表示:没有箭头(其实就是线段),箭头在终点,箭头在起点,双向都有箭头。 size length,angle,backangle 箭头尺寸,默认长度单位为 first 坐标单位长度。 下图中 A,B,C 分别代表 length , angle , backangle 。 filled, empty, nofilled 箭头的三种填充风格: 下面我们看例子,还是画 sinc(x) 函数: gnuplot set term wxt font "DejaVu Sans,12" gnuplot sinc(x) = sin(pi*x)/(pi*x) gnuplot set xlabel "X" gnuplot set ylabel "Y" gnuplot set yrange gnuplot set title "sinc(x) 函数" gnuplot unset key gnuplot set samples 500 gnuplot set arrow 1 from 2,1.05 to 0.3,1 filled size 0.5,15,60 lw 2 gnuplot set label 1 at 0,1 point pt 7 ps 1.5 lc rgb "#F87217" gnuplot set label 2 "最大值在(0, 1)" at 2.1,1.05 gnuplot plot sinc(x) lw 2
个人分类: 开源软件|7935 次阅读|0 个评论
[转载]matlab动态画图(2)
onewaystreet 2011-10-28 20:47
Matlab除了强大的矩阵运算,仿真分析外,绘图功能也是相当的强大,静态画图没什么问题,由于Matlab本身的多线程编程缺陷,想要动态的画图,并且能够很好的在GUI中得到控制,还不是一件很容易的事情,下面总结几种方法。 一. AXIS 移动坐标系 这种方法是最简单的一种方法,适合于数据已经全部生成的场合,先画图,然后移动坐标轴。实例代码如下: %% %先画好,然后更改坐标系 %在命令行中 使用 Ctrl+C 结束 t=0:0.1:100*pi; m=sin(t); plot(t,m); x=-2*pi; axis( ); grid on while 1 if xmax(t) break; end x=x+0.1; axis( ); %移动坐标系 pause(0.1); end 二. Hold On 模式 此种方法比较原始,适合于即时数据,原理是先画上一帧,接着保留原始图像,追加下一幀图像,此种方式比较繁琐,涉及画图细节,并且没有完整并连续的Line对象数据。 例如: %% % Hold On 法 % 此种方法只能点,或者 分段划线 hold off t=0; m=0; t1= ; %要构成序列 m1= ; p = plot(t,m,'*',t1,m1(1,:),'-r',t1,m1(2,:),'-b','MarkerSize',5); x=-1.5*pi; axis( ); grid on; for i=1:100 hold on t=0.1*i; %下一个点 m=t-floor(t); t1=t1+0.1; %下一段线(组) m1= ; p = plot(t,m,'*',t1,m1(1,:),'-r',t1,m1(2,:),'-b','MarkerSize',5); x=x+0.1; axis( ); pause(0.01); end 三. Plot 背景擦除模式 这种模式比较适合画动画,效率比较高,刷新闪烁小,适合即时数据,最终的Line结构数据完整。 了解此方法之前要搞清楚 Plot函数的原型是什么: Plot函数,输入为 X-Y (-X)坐标元组、以及“属性”-“值对,输出为一个列向量(每条曲线岁对应的Line结构 Handle,每一行代表一个 线条的handles), 每一线条都有 XData,YData 向量。如果你画了2条线,那么会返回 2×1的向量。 重新画图不需要 重新书写 Plot,只需要 刷新图像即可,使用drawnow函数。 完整实例如下: 1. 画一个点的动画: %% %采用背景擦除的方法,动态的划点,并且动态改变坐标系 % t,m 均为一行 ,并且不能为多行 t=0; m=0; p = plot(t,m,'*',... 'EraseMode','background','MarkerSize',5); x=-1.5*pi; axis( ); grid on; for i=1:1000 t=0.1*i;%两个变量均不追加 m=sin(0.1*i); set(p,'XData',t,'YData',m) x=x+0.1; drawnow axis( ); pause(0.1); end 2. 动态多条曲线(即时数据) %% %采用背景擦除的方法,动态的划线,并且动态改变坐标系 % 多行划线 t= m= p = plot(t,m,... 'EraseMode','background','MarkerSize',5); x=-1.5*pi; axis( ); grid on; for i=1:1000 t= ;%Matrix 1*(i+1) m= ]; %Matrix 2*(i+1) set(p(1),'XData',t,'YData',m(1,:)) set(p(2),'XData',t,'YData',m(2,:)) drawnow x=x+0.1; axis( ); pause(0.5); end 上面的这几个画图方式的示例只是简单的for循环,是单线程的,如果是涉及到GUI的编程,那么请使用Timer来完成这件事情,Timer是我在Matlab中实现多线程唯一方法(没有找到别的方法)。
个人分类: matlab|4074 次阅读|0 个评论
[转载]matlab动态画图
onewaystreet 2011-10-28 20:45
一般来说,matlab制作动画有四种方式。 第一 、 以质点运动轨迹的方式显示 使用comet、comet3函数,前者是二维,后者是三维 comet(y)显示质点绕向量y,comet(x,y)显示质点绕向量y与x,comet(x,y,p),其中为轨迹尾巴的长度 以comet(x,y)为例, 显示平抛运动 vx = 40; t = 0:0.001:10; x = vx*t; y = -9.8*t.^2/2; comet(x,y) 显示导弹发射 vx = 100*cos(1/4*pi); vy = 100*sin(1/4*pi); t = 0:0.001:15; x = vx*t; y = vy*t-9.8*t.^2/2; comet(x,y) 匀速圆周运动 sita = 0:0.0001:2*pi; r = 10; x=r*cos(sita); y=r*sin(sita); comet(x,y) comet3与comet的用法相类似,可以在帮助文件里的例子 t = -10*pi:pi/250:10*pi; comet3((cos(2*t).^2).*sin(t),(sin(2*t).^2).*cos(t),t) 第二、以电影播放的方式显示 保存想要产生动画的图片,存储为一系列各种类型的二维、三维图,再像放电影的方式按次序播放出来。步骤由getframe函数将当前的图片抓取为电影的画面,再由movie函数将动画显示出来。 如: = meshgrid( ); z = x.*exp(-x.^2-y.^2); axis tight; set(gca,'nextplot','replacechildren'); for j = 1:40 surf(x*sin(pi*j/100),y*sin(pi*j/100),z*sin(-pi*j/100)); m(j) = getframe end movie(m) 第三、以对象方式显示 设置对象的属性EraseMode,更新对象来产生新图,drawnow()函数进而覆盖旧图,从而使得图形不断发生变化。 例: x = -pi:pi/30:pi; h = plot(x,cos(x),'o','MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',8,'EraseMode','Xor') for j = 1:10000 y = 1/2*sin(3*x+0.006*j); set(h,'ydata',y); drawnow; end 第四、以旋转颜色的方式显示
个人分类: matlab|15355 次阅读|0 个评论
谈谈gnuplot(十五):gnuplot 的坐标系统及标签
热度 1 yusufma 2011-10-28 07:10
我们现在知道了 gnuplot 有第一( first )和第二( second )两套坐标系统,但是 gnuplot 的坐标系统还不止于此。除此之外,它还有 graph , screen 和 character 三套坐标系统。 graph 和 screen 都是归一化的坐标系统。 graph 以坐标轴包围区域为界,左下角为 0,0 ,右上角为 1,1 ; screen 以整个图片区域为界,左下角为 0,0 ,右上角为 1,1 。 character 顾名思义,是以字符大小为单位长度的坐标系统,因此它的单位长度依赖于字体大小。它的原点位置和 screen 相同。 下面我们结合 label 命令来了解一下这几个坐标系统。我们之前讲过 xlabel 和 ylabel 。而这里的 label 命令,是在图中任何地方插入文字标签。还是来看例子: gnuplot sinc(x) = sin(pi*x)/(pi*x) gnuplot set xlabel "X" gnuplot set ylabel "Y" gnuplot unset key gnuplot set samples 500 gnuplot set xrange gnuplot set xtics 1 gnuplot set x2range gnuplot set x2tics 1 gnuplot set y2range gnuplot set y2tics 1 gnuplot set grid gnuplot set label 1 "Hello first" at 2,0.5 gnuplot set label 2 "Hello second" at second 2,0.5 gnuplot set label 3 "Hello graph" at graph 0.2,0.5 gnuplot set label 4 "Hello screen" at screen 0.2,0.5 gnuplot set label 5 "Hello character" at character 10,5 gnuplot plot sinc(x) 这里我们画一个 sinc 函数图像。为了说明问题,我们把第二坐标系也都标示了出来,虽然函数图像并没有用到第二坐标。其他命令前面都讲过了,这里只看五个 set label 命令。 set label 之后紧跟的那个整数,就是一个标识符,用以区别各个 label ,可以随便选个整数。在字符串之后, at 参数指定标签坐标。默认为 first 坐标系统,也可以使用其它坐标系统。下面是生成的图片: 为了帮助大家理解,我们把 graph 和 screen 各自的坐标区域分别用绿色和橙色表示了出来。 标签文字的默认对齐方式为居左,也就是指定的坐标位置在文字的左边。我们也可以在 label 命令里选择其他对齐方式。除此之外,我们还可以在 label 命令里指定文字颜色,旋转文字,或者在指定坐标位置处加一个点。下面例子中的每个参数不必一一解释了,因为和我们前面接触过的命令都是一致的: gnuplot set label 1 "Hello red left" at 2,0.4 left textcolor rgb "#FF0000" gnuplot set label 2 "Hello green center" at 2,0.5 center textcolor rgb "#00FF00" gnuplot set label 3 "Hello blue right" at 2,0.6 right textcolor rgb "#0000FF" gnuplot set label 4 "Hello rotate" at -2,0.4 rotate by 45 gnuplot set label 5 "Hello point" at -3,0.2 point pt 7 lc rgb "#FF9900" gnuplot replot
个人分类: 开源软件|13039 次阅读|1 个评论
谈谈gnuplot(十四):第二坐标轴
yusufma 2011-10-27 06:46
回首看看我们以前所有的作图,横坐标都标示在底部,而纵坐标都标示在左侧。其实,在图像顶部和右侧,还隐藏着一对不太引人注意的坐标轴,我们可以管它们叫做“第二坐标轴”。平时,它们只是第一对坐标轴的镜像;在我们需要的时候,它们可以用来表示不同的物理量。有时候,我们会有两组性质不同但是又相互关联的数据,这时候我们或许想把他们画在同一副图上,以便比较。 还拿北京市月平均降水量举例,但是这次,我们把温度也加上。下面是我们的数据文件 weather_beijing.dat : ### 文件开始 ### # 北京月平均降水量(毫米)及气温(摄氏度) # # 月份 降水量 气温 # =================== 1 2.5 -4 2 5.1 -2 3 10.2 6 4 25.4 13 5 27.9 20 6 71.1 24 7 175.3 26 8 182.9 25 9 48.3 20 10 17.8 13 11 5.1 5 12 2.5 -2 ### 文件结束 ### 我们之前讲过的所有有关坐标的参数,在第二坐标轴上均适用,只不过相应的名字起始字母改为 x2 或者 y2 ,例如 ylabel 改为 y2label 。另外, plot 命令有一个新的参数 axis ,用来控制使用哪个坐标轴,例如 axis x1y2 就表示使用第一横轴和第二纵轴。现在我们来看用上面数据作图的例子: gnuplot set xlabel "月份" gnuplot set ylabel "降水量(毫米)" gnuplot set y2label "气温(摄氏度)" gnuplot set title "北京市月平均降水量及气温" gnuplot set xrange gnuplot set xtics 1,1,12 gnuplot plot "weather_beijing.dat" u 1:2 w lp pt 5 lc rgbcolor "#2B60DE" axis x1y1 t "降水量", "weather_beijing.dat" u 1:3 w lp pt 7 lc rgbcolor "#F62817" axis x1y2 t "气温" 这里的气温数据使用了图像右边的第二纵轴 y2 ,但是 y2 轴上的刻度并没有变化,依然是左边 y1 轴的镜像。我们在这里有两件事要做: 去除右边纵轴上的 y1 刻度镜像,否则这些刻度标记将和新的 y2 刻度标记混起来,导致无法识别; 在右边纵轴上加上 y2 刻度标记。 我们执行下面的命令: gnuplot set ytics nomirror gnuplot set y2tics gnuplot replot 好了,现在降水量和温度数据分别对应于左侧和右侧的纵坐标。 看到这里,我们可能有点怀念我们上一讲谈到的 grid 。如果能加上栅格,数据图示就更清楚了。但是现在我们有两组不同的纵坐标,如果都开启栅格,还不乱套了? set grid 命令允许我们在开启栅格时,选择使用哪一组坐标。例如: gnuplot set grid xtics y2tics 会开启 x1 和 y2 的栅格。但是这还是不能兼顾两组数据。最好的解决方案是,让两个纵轴有相同数目的分格,这样两套 grid 也就重合了,开启任何一个就可以了。例如,上面的图中左侧纵轴有 10 个分格,我们让右侧纵轴也有 10 个分格: gnuplot set y2range gnuplot set y2tics 5 gnuplot set grid gnuplot replot 现在看起来好多了。 最后,不知道大家注意到没有,在开始的 plot 命令里,我们用了新的方式定义图线颜色。在 第七讲“点线风格” 里,我们提到过,可以用预定义的数字代码来定义图线颜色。但是在这里,我们使用了 rgbcolor 来定义颜色,这很大程度上增加了颜色选择范围,允许更好的显示效果。而其用法也很简单,就是在 rgbcolor 之后,加上颜色的 RGB 代码,了解 HTML 的朋友应该对这个不陌生。
个人分类: 开源软件|10094 次阅读|0 个评论
谈谈gnuplot(十三):栅格以及方程数值解估算
yusufma 2011-10-26 03:18
我们现在来画一个 0 阶贝塞尔函数 J 0 (x) : gnuplot set term wxt enhanced gnuplot set xlabel "X" gnuplot set ylabel "Y" gnuplot set xrange gnuplot set xtics 0,1,10 gnuplot unset key gnuplot set title "0阶贝塞尔函数 J_0(x)" gnuplot plot besj0(x) 这里的 besj0(x) 就是 gnuplot 里面预定义的 0 阶贝塞尔函数。如果现在请您从这个图上估计出 内 J 0 (x) 的零点数值,也就是方程 J 0 (x)=0 的解,恐怕您很难说的准确。但是如果为这个图加上栅格( grid ),就容易多了: gnuplot set grid gnuplot replot 这时我们很容易估计出三个零点的数值: 2.4, 5.5, 8.6 。通过查表我们可以知道,这三个零点比较精确的数值分别为 2.4048, 5.5201, 8.6537 。这和我们的估计值差不太多。如果我们想更精确的估计数值,可以尝试改一下 xrange : gnuplot set xrange gnuplot set xtics 8,0.1,9 gnuplot replot 这相当于把图像在零点附近放大了。把鼠标放在画图区域,画图框左下角就会显示出鼠标所在位置的坐标。现在我们把鼠标放在函数图线和 X 轴的交叉点上,左下角显示的横坐标为 8.65243 ,这和我们查表所得的数值更接近了。 如果想进一步让结果精确一些,我们可以利用 gnuplot 的计算功能。我们可以通过尝试计算的方法获得方程的数值解: gnuplot print besj0(8.65) 0.00101216621937318 gnuplot print besj0(8.66) -0.0017019446057587 gnuplot print besj0(8.6537) 7.5770361108123e-06 gnuplot print besj0(8.6536) 3.47225104115535e-05 gnuplot print besj0(8.6538) -1.95681245811775e-05 所以在 8.6 附近, J 0 (x)=0 精确到小数点后 4 位的数值解为 8.6537 ,这和我们查表的结果一模一样。由于我们已经通过图像知道了数值解的大概位置,再加上合理利用线性插值,我们可以很快得到精确的结果。
个人分类: 开源软件|6941 次阅读|0 个评论
谈谈gnuplot(十二):插入 LaTeX 公式
热度 4 yusufma 2011-10-25 03:32
上一次我们谈过,在 gnuplot 里使用 enhanced 模式虽然可以生成一些简单的数学表达式,但是对于稍复杂的数学公式来说, enhanced 模式没办法生成令人满意的结果。这里我们介绍 gnuplot 的另外一个 terminal : epslatex 。 epslatex 和我们之前介绍过的 postscript eps 输出方式非常接近,因此它们很多参数都是相同的。区别在于, epslatex 使用 postscript eps 仅生成图形存于 eps 文件,而所有文字标签包含在另外一个 LaTeX 文件中。在 gnuplot 完成输出之后,使用 LaTeX 命令最终生成完整的图片。这种做法的好处是不言而喻的,即使在输出完成后,我们仍然可以编辑 LaTeX 文件获得我们想要的显示效果。 下面我们看例子: gnuplot set xlabel 'X' gnuplot set ylabel 'Y' gnuplot set title 'Error function $\displaystyle\mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2}\,\mathrm{d}t$' gnuplot set xrange gnuplot set yrange gnuplot unset key gnuplot set term epslatex standalone lw 2 color 11 gnuplot set output "erf.tex" gnuplot plot erf(x) lw 2 gnuplot set output 前三行的标签文字,我们使用了单引号,避免了双引号内需要连着两个反斜杠的麻烦。在 title 里面,我们使用了 LaTeX 数学公式。在 set term 命令里, standalone 是一个新的参数,它表示生成完整的 LaTeX 文件。如果没有这个参数,生成的 LaTeX 文件将不能单独编译,必须把代码插入其它的 LaTeX 文件中编译。 set term 最后的数字 11 代表字体大小。 set output 只需指定 LaTeX 文件名,而不需要指定 eps 文件名。我们谈到过 gnuplot 包含很多预定义的函数,这里的 erf 就是其中之一,表示误差函数。 我们通过 epstopdf 命令将生成的 eps 图片转为 pdf,然后用 pdflatex 命令可以把图片直接编译为 pdf 文件,下面是显示的效果: 怎么样?LaTeX 的数学公式效果真不是盖的。
个人分类: 开源软件|13725 次阅读|11 个评论
谈谈gnuplot(十):输出 pdf 和 png 图片
热度 2 yusufma 2011-10-21 02:10
这次讲讲怎样把图片输出为 pdf 和 png 格式。 上次讲过的 eps 文件其实很容易转换为 pdf,为什么我们还专门讲一下 pdf 格式输出呢?因为在 postscript terminal 下很难使用中文字体,而在 pdf 下面就容易多了,而 pdf 文件又很容易转换回 eps,这样就绕开了 eps 中文支持不好的问题。 png 是一种无损压缩位图格式,图形质量优于 jpg 等有损压缩格式,支持透明效果,可以生成非常小体积的文件,适于放在网上交流。通过各种图像处理软件,png 文件也很容易转换为其他位图格式。 下面首先看一个 pdf 输出的例子,咱们还是用之前用过的城市降水量数据文件: gnuplot set xlabel "月份" gnuplot set ylabel "降水量(毫米)" gnuplot set title "各城市月平均降水量" gnuplot set xrange gnuplot set xtics 1,1,12 gnuplot set term pdfcairo lw 2 font "Times New Roman,8" gnuplot set output "precipitation.pdf" gnuplot plot "precipitation.dat" u 1:2 w lp pt 5 title "北京", "precipitation.dat" u 1:3 w lp pt 7 title "上海" gnuplot set output 这里我们用的 terminal 是 pdfcairo ,而不是简单的 pdf 。区别是 pdfcairo 使用了 cairo(一个2D图形程序库)和 pango(一个字体渲染程序库)来生成 pdf 文件,优点是更好的国际支持。有了之前的经验,这里的 terminal 参数不需要多解释了。这里我们使用了“Times New Roman” 8号字体。和 eps 下使用 postscript 字体不同,这里可以是电脑系统里安装的任何字体。在 Linux 下,可以使用 fc-list 命令察看系统里到底有哪些字体可用。 下面我们来看生成的 pdf 图片: 这里有一个小问题:虽然数字使用了Times New Roman字体,但是汉字使用了其他字体(这里是我的系统默认的“文泉驿正黑”)。这是因为Times New Roman本来就不是中文字体。如果我们想让中英文混排时字体统一,必须使用支持中文的字体。 下面我们来看 png 输出的例子: gnuplot set term pngcairo lw 2 font "AR PL UKai CN,14" gnuplot set output "precipitation.png" gnuplot replot gnuplot set output gnuplot set term wxt 基于和上面同样的原因,这里使用的 terminal 是 pngcairo 而不是简单的 png ,而字体是 AR PL UKai CN (文鼎简中楷)。下面是生成的 png 图片:
个人分类: 开源软件|27855 次阅读|3 个评论
谈谈gnuplot(七):点线风格
yusufma 2011-10-18 04:25
我们接着上次的数据图谈起。上次我们得到了这样一个“点线”图: 这里的数据点是由小“十”字表示的,但是似乎太小了,有点看不清楚。另外,如果我们想在做报告时把这个图用到幻灯片中去,小“十”字很不醒目,这时候我们可能想用其他的标志。gnuplot里面有几个控制点和线画法风格的参数: linestyle 连线风格(包括 linetype , linewidth 等) linetype 连线种类 linewidth 连线粗细 linecolor 连线颜色 pointtype 点的种类 pointsize 点的大小 我们看下面的例子: gnuplot plot "datafile.dat" with linespoints linecolor 3 linewidth 2 pointtype 7 pointsize 2 这几个参数的用法不难理解,直接跟在 with 命令之后就可以了,但是2、3、7这些数字都代表什么意思呢?这些数字是代表不同画法风格的代码,具体某个数字代表什么意思,这个依赖于我们使用的 terminal (还记得我们在第二讲里曾经讲过的 terminal 吗?)拿我们现在正在使用的 wxt terminal 举例,如果想知道这些数字究竟代表什么意思,可以输入命令: gnuplot test 这样当前 terminal 会输出一个测试图: 测试图中包含当前 terminal 的风格代码实例。例如,左下角显示的是连线粗细,右边显示的是色彩和数据点显示风格对应代码。 最后,告诉大家一个好消息:gnuplot里面很多命令有缩写形式。例如上面例子中的绘图命令可以简写为: gnuplot plot "datafile.dat" w lp lc 3 lw 2 pt 7 ps 2 至于其他参数命令的缩写形式,相信不难猜出来,大家可以试验一下猜猜看。
个人分类: 开源软件|22809 次阅读|0 个评论
谈谈gnuplot(三):数学表达式
yusufma 2011-10-14 01:38
在我们开始画图之前,需要知道gnuplot里面是如何表达数学公式的。 加、减、乘、除、乘方 分别用 + , - , * , / , ** 表示 整数和浮点数 和C语言类似,gnuplot对整数和浮点数(实数)区别对待,整数的运算结果还是整数。所以在处理整数除法时要尤其小心,例如7/2的结果是3而不是3.5 复数 gnuplot支持复数运算,复数用包含在花括号内的一对实数表示,例如 {3,5} 表示 3+5i 数学函数 gnuplot含有丰富的数学函数,格式和C语言几乎相同。对于实数和复数,函数名是一样的。下面的链接可以看到预定义的函数列表: http://www.gnuplot.info/docs_4.2/gnuplot.html#x1-5300013.1 自定义函数 自定义函数很容易,例如 f(x)=x+1 定义一个一元函数, f(x,y)=x+y 定义一个二元函数。 π(圆周率) π 在gnuplot里用 pi 表示。 这里是一些例子: 这里用到了 print 命令,就是把结果输出到屏幕上。 有了这些知识做准备,我们就可以正式开始画图了。
个人分类: 开源软件|9981 次阅读|0 个评论
[转载]图说“博士”研究
liuj589 2010-12-20 16:51
什么是博士? 原文地址: http://www.ruanyifeng.com/blog/2010/08/illustrated_guide_to_a_phd.html 今天,我看到 一组图 。 美国犹他大学的助理教授 Matt Might ,用这组图解释,博士学位到底是什么意思。他说,每年都有新生的入学教育,但是有些观点语言说不清楚,不如画图。 我觉得,这组图真的很好懂,而且一点没错,博士就应该是图中的意思。老子说大道至简,可是真的要很简单地表达出来,却是非常难的一件事。 ======================================== 什么是博士? 作者:Matt Might 译者:阮一峰 1. 假设人类所有的知识,就是一个圆。圆的内部代表已知,圆的外部代表未知。 2. 读完小学,你有了一些最基本的知识。 3. 读完中学,你的知识又多了一点。 4. 读完本科,你不仅有了更多的知识,而且还有了一个专业方向。 5. 读完硕士,你在专业上又前进了一大步。 6. 进入博士生阶段,你大量阅读文献,接触到本专业的最前沿。 7. 你选择边界上的一个点,也就是一个非常专门的问题,作为自己的主攻方向。 8. 你在这个点上苦苦思索,也许需要好几年。 9. 终于有一天,你突破了这个点。 10. 你把人类的知识向前推进了一步,这时你就成为博士了。 11. 现在你就是最前沿,其他人都在你身后。 12. 但是,不要陶醉在这个点上,不要把整张图的样子忘了。 继续努力向前推进吧! (完)
个人分类: 重要新闻|2582 次阅读|0 个评论
小记IDL画全天中性氢分布图
qianlivan 2010-10-15 17:34
和地理学差不多,天文的研究对象通常也可看成分布在一个球面上,此球面即为天球。天球的半径并不重要,定为多少其实都没有关系,教科书上说是无穷大,也并非不可,但是此说法通常误导人。在描述天体的时候,如果不需要考虑距离,就把它们投影到天球上(就是取其球坐标中的俯仰角和方位角,而不考虑半径,这才是本质)。 通常我们的书和电脑屏幕都是平面,我们经常需要在上面表示整个世界或者整个天球。不过,即使从拓扑学的角来讲,平面(特别是有限大小的平面区域)和球面也有本质的不同。所以如果要在一个平面区域上表示出球面的一部分或者整个球面就需要一些变换,就是投影。 关于天文里常用的FITS文件中对投影的规范,可以参考http://www.cv.nrao.edu/fits/documents/wcs/wcs.html,里面有一篇文章专门叙述了FITS文件所接受的投影变换的细节。 这里想记的其实只是怎么把一个包含了全体中性氢数据的文件里的数据画出来,并标上坐标。首先要做的是读出数据,目前比较好的一个全体的中性氢分布数据是LAB巡天的数据。假设文件名为lab.fit,那么首先可以用mrdfits读出 IDLa=mrdfits(lab.fit) 这时a是一个三维数组,为了画全天的分布图,可以画一个频率通道,也可以画总强度。假设画总强度 IDLb=total(a,3) 下面就是画图了 IDLLoadCT, 39, NColors=colors, Bottom=1, /Silent IDLdevice,decompose=0 这是定义颜色,然后对数据作投影,画坐标线 IDLmap_set,0,0,/MOLLWEIDE,/ISOTROPIC,/HORIZON,/GRID IDLresult=map_image(b,Startx,Starty,Xsize,Ysize,compress=1,LATMIN=latmin,LONMIN=lonmin,LATMAX=latmax,LONMAX=lonmax,scale=0.1) IDLresult=bytscl(result) IDLtvscl,result,Startx,Starty,XSIZE=Xsize,YSIZE=Ysize IDLlons=indgen(360/45+1)*45-180 IDLlonnames=strtrim(-lons) IDLmap_grid,latdel=10,lonnames=lonnames,lons=lons,color=0.30*!d.n_colors,/LABEL,/HORIZON 其中map_image做的事就是对数据作投影,这里用的是Mollweide投影(用map_set设置/MOLLWEIDE), 两个重要的返回值是Xsize和Ysize ,如果不将这两个值传递给tvscl,那么在 PS图 里就无法将坐标线和投影图对齐,但是对于Xwindow的输出没有影响,原因是PS图像素的大小是可以变的。在标注坐标线的经度的时候要注意,星图和地图正好是反的,所以有lonnames=strtrim(-lons)一句。 先草记这些,要画一幅漂亮的全天分布图其实有很大学问,等我会了再写。
个人分类: 总结|9656 次阅读|0 个评论
[转载]Matlab学习笔记——matlab画图指令
arcechen 2010-8-20 14:52
1.基本画图 程序如下: view plaincopy to clipboardprint? x=0:pi/1000:2*pi; y1=sin(2*x); y2=2*cos(2*x); %输出图像 plot(x,y1,'k-',x,y2,'b--'); title(' Plot of f(x)=sin(2x) and its derivative'); %设置X坐标和Y坐标的标签 xlabel('x'); ylabel('y'); %制作图例 legend('f(x)=sin(2x)','d/dx f(x)') %显示网格 grid on; x=0:pi/1000:2*pi; y1=sin(2*x); y2=2*cos(2*x); %输出图像 plot(x,y1,'k-',x,y2,'b--'); title(' Plot of f(x)=sin(2x) and its derivative'); %设置X坐标和Y坐标的标签 xlabel('x'); ylabel('y'); %制作图例 legend('f(x)=sin(2x)','d/dx f(x)') %显示网格 grid on; 2.利用有限个点画出平滑曲线 view plaincopy to clipboardprint? x=1:1:10; y= ; xx=1:0.01:10; %使用了函数interp1 yy=interp1(x,y,xx,'cublic'); plot(xx,yy,x,y,'.'); grid on; x=1:1:10; y= ; xx=1:0.01:10; %使用了函数interp1 yy=interp1(x,y,xx,'cublic'); plot(xx,yy,x,y,'.'); grid on; 3.绘制三维曲线 view plaincopy to clipboardprint? t=0:pi/100:20*pi; x=sin(t); y=cos(t); z=t.*sin(t).*cos(t); plot3(x,y,z); title('Line in 3-D Space'); xlabel('X');ylabel('Y');zlabel('Z'); grid on; t=0:pi/100:20*pi; x=sin(t); y=cos(t); z=t.*sin(t).*cos(t); plot3(x,y,z); title('Line in 3-D Space'); xlabel('X');ylabel('Y');zlabel('Z'); grid on; 4.绘制三维曲面 view plaincopy to clipboardprint? =meshgrid(-8:0.5:8); z=sin(sqrt(x.^2+y.^2))./sqrt(x.^2+y.^2+eps); subplot(2,2,1); mesh(x,y,z); title('mesh(x,y,z)') subplot(2,2,2); meshc(x,y,z); title('meshc(x,y,z)') subplot(2,2,3); meshz(x,y,z) title('meshz(x,y,z)') subplot(2,2,4); surf(x,y,z); title('surf(x,y,z)') =meshgrid(-8:0.5:8); z=sin(sqrt(x.^2+y.^2))./sqrt(x.^2+y.^2+eps); subplot(2,2,1); mesh(x,y,z); title('mesh(x,y,z)') subplot(2,2,2); meshc(x,y,z); title('meshc(x,y,z)') subplot(2,2,3); meshz(x,y,z) title('meshz(x,y,z)') subplot(2,2,4); surf(x,y,z); title('surf(x,y,z)') 标准三维曲面: view plaincopy to clipboardprint? t=0:pi/20:2*pi; = cylinder(2+sin(t),30); subplot(2,2,1); surf(x,y,z); subplot(2,2,2); =sphere; surf(x,y,z); subplot(2,1,2); =peaks(30); surf(x,y,z); t=0:pi/20:2*pi; = cylinder(2+sin(t),30); subplot(2,2,1); surf(x,y,z); subplot(2,2,2); =sphere; surf(x,y,z); subplot(2,1,2); =peaks(30); surf(x,y,z); 其他函数: view plaincopy to clipboardprint? subplot(2,2,1); bar3(magic(4)) subplot(2,2,2); y=2*sin(0:pi/10:2*pi); stem3(y); subplot(2,2,3); pie3( ); subplot(2,2,4); fill3(rand(3,5),rand(3,5),rand(3,5), 'y' ) subplot(2,2,1); bar3(magic(4)) subplot(2,2,2); y=2*sin(0:pi/10:2*pi); stem3(y); subplot(2,2,3); pie3( ); subplot(2,2,4); fill3(rand(3,5),rand(3,5),rand(3,5), 'y' ) 绘制多峰函数的瀑布图和等高线图: view plaincopy to clipboardprint? subplot(1,2,1); =peaks(30); waterfall(X,Y,Z) xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis'); subplot(1,2,2); contour3(X,Y,Z,12,'k'); %其中12代表高度的等级数 xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis'); subplot(1,2,1); =peaks(30); waterfall(X,Y,Z) xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis'); subplot(1,2,2); contour3(X,Y,Z,12,'k'); %其中12代表高度的等级数 xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis'); 5.图形修饰处理 三种图像着色方式效果显示: view plaincopy to clipboardprint? =sphere(20); colormap(copper); subplot(1,3,1); surf(x,y,z); axis equal subplot(1,3,2); surf(x,y,z);shading flat; axis equal subplot(1,3,3); surf(x,y,z);shading interp; axis equal =sphere(20); colormap(copper); subplot(1,3,1); surf(x,y,z); axis equal subplot(1,3,2); surf(x,y,z);shading flat; axis equal subplot(1,3,3); surf(x,y,z);shading interp; axis equal 光照处理后的球面: view plaincopy to clipboardprint? =sphere(20); subplot(1,2,1); surf(x,y,z);axis equal; light('Posi', ); shading interp; hold on; plot3(0,1,1,'p');text(0,1,1,' light'); subplot(1,2,2); surf(x,y,z);axis equal; light('Posi', ); shading interp; hold on; plot3(1,0,1,'p');text(1,0,1,' light'); =sphere(20); subplot(1,2,1); surf(x,y,z);axis equal; light('Posi', ); shading interp; hold on; plot3(0,1,1,'p');text(0,1,1,' light'); subplot(1,2,2); surf(x,y,z);axis equal; light('Posi', ); shading interp; hold on; plot3(1,0,1,'p');text(1,0,1,' light'); 图形的裁剪处理: view plaincopy to clipboardprint? =meshgrid(-5:0.1:5); z=cos(x).*cos(y).*exp(-sqrt(x.^2+y.^2)/4); subplot(1,2,1); surf(x,y,z);shading interp; title('裁剪之前'); i=find(x=0y=0); z1=z;z1(i)=NaN; subplot(1,2,2); surf(x,y,z1);shading interp; title('裁剪之后'); =meshgrid(-5:0.1:5); z=cos(x).*cos(y).*exp(-sqrt(x.^2+y.^2)/4); subplot(1,2,1); surf(x,y,z);shading interp; title('裁剪之前'); i=find(x=0y=0); z1=z;z1(i)=NaN; subplot(1,2,2); surf(x,y,z1);shading interp; title('裁剪之后'); 从文件载入图像: view plaincopy to clipboardprint? =imread('flower.jpg'); %读取图像的数据阵和色图阵 image(x);colormap(cmap); axis image off %保持宽高比并取消坐标轴 =imread('flower.jpg'); %读取图像的数据阵和色图阵 image(x);colormap(cmap); axis image off %保持宽高比并取消坐标轴 制作动画: view plaincopy to clipboardprint? =peaks(30); surf(X,Y,Z) axis( ) axis off; shading interp; colormap(hot); m=moviein(20); %建立一个20列大矩阵 for i=1:20 view(-37.5+24*(i-1),30) %改变视点 m(:,i)=getframe; %将图形保存到m矩阵 end movie(m,2); %播放画面2次 注意 在同时绘制多条曲线时,如果没有指定曲线属性,plot按顺序循环使用当前坐标系中ColorOrder和LineStyleOrder两个属性。 默认情况,MATLAB在每次调用plot函数时将ColorOrder和LineStyleOrder自动重置为DefaultAxesColorOrder和DefaultAxesLineStyleOrde r。Default**属性我们可以自定义,有效期至MATLAB关闭,Matlab下次启动时将Default**属性重置为厂家设置(Factory) set(0,'DefaultAxesColorOrder',r|g|b|k,... 'DefaultAxesLineStyleOrde r','-|-.|--|:') 使用hold all命令可以阻止调用plot函数时自动重置ColorOrder和LineStyleOrder属性,而是循环使用。注意hold on只是使多次绘制的图形叠加(相当于NextPlot),但不能阻止属性重置。 另外我们可以通过下面四个属性设置标识符的颜色和大小 LineWidth指定线宽 MarkerEdgeColor指定标识符的边缘颜色 MarkerFaceColor指定标识符填充颜色 MarkerSize指定标识符的大小 注意上面四个属性是针对当前坐标系中所有曲线的 实例 % by dynamic % 2009.8.20 % X=1:10; % 两个都是数组,必须具有相同的尺寸 X1= ';%103 Y1=rand(10,3)+1;%103 % 其中一个为向量,另一个为数组,自动匹配尺寸相等方向 X2=1:0.1:10;%191 Y2= ';%912 % 其中一个是标量,另一为矢量,绘制垂直坐标轴的离散点 X3=1:10; Y3=-0.5; fh=figure('numbertitle','off','name','PLOT Usability Demo');%创建figure对象 ah=axes;%创建axes对象 h=plot(...%返回所有曲线句柄 ah,...%指定坐标系,可以省略,此时默认gca X1,Y1,...%坐标数据 '-.^',...%曲线属性,可以省略或部分省略,此时自动选择 X2,Y2,... 'm-',... X3,Y3,... 'o',...%注意此组数据设置线型和颜色无效,因为默认绘制离散点 'LineWidth',2,...%线宽 'MarkerEdgeColor','k',...%标识符边缘颜色 'MarkerFaceColor','r',...%标识符填充颜色 'MarkerSize',8)%标识符大小 set (gca,'XTick',-pi:pi/2:pi) set (gca,' XTickLabel ',{'-pi','-pi/2','0','pi/2','pi'})
个人分类: 未分类|9848 次阅读|0 个评论
关于画图软件Supermongo的一点总结
热度 1 qianlivan 2010-5-6 16:20
Supermongo是一些天文工作者常用的画图工具,曾经可以免费获取,但现在好像要收费了。我原来下载过一个版本,现在也不时用一下。 Supermongo的安装可以参考安装文件夹里的sm.install文件。基本就是 $ set_opts 这一步要选一些参数,设定安装路径(包括可被C和Fortran调用的库的路径)。 $ make $ make install 至此就安装完毕了。具体的使用没有太多可说的,用C和Fortran调用Supermongo的库有时候对计算有很大帮助,这样一来就可以在计算的同时可视化了,也就是可以监控计算的过程,出错了可以及时发现。我经常求解微分方程,这个功能对我比较有用。 下面就放一段程序作为示例。编译的时候要做如下操作(路径根据库的位置更改) $ g77 -c bs.f $ g77 -o g77 bs.o -L/home/lqian/local/lib/ -lplotsub -ldevices -lutils -L/usr/X11R6/lib/ -lX11 -o bs 然后就生成了可执行文件bs。 PROGRAM bs INTEGER nvar,nok,nbad REAL*8 EPS PARAMETER (EPS=1.d-12,nvar=4) REAL*8 hmin,h1,x1,x2,ystart(nvar) REAL smx,smmu,smnu,smsigma,smsigmap REAL*8 Omega,Lambda EXTERNAL bsstep,drv COMMON/para/ Omega, Lambda Omega=2.0d0 Lambda=0.0d0 cccccccccccc for plot on screen: if(sm_device(x11 -bg 0).lt.0)then cccccccccccc for plot in file test.ps: c if(sm_device(postfile :SY@: :OF@: equip.ps ). c 1 lt.0)then print*,'Cannot open required device' stop endif call sm_graphics call sm_defvar(TeX_strings,1) call sm_expand(1.2) call sm_limits(0.0,100.0,-0.1,0.1) ! set x,y axises call sm_location(3500,32000,3500,32000) call sm_box(1,2,0,0) call sm_xlabel('x') call sm_ylabel('...') call sm_alpha mu0=0.0d0 nu0=0.0d0 sigma0=0.1d0 sigmap0=0.0d0 x1=0.0D0 x2=0.01D0 ystart(1)=mu0 ystart(2)=nu0 ystart(3)=sigma0 ystart(4)=sigmap0 hmin=0.0D0 OPEN(1,FILE='bs.txt') DO 5 WHILE(x2.LT.1000.0D0) h1=(x2-x1)/10000.0 CALL odeint(ystart,nvar,x1,x2,EPS,h1,hmin,nok,nbad,drv,bsstep) PRINT *,x2,ystart(1),ystart(2),ystart(3),ystart(4) c WRITE(1,*)ystart(1),ystart(2) smx=x2 smmu=ystart(1) smnu=ystart(2) smsigma=ystart(3) smsigmap=ystart(4) call sm_ctype('blue') call sm_ptype(41.,0) call sm_points(smx,smmu,1) call sm_ctype('red') call sm_ptype(41.,0) call sm_points(smx,smnu,1) call sm_ctype('yellow') call sm_ptype(41.,0) call sm_points(smx,smsigma,1) call sm_ctype('green') call sm_ptype(41.,0) call sm_points(smx,smsigmap,1) x1=x2 x2=x1+0.01 5 CONTINUE CLOSE(1) PAUSE END SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs, *rkqs) INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX REAL*8 eps,h1,hmin,x1,x2,ystart(nvar),TINY EXTERNAL derivs,rkqs PARAMETER (MAXSTP=100000,NMAX=50,KMAXX=200,TINY=1.e-30) INTEGER i,kmax,kount,nstp REAL*8 dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),xp(KMAXX),y(NMAX), *yp(NMAX,KMAXX),yscal(NMAX) COMMON /path/ kmax,kount,dxsav,xp,yp x=x1 h=sign(h1,x2-x1) nok=0 nbad=0 kount=0 do 11 i=1,nvar y(i)=ystart(i) 11 continue if (kmax.gt.0) xsav=x-2.*dxsav do 16 nstp=1,MAXSTP call derivs(x,y,dydx) do 12 i=1,nvar yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY 12 continue if(kmax.gt.0)then if(abs(x-xsav).gt.abs(dxsav)) then if(kount.lt.kmax-1)then kount=kount+1 xp(kount)=x do 13 i=1,nvar yp(i,kount)=y(i) 13 continue xsav=x endif endif endif if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs) if(hdid.eq.h)then nok=nok+1 else nbad=nbad+1 endif if((x-x2)*(x2-x1).ge.0.)then do 14 i=1,nvar ystart(i)=y(i) 14 continue if(kmax.ne.0)then kount=kount+1 xp(kount)=x do 15 i=1,nvar yp(i,kount)=y(i) 15 continue endif return endif if(abs(hnext).lt.hmin) pause *'stepsize smaller than minimum in odeint' h=hnext 16 continue pause 'too many steps in odeint' return END SUBROUTINE bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext,derivs) INTEGER nv,NMAX,KMAXX,IMAX REAL*8 eps,hdid,hnext,htry,x,dydx(nv),y(nv),yscal(nv),SAFE1,SAFE2, *REDMAX,REDMIN,TINY,SCALMX PARAMETER (NMAX=50,KMAXX=8,IMAX=KMAXX+1,SAFE1=.25D0,SAFE2=.7D0, *REDMAX=1.d-5,REDMIN=.7D0,TINY=1.d-30,SCALMX=.1D0) CU USES derivs,mmid,pzextr INTEGER i,iq,k,kk,km,kmax,kopt,nseq(IMAX) REAL*8 eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest,xnew, *a(IMAX),alf(KMAXX,KMAXX),err(KMAXX),yerr(NMAX),ysav(NMAX), *yseq(NMAX) LOGICAL first,reduct SAVE a,alf,epsold,first,kmax,kopt,nseq,xnew EXTERNAL derivs DATA first/.true./,epsold/-1.D0/ DATA nseq /2,4,6,8,10,12,14,16,18/ if(eps.ne.epsold)then hnext=-1.d29 xnew=-1.d29 eps1=SAFE1*eps a(1)=nseq(1)+1 do 11 k=1,KMAXX a(k+1)=a(k)+nseq(k+1) 11 continue do 13 iq=2,KMAXX do 12 k=1,iq-1 alf(k,iq)=eps1**((a(k+1)-a(iq+1))/((a(iq+1)-a(1)+1.)*(2*k+ *1))) 12 continue 13 continue epsold=eps do 14 kopt=2,KMAXX-1 if(a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt))goto 1 14 continue 1 kmax=kopt endif h=htry do 15 i=1,nv ysav(i)=y(i) 15 continue if(h.ne.hnext.or.x.ne.xnew)then first=.true. kopt=kmax endif reduct=.false. 2 do 17 k=1,kmax xnew=x+h if(xnew.eq.x)pause 'step size underflow in bsstep' call mmid(ysav,dydx,nv,x,h,nseq(k),yseq,derivs) xest=(h/nseq(k))**2 call pzextr(k,xest,yseq,y,yerr,nv) if(k.ne.1)then errmax=TINY do 16 i=1,nv errmax=max(errmax,abs(yerr(i)/yscal(i))) 16 continue errmax=errmax/eps km=k-1 err(km)=(errmax/SAFE1)**(1./(2*km+1)) endif if(k.ne.1.and.(k.ge.kopt-1.or.first))then if(errmax.lt.1.)goto 4 if(k.eq.kmax.or.k.eq.kopt+1)then red=SAFE2/err(km) goto 3 else if(k.eq.kopt)then if(alf(kopt-1,kopt).lt.err(km))then red=1./err(km) goto 3 endif else if(kopt.eq.kmax)then if(alf(km,kmax-1).lt.err(km))then red=alf(km,kmax-1)*SAFE2/err(km) goto 3 endif else if(alf(km,kopt).lt.err(km))then red=alf(km,kopt-1)/err(km) goto 3 endif endif 17 continue 3 red=min(red,REDMIN) red=max(red,REDMAX) h=h*red reduct=.true. goto 2 4 x=xnew hdid=h first=.false. wrkmin=1.d35 do 18 kk=1,km fact=max(err(kk),SCALMX) work=fact*a(kk+1) if(work.lt.wrkmin)then scale=fact wrkmin=work kopt=kk+1 endif 18 continue hnext=h/scale if(kopt.ge.k.and.kopt.ne.kmax.and..not.reduct)then fact=max(scale/alf(kopt-1,kopt),SCALMX) if(a(kopt+1)*fact.le.wrkmin)then hnext=h/fact kopt=kopt+1 endif endif return END SUBROUTINE mmid(y,dydx,nvar,xs,htot,nstep,yout,derivs) INTEGER nstep,nvar,NMAX REAL*8 htot,xs,dydx(nvar),y(nvar),yout(nvar) EXTERNAL derivs PARAMETER (NMAX=50) INTEGER i,n REAL*8 h,h2,swap,x,ym(NMAX),yn(NMAX) h=htot/nstep do 11 i=1,nvar ym(i)=y(i) yn(i)=y(i)+h*dydx(i) 11 continue x=xs+h call derivs(x,yn,yout) h2=2.*h do 13 n=2,nstep do 12 i=1,nvar swap=ym(i)+h2*yout(i) ym(i)=yn(i) yn(i)=swap 12 continue x=x+h call derivs(x,yn,yout) 13 continue do 14 i=1,nvar yout(i)=0.5*(ym(i)+yn(i)+h*yout(i)) 14 continue return END SUBROUTINE pzextr(iest,xest,yest,yz,dy,nv) INTEGER iest,nv,IMAX,NMAX REAL*8 xest,dy(nv),yest(nv),yz(nv) PARAMETER (IMAX=13,NMAX=50) INTEGER j,k1 REAL*8 delta,f1,f2,q,d(NMAX),qcol(NMAX,IMAX),x(IMAX) SAVE qcol,x x(iest)=xest do 11 j=1,nv dy(j)=yest(j) yz(j)=yest(j) 11 continue if(iest.eq.1) then do 12 j=1,nv qcol(j,1)=yest(j) 12 continue else do 13 j=1,nv d(j)=yest(j) 13 continue do 15 k1=1,iest-1 delta=1./(x(iest-k1)-xest) f1=xest*delta f2=x(iest-k1)*delta do 14 j=1,nv q=qcol(j,k1) qcol(j,k1)=dy(j) delta=d(j)-q dy(j)=f1*delta d(j)=f2*delta yz(j)=yz(j)+dy(j) 14 continue 15 continue do 16 j=1,nv qcol(j,iest)=dy(j) 16 continue endif return END SUBROUTINE drv(X,Y,DYDX) IMPLICIT NONE INTEGER NMAX PARAMETER (NMAX=50) REAL*8 X,DYDX(NMAX),Y(NMAX) REAL*8 mu,nu,sigma,sigmap REAL*8 dmudx,dnudx,dsigmadx,dsigmapdx REAL*8 Omega,Lambda COMMON/para/ Omega, Lambda mu=Y(1) nu=Y(2) sigma=Y(3) sigmap=Y(4) IF(DABS(X).LT.1.0d-8) THEN dmudx=(Omega**2*dexp(-nu)*sigma**2+sigma**2+Lambda/2.0*sigma**4) 1*x*dexp(mu) dmudx=(Omega**2*dexp(-nu)*sigma**2-sigma**2-Lambda/2.0*sigma**4) 1*x*dexp(mu) dsigmadx=sigmap dsigmapdx=-dexp(mu-nu)*Omega**2*sigma+dexp(mu)*(1+Lambda*sigma** 12)*sigma ELSE dmudx=(Omega**2*dexp(-nu)*sigma**2+dexp(-mu)*sigmap**2+sigma**2 1+Lambda/2*sigma**4)*x*dexp(mu)-(dexp(mu)-1)/x dnudx=(Omega**2*dexp(-nu)*sigma**2+dexp(-mu)*sigmap**2-sigma**2 1-Lambda/2*sigma**4)*x*dexp(mu)+(dexp(mu)-1)/x dsigmadx=sigmap dsigmapdx=((sigma**2+Lambda/2*sigma**4)*x*dexp(mu)-(dexp(mu)-1) 1/x-2/x)*sigmap-dexp(mu-nu)*Omega**2*sigma+dexp(mu)*(1+Lambda* 2sigma**2)*sigma; ENDIF DYDX(1)=dmudx DYDX(2)=dnudx DYDX(3)=dsigmadx DYDX(4)=dsigmapdx RETURN END 有耐性看完此文的人的福利: http://www.lamost.org/~yzhao/soft/plotsoft/sm-2.3.1.tar.gz
个人分类: 总结|8403 次阅读|2 个评论
R1.8 画图参数-pch
anny424 2010-1-28 14:48
查看pch各参数形状: plot(rep(1,10),ylim=c(-2,1.2),pch=1:10,cex=3,axes=F,xlab=,ylab=) text(rep(0.6,10),as.character(1:10)) points(rep(0,10),pch=11:20,cex=3) text(rep(-0.4,10),as.character(11:20)) points(rep(-0.8,5),pch=21:25,cex=3) text(rep(-1.2,5),as.character(21:25)) points(6:10, rep(-0.8,5),pch=c(*,?,X,x,),cex=3) text(6:10,rep(-1.2,5),c(*,?,X,x,)) 参考资料: 台北大学 林建甫 C.F. Jeff Lin, PhD. R语言课件及代码
个人分类: 统计学习笔记|13490 次阅读|1 个评论
R1.2 画坐标、线和矩形图
anny424 2009-8-27 22:57
plot(c(0, 250), c(0, 350), type = n, xlab=, ylab=, main = plot) #画坐标 x=c(0,250) y=rep(200,length(x)) lines(x,y) #画第一条横线 z1=c(150,250) x1=c(100,100) x2=c(200,200) lines(x1,z1) #画第二条竖线 lines(x2,z1) #画第三条竖线 rect(125,175,175,225,col=black) #画矩形 画热度图矩形是最基本的 #画坐标,6个时间点,10个功能分类 plot(c(0,600),c(0,1000),type = n,xlab=percentage,ylab=,main = Function catagory) #第一个时间点上各功能分类 rect(0,0,100,100,col=red) rect(0,100,100,200,col=green) rect(0,200,100,300,col=blue) rect(0,300,100,400,col=yellow) rect(0,400,100,500,col=red) rect(0,500,100,600,col=red) rect(0,600,100,700,col=red) rect(0,700,100,800,col=red) rect(0,800,100,900,col=red) rect(0,900,100,1000,col=red) #第一个功能分类在各时间点上的情况 rect(0,0,100,100,col=red) rect(100,0,200,100,col=red) rect(200,0,300,100,col=red) rect(300,0,400,100,col=red) rect(400,0,500,100,col=red) rect(500,0,600,100,col=red) 基本的元素可以完成了,下一步学习循环语句,用循环来完成
个人分类: 统计学习笔记|6693 次阅读|1 个评论
[原创ZHOU Feng]ODV中定制海岸线、水深的方法
opensesame 2008-12-18 17:07
在Ocean Data View (ODV) 中定制水深、海岸线的方法 by ZHOU feng SOED, 2nd Institute of Oceanography, Hangzhou E-Mail: zhoufeng.hz@gmail.com ODV (Ocean Data View)是海洋科研中用处比较大的免费的画图软件,成图速度快,兼有数据库功能。但是软件提供的海岸线、水深资料很难满足近海画图的需要。需要自己定制。以前费了很多次没弄明白怎么定制,今天突然开窍了,赶紧记录下来。 本方案根据ODV中的说明形成,并有渤海成功案例 准备材料: ODV3.4, Surfer (一)需要有较高精度的水深或者岸线(岸线即水深为0处,因此以下论述均只谈水深的处理) 从地形文件中提取,如Etopo2m.cdf; 从海图中提取, 适合GIS 例子: (1)从Etopo2m.cdf中提取渤海水深资料,并通过Surfer画图软件画出地形,如10、20、25、30、50m等; (2)把上述地形图套上岸线,岸线资料如 bohuangdonghai.bln,并把区域集中到渤海海域,如117-122.5E, 37-41N; (3)利用surfer提供的digitize的功能,逐点描出10m等深线,由此可获得10m等深线的经度、纬度的bln文件, 姑且称为10m.bln, 等深线最好封闭,便于画阴影图,记录个数不要超过1500个(ODV的要求); (二)把水深文件转换成为.cdt (4)把10m.bln文件中的经度、纬度记录位置互换, ODV要求纬度在前、经度在后; 27 27 38.748741944 122.449183275 38.7104663748 122.3726331 38.7040869527 122.289703503 (5)10m.bln结尾加上ODV认可的结束符号0 0,并改名为10m.coa, 如下: 27 27 38.748741944 122.449183275 38.7104663748 122.3726331 38.7040869527 122.289703503 0 0 (6)转换成.cdt, 把10m.coa拷贝至ODV安装目录中的bin_w32\(其中含coa2cdt.exe转换程序,如果系统能找到该程序路径,则无所谓那个目录) (7)建立coa2cdt.inp文件,内容为需要转换呈.cdt的所有.coa文件,每个文件占一行,同时不加.coa的后缀,如: 10m (8)确保10m.coa, coa2cdt.inp和coa2cdt.exe同一目录, 然后执行coa2cdt.exe(windows下双击即可),即生成10m.cdt (9)在ODV的安装目录下的coast\目录下,分别建立目录并拷贝文件: coast\world.cdt (这个文件是岸线文件,ODV中自带GlobHR\下, 如果已制作也可以换) bathymetry\10m.cdt overlays\ topography\10m.cdt (10)然后在ODV的Map- right click-display options- layers-Ocean bathymetry-能看到定制的10m, 修改线条、阴影属性,然后选择, ok! 就成功了 (11)例子结果如图(太大,格式有限制,传不上) (卫星海洋环境动力学国家重点实验室, 杭州市保俶北路36号)
个人分类: 软件技巧|13426 次阅读|0 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-5-11 19:16

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部