GMT学习记录

最新推荐文章于 2024-07-15 12:01:59 发布

weixin_30663471

最新推荐文章于 2024-07-15 12:01:59 发布

阅读量638

点赞数 1

原文链接：http://www.cnblogs.com/eroszhou/p/4855721.html

版权

缘起
固体力学老师留了一个作业，用gmt画出洋壳俯冲时地幔楔的应力分布。相应的公式Donald L. Turcotte和Gerald Schubert写的Geodynamics中都有，要完成作业重点在使用gmt。
开始
1. 下载安装gmt5
  gmt是开源软件，主页在http://gmt.soest.hawaii.edu，里面有软件包下载，安装指南，以及gmt的使用手册。因为官方强烈建议使用gmt5，所以就下了gmt5。按照安装指南在我的ubuntu15.04上安装，没遇到问题，就是不知道怎么安装手册。
2. 熟悉画图原理
  官网上有cookbook，看了一会，有点烦，了解了一下工作方式。大概是给了一个模板，用户通过调整参数使得模板变成自己期待的样子，我用origin的时候也是存了好几个模板的。不过这显然不是重点，重点是如何表达数据。
  
  现在不想看cookbook，从其他地方找了个教程http://web.ics.purdue.edu/~ecalais/teaching/gmt/gmt.html，刚开始看（2015.10.05）。
  
  几张PPT看了一下，脚本用的是csh，与bash还是有点区别，gmt也不是我安装的版本，有些测试不能通过，比如彩色一直调不出来，还是回来看cookbook。
  1. 坐标纸设置

未完待续...

看得太烦了，暂时不想记录学习过程中开心或者不开心的事情了，赶在deadline之前把作业做了，脚本写得很乱。附上脚本和结果中的一张图吧，

#!/bin/bash
# Subducted Oceanic Crust
#
# Purpose: Stress distrubution
# GMT progs: grdmath, grdvector, pstext
#
theta=60 # The deg of subduction
velocity=5000 # The velocity of subduction
mu=10000000 # The viscosity of mantle
pi=3.14159265359 # The ratio of a circle's circumference to its diameter
unitv=3.1688955410 # The unit of velocity
g=9.80665 # The gravity acceleration speed
ru=3300000 # The density of mantle
XMAX=300; YMAX=300 # The deep we consider
XTICK=$(($XMAX/10)); YTICK=$(($YMAX/10))
YH=`echo "scale=1; $YMAX/$XMAX*6" | bc -l`i
SH=`echo "scale=1; $YMAX/$XMAX*3" | bc -l`i
GAP=15.0 # The gap of X-Y mesh

C=`echo "scale=20; $velocity*$unitv*s($theta/180*$pi)*$theta/180*$pi/(e(l($theta/180*$pi)*2.0)-e(l(s($theta/180*$pi))*2.0))" | bc -l`
D=`echo "scale=20; $velocity*$unitv*(s($theta/180*$pi)-$theta/180*$pi*c($theta/180*$pi))/(e(l($theta/180*$pi)*2.0)-e(l(s($theta/180*$pi))*2.0))" | bc -l`
#echo $C; echo $D
ps=homework.ps
gmt grdmath -R0/$XMAX/0/$YMAX -I$GAP $ru $g Y MUL MUL 2 $mu $C X MUL $D Y MUL ADD MUL MUL X Y R2 DIV SUB 2 $mu $D X MUL $C Y MUL SUB ABS MUL MUL X Y R2 DIV ADD 1000000000 DIV \
= sigma1.nc
gmt grdmath -R0/$XMAX/0/$YMAX -I$GAP $ru $g Y MUL MUL 2 $mu $C X MUL $D Y MUL ADD MUL MUL X Y R2 DIV SUB 2 $mu $D X MUL $C Y MUL SUB ABS MUL MUL X Y R2 DIV SUB 1000000000 DIV \
= sigma2.nc
gmt grdmath -R0/$XMAX/0/$YMAX -I1.0 2 $mu $D X MUL $C Y MUL SUB ABS MUL MUL X Y R2 DIV 1000000 DIV $theta 180 DIV PI MUL TAN Y X DIV GT MUL = tau.nc
gmt grdmath -R0/$XMAX/0/$YMAX -I$GAP X SQR Y SQR SUB 2 X Y MUL MUL DIV ATAN 2 DIV $D X MUL $C Y MUL GT PI 2 DIV MUL ADD = alpha.nc
#gmt grdmath -R0/$XMAX/0/$YMAX -I$GAP X SQR Y SQR SUB 2 X Y MUL MUL DIV ATAN 2 DIV = alpha.nc
gmt grdmath sigma1.nc alpha.nc COS MUL $theta 180 DIV PI MUL TAN Y X DIV GT MUL = sigma1_x.nc
gmt grdmath sigma1.nc alpha.nc SIN MUL $theta 180 DIV PI MUL TAN Y X DIV GT MUL = sigma1_y.nc
gmt grdmath sigma1_x.nc NEG = sigma1_x_.nc
gmt grdmath sigma1_y.nc NEG = sigma1_y_.nc
gmt grdmath sigma2.nc alpha.nc PI 2 DIV ADD COS MUL $theta 180 DIV PI MUL TAN Y X DIV GT MUL = sigma2_x.nc
gmt grdmath sigma2.nc alpha.nc PI 2 DIV ADD SIN MUL $theta 180 DIV PI MUL TAN Y X DIV GT MUL = sigma2_y.nc
gmt grdmath sigma2_x.nc NEG = sigma2_x_.nc
gmt grdmath sigma2_y.nc NEG = sigma2_y_.nc

makecpt -T-0/10000/1 -Z -Crainbow > tau.cpt
gmt grdimage tau.nc -R0/$XMAX/0/$YMAX -JX-6i/-$YH -nc -P -K -Bxa$XTICK+lkm -Bya$YTICK+lkm -BWSne -Ctau.cpt > $ps
psscale -D6.5i/$SH/$YH/0.15i -Bxa1000 -By+lMPa -Ctau.cpt -O -K >> $ps
psscale -D5.0i/0.5i/0.2i/0.5ph -Bxf10 -By+l10GPa -M -Ctau.cpt -O -K >> $ps

gmt grdvector sigma1_x.nc sigma1_y.nc -R -I$GAP -J -K -O -BWSne -W0.5p -S50i >> $ps
gmt grdvector sigma1_x_.nc sigma1_y_.nc -R -I$GAP -J -K -O -W -S50i >> $ps
gmt grdvector sigma2_x.nc sigma2_y.nc -R -I$GAP -J -K -O -W -S50i >> $ps
gmt grdvector sigma2_x_.nc sigma2_y_.nc -R -I$GAP -J -K -O -W -S50i >> $ps

TITLE="The primary stress and maximum shear stress | U0 = $velocity cm/year @%12%\\161@%% = $theta"
echo "150 -10 $TITLE" | gmt pstext -R -J -O -U -N -F+f16p,Times-Roman+jCB >> $ps
rm -f *.nc *.cpt